Modeling in R: lm(), glm() and Formula Syntax Step by Step

Who this is for: researchers and analysts who need to fit a regression-type model in R with the built-in lm() and glm() functions and read the result. After this page you can write a model formula, fit a linear or logistic model, extract the pieces you need, and compare two models.

Short answer: in R you describe a model with a formula (response ~ terms), pass it to a fitting function such as lm() for linear models or glm() for generalized linear models, store the result in an object, and then ask that object questions with extractor functions such as summary(), coef(), anova() and predict(). The worked output below was produced in R 4.6.1.

What it is used for in research

Model fitting in R covers the everyday inferential workload of empirical research: estimating a group difference, adjusting an effect for covariates, relating a binary outcome (for example, event or no event) to predictors, and testing whether an added term improves a model. The same formula syntax is shared across lm(), glm() and many add-on modeling functions, so learning it once carries over. The R manual An Introduction to R treats these under “Statistical models in R” and describes glm() families as a combination of a response distribution, a link function and other information (read 2026-11-08, manual for R 4.6.2).

Install or access

R is distributed through the Comprehensive R Archive Network at cran.r-project.org, which also hosts the official manuals cited on this page. lm() and glm() live in the stats package, which is part of base R and attached at startup, so no extra package is needed for anything in this guide.

Step by step: fit, extract and compare models with lm() and glm()

Step 1. Write the formula

The R reference for lm says a typical model has the form response ~ terms, and that a formula has an implied intercept. Remove the intercept with y ~ x - 1 or y ~ 0 + x. The operators below are from An Introduction to R, chapter 11.

Formula piece	Meaning (R manual wording, condensed)
Y ~ M	Y is modeled as M
M1 + M2	Include both M1 and M2
M1 - M2	Include M1 but leave out M2
M1 : M2	The tensor product (interaction) of M1 and M2
M1 * M2	Shorthand for M1 + M2 + M1:M2
M1 %in% M2	Similar to : but with a different coding
M1 / M2	Shorthand for M1 + M2 %in% M1 (nesting)
M^n	All terms in M together with interactions up to order n
I(M)	Insulates M so that ordinary arithmetic applies inside it, such as I(wt^2)

Step 2. Fit a linear model with lm()

This is the group-comparison example from the lm help page (Dobson, 1990): ten control and ten treatment plant weights.

ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
group <- gl(2, 10, 20, labels = c("Ctl","Trt"))
weight <- c(ctl, trt)
lm.D9  <- lm(weight ~ group)
lm.D90 <- lm(weight ~ group - 1)
summary(lm.D9)

Data can also be supplied through the data argument (a data frame, list or environment); otherwise variables are taken from the environment of the formula. Optional arguments include subset, weights (a non-NULL value gives weighted least squares, otherwise ordinary least squares) and na.action.

Step 3. Fit a generalized linear model with glm()

For a binary outcome, give glm() a family. In the manual, binomial uses the logit link by default and binomial(link = probit) switches to probit. The example uses R’s built-in mtcars data to model whether a car has a manual transmission (am) from its weight (wt, in 1000 lb).

fit <- glm(am ~ wt, family = binomial, data = mtcars)
summary(fit)
exp(coef(fit))

Step 4. Extract what you need

The manual lists these extractor functions for fitted model objects.

Function	Returns	Related CASRAI guide
coef(object)	Coefficients	How to read a regression table
summary(object)	A comprehensive regression summary	How to read a regression table
anova(object1, object2)	Comparison of a submodel with an outer model	Likelihood ratio test for nested models
predict(object, newdata = ...)	Predictions for new data	predict() in R
residuals(object)	Residuals	Assumptions of simple regression
deviance(object)	Residual sum of squares (weighted if appropriate)
vcov(object)	Variance-covariance matrix of the coefficients
step(object)	Stepwise add/drop; the returned model has the smallest AIC found in that search	Stepwise regression steps

The manual’s index also lists plot, update, add1 and drop1.

Step 5. Compare two models

Fit a smaller and a larger model on the same data and pass both to anova():

fit2 <- lm(mpg ~ wt + hp, data = mtcars)
fit3 <- lm(mpg ~ wt * hp, data = mtcars)
anova(fit2, fit3)

Output and how to read it

lm(): intercept coding versus no-intercept coding

Output of summary(lm.D9) (R 4.6.1):

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   5.0320     0.2202  22.850 9.55e-15 ***
groupTrt     -0.3710     0.3114  -1.191    0.249

Residual standard error: 0.6964 on 18 degrees of freedom
Multiple R-squared:  0.07308,	Adjusted R-squared:  0.02158
F-statistic: 1.419 on 1 and 18 DF,  p-value: 0.249

With the default intercept, the (Intercept) row is the mean of the reference level (Ctl, 5.032) and groupTrt is the difference of Trt from Ctl (-0.371). That difference is the quantity tested, with p = 0.249. Refitting with group - 1 removes the intercept and the two coefficients become the group means themselves: coef(lm.D90) returns groupCtl 5.032 and groupTrt 4.661. Same fit, different parametrisation; the t test in the second version tests each mean against zero, which is rarely the research question. anova(lm.D9) gives the same test as an F table: F = 1.4191 on 1 and 18 degrees of freedom, p = 0.249.

glm(): coefficients are on the link scale

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   12.040      4.510   2.670  0.00759 **
wt            -4.024      1.436  -2.801  0.00509 **

(Dispersion parameter for binomial family taken to be 1)
    Null deviance: 43.230  on 31  degrees of freedom
Residual deviance: 19.176  on 30  degrees of freedom
AIC: 23.176

For a binomial model with the default logit link, the wt coefficient of -4.024 is a change in log-odds per 1000 lb. Exponentiating gives an odds ratio of 0.0179 (exp(-4.024)). The summary reports z values rather than t because the dispersion is fixed at 1 for the binomial family. For predictions on the probability scale rather than the link scale, use predict(fit, newdata = data.frame(wt = 3), type = "response"), which returns 0.4921 here. See logistic regression in R for the full workflow.

anova(): does the interaction earn its place?

Model 1: mpg ~ wt + hp
Model 2: mpg ~ wt * hp
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)
1     29 195.05
2     28 129.76  1    65.286 14.088 0.0008108 ***

wt * hp expands to wt + hp + wt:hp, so the two models differ by exactly one term and the F test (1 and 28 degrees of freedom) addresses that term alone. Note that once an interaction is in the model, the wt and hp coefficients are no longer “main effects” at average values; they are the slopes when the other variable is zero (49.808 intercept, wt -8.217, hp -0.120, wt:hp 0.028 in this fit). For squared terms, wrap the arithmetic in I(): lm(mpg ~ wt + I(wt^2), data = mtcars) gives coefficients 49.9308, -13.3803 and 1.1711, whereas wt + wt^2 without I() would collapse to just wt because ^ is a formula operator.

Troubleshooting

• A coefficient name you did not expect (groupTrt). For a factor, R reports the non-reference levels as differences from the reference level. Check the reference with levels().
• Squared or arithmetic terms seem to vanish. ^, *, / and - are formula operators inside a formula. Wrap arithmetic in I().
• Rows silently drop. Missing values are handled by na.action, taken from R’s options and, if unset, na.fail. Compare nobs(fit) with nrow(data) before reporting a sample size.
• Intercept-free fit gives a huge R-squared. With - 1 or 0 +, R-squared is computed against zero rather than the mean, so it is not comparable to the with-intercept value.
• Checking assumptions afterwards. Run plot(fit) for the standard diagnostics, and see variance inflation factors in R and the Breusch-Pagan test in R for collinearity and heteroscedasticity checks.

Other software that fits the same models

The formula-and-extractor workflow here has counterparts in other tools; these CASRAI pages cover them on the site.

Software	Where to read about it on CASRAI
SPSS	SPSS vs R
Python (statsmodels)	statsmodels API; ANOVA in Python with statsmodels

For models that go beyond lm() and glm(), see generalised additive models in R, linear mixed-effects models in R and the overview of the generalized linear model.

FAQ

What is the difference between lm() and glm() in R?

lm() fits linear models by least squares (ordinary, or weighted when weights is supplied). glm() fits generalized linear models, where a family specifies a response distribution and a link function; binomial defaults to the logit link.

How do I remove the intercept in an R model?

Write y ~ x - 1 or y ~ 0 + x. A formula otherwise has an implied intercept.

What does * mean in an R formula?

M1 * M2 is shorthand for M1 + M2 + M1:M2: both terms plus their interaction.

How do I use a squared term in an R model?

Use I(x^2), because I() insulates its argument so arithmetic keeps its ordinary meaning inside the formula.

How do I choose between two nested models?

Pass both to anova(smaller, larger); the manual describes this as comparing a submodel with an outer model. step() automates add/drop and returns the model with the smallest AIC found in its search.

Sources

• R Core Team, An Introduction to R, chapter 11 “Statistical models in R” (manual for R 4.6.2), read 2026-11-08.
• R documentation, stats::lm reference, read 2026-11-08.
• All printed output on this page was generated by running the shown code in R 4.6.1 (2026-06-24) on 2026-11-08.