Linear models in R are rarely enough on their own.
You fit an lm(), look at the output, and realize quickly that your data doesn't cooperate with straight lines and constant variance. This is where extending the linear model becomes a daily practice rather than an academic exercise. R gives you several tools for this, and they are straightforward once you stop treating them like magic. The most common extension is glm() for generalized linear models. You keep the same formula syntax you already know from lm(), but you add a family argument that tells R what distribution your response variable follows and what link function to use. A binary outcome uses binomial with logit link, a count variable uses poisson with log link, and a strictly positive continuous variable with right-skewed errors might use Gamma with inverse link. The syntax change is minimal, but the interpretation shifts entirely because you are no longer modeling the raw mean of y. Predictions from a glm use binomial type="response" to get probabilities back on the original scale. Without that argument, you get logits, which are useless in a report.
I worked on a project last year where we modeled accident counts at traffic intersections using poisson regression. The residual plot showed clear overdispersion, which means the variance was substantially larger than the mean. The poisson model was technically mis-specified, and the standard errors were too small, making everything look more significant than it actually was. I switched to glm.nb() from the MASS package, which fits a negative binomial model with a dispersion parameter. The coefficient estimates barely changed, but the confidence intervals widened to something reasonable. This is one of those cases where the model looks fine until you check the assumptions properly, and the fix is almost always one line of code.
Polynomial and interaction terms
Before reaching for a generalized model, you should check whether simple transformations of your predictors solve the problem. Adding I(x^2) or using poly(x, degree=2) in your formula lets lm() fit curvature without any extra packages. Interaction terms use the * and : operators. x*z includes main effects and the interaction, while x:y includes only the interaction. This is obvious to anyone who has read a textbook, but people still write x*z when they only want the interaction term, and then wonder why their model looks inflated. A practical note about poly(): orthogonal polynomials are numerically safer than raw powers, especially when degree goes above 2. The coefficients are harder to interpret directly, but the model fitting is more stable and the Anova() output gives clean sequential tests.
Mixed effects models
When your data has a hierarchical structure, like students nested inside schools or repeated measurements on the same subject, ordinary lm() breaks the independence assumption. The lme4 package handles this with lmer() for Gaussian responses and glmer() for non-Gaussian responses. You specify random intercepts with (1|group) and random slopes with (variable|group). The syntax is compact, but getting the model to converge is not guaranteed. I ran into a boundary convergence warning recently on a glmer model with a binary outcome and a random slope for a continuous predictor across about forty groups. The variance component for that slope was estimated near zero, which triggered the warning. I simplified the random effects structure to a random intercept only, and the warning disappeared. The fixed effects remained essentially unchanged. If you encounter this, check the correlation between random effects with VarCorr() before dropping terms blindly, because sometimes the issue is a near-perfect correlation rather than a zero variance.
Nonlinear models
nls() fits nonlinear least squares models. You write the actual functional form of your model, not a linear approximation. A classic example is the Michaelis-Menten enzyme kinetics equation: rate ~ Vmax*substrate/(Km + substrate). You need starting values for Vmax and Km, and if your guesses are far off, nls() will fail or converge to a local minimum. The SSmicmen self-starting function can generate reasonable initial values automatically, which saves a lot of fiddling. I once spent three hours debugging an nls() model that appeared to converge but produced nonsensical parameter estimates. The issue was that my starting values were close enough for the optimizer to find a local minimum on a poorly scaled surface. Centering and scaling the predictor variable before fitting resolved it completely. The model converged in two seconds after that.
Generalized additive models
gam() from the mgcv package extends the linear framework by allowing smooth, nonparametric functions of predictors. You write s(x) in the formula instead of x, and R estimates the shape of the relationship from the data. This is powerful when you suspect a nonlinear pattern but have no theory about its exact form. The gam.check() function evaluates whether your smoothers are overfitting or underfitting using residual diagnostics and effective degrees of freedom. A limitation many people overlook is that GAMs can be slow with large datasets. Fitting a model with ten smooth terms on two million rows might take twenty minutes or more, while an equivalent glm runs in under a second. If your sample is that large, consider sampling a subset for model development and then refitting on the full data, or use the bam() function in mgcv, which is designed for big data and uses P-splines with integrated nested Laplace approximation.
Model comparison and selection
AIC and BIC from the AIC() function are the standard tools for comparing non-nested models. For nested models, anova() performs a likelihood ratio test. Neither approach is perfect, and both will happily suggest complex models that overfit your data. I have seen people run stepAIC() forward and backward and end up with a model that looks impressive on their training data but fails immediately on new observations. Regularization methods like glmnet, which fit lasso or elastic net penalties, are often more reliable for high-dimensional data where the number of predictors approaches or exceeds the sample size. Diagnostic plots for glm and gam objects work differently than for lm(). plot(glm_model) gives deviance residuals against fitted values, not ordinary residuals. For binary models, the deviance residuals are bounded and clustering patterns in the residual plot often point to missing predictors or incorrect link functions. The DHARMa package simulates standardized residuals that are uniformly distributed under the null hypothesis, making diagnosis much more intuitive. I use this for almost every glm and gam I fit because the default residual plots can be misleading, especially with small datasets or discrete outcomes. The key thing most tutorials skip is that extending the linear model is not about finding the best-fitting model. It is about specifying a model that matches the data generating process well enough to make reliable predictions or valid inferences. R makes it easy to add complexity. It does not make it easy to justify that complexity, and you should be uncomfortable when your extended model has more parameters than you can explain in a single sentence.
