Statistical Models
|
This content was generated with the assistance of AI and should be verified against the official documentation before being relied on in production. The reference material consulted while preparing these pages is listed in the section bibliography. |
Checked against R 4.6.1 (2026-06-24) on 2026-10-20; R 4.6.2 is scheduled for 2026-10-29. Almost everything on this page is base R: the stats package (lm(), glm(), nls(), optim() and the generics that read their results), stats4 (mle()), graphics and grDevices (the diagnostic plots, drawn to a device that is thrown away) and utils; all of them ship with R. One add-on package is run once, in a short section: broom 1.0.13. The recommended packages nlme 3.1-169, MASS 7.3-65, mgcv 1.9-4 and survival 3.8-6, which ship with R, are named but not run. The lintr 3.4.0 default linters were used only to check the scripts. A statistical model in R is an ordinary object: a function such as lm() takes a formula and a data frame, fits the model and returns a list with a class, and a small set of generic functions (print(), summary(), coef(), residuals(), predict(), anova(), confint(), plot(), update()) read that list so that you never have to. The same design carries from straight-line regression to generalised linear models, nonlinear least squares and hand-written maximum likelihood, which is why this page is organised around the fitting functions and one shared vocabulary of extractors. It covers the formula interface of lm() and the full summary.lm() output, every extractor generic and what each returns, missing values and na.action, weights, offsets and subsets, influence measures and the diagnostic plots, update() and stepwise selection, factor predictors and contrasts, interactions and polynomials, glm() with its families and links (a logistic model for overdue loans, counts, positive durations, grouped binomial data), nls() and its failures, optim(), optimize(), nlm(), nlminb() and stats4::mle() for maximum likelihood, the model object as an S3 list with a wrapper class of your own, and broom.
Every example is a script that was run with Rscript on R 4.6.1 with warnings printed as they happen (options(warn = 1)), and the lines that start with #> are what it printed. The scripts also pass the default lintr linters. The names books, loans, shelf_today, loan_days and the book() constructor come from the running library-shelf scenario described on the section index: a script that needs them starts with source("scenario.R"). Each example repeats the preparation it needs, so every listing runs on its own; the scripts that model loan durations all start with the same six statements, explained in The data for the examples, which build the data frame dur (one row per returned loan, with the book’s pages, year and genre). Where a message contains quotation marks, they are the typographic quotes R prints in a UTF-8 locale. A formula created at the top level of a script has the global environment as its environment, so printing it shows no <environment: ...> line; the one place where a model object holds an environment is its terms component, and the scripts never print it. Diagnostic plots are drawn to svg() or pdf(NULL) devices that are thrown away, so that the scripts run without a screen: the plots are not reproduced here, and the numbers behind each one are computed and printed instead.
This page follows Formulas, which explains the formula language itself (the operators +, -, :, *, ^, %in%, . and I(), terms, model frames and model matrices) and is not repeated here; the small amount of formula syntax these models need is explained where it appears. A fitted model is an S3 object, so S3 Classes and Generics explains how summary() and predict() find their methods, and Factors explains the genre predictor. The classical tests (t.test(), chisq.test(), aov() with TukeyHSD()) belong to Descriptive Statistics and Tests, and the distribution functions used by the likelihood examples (dnorm(), dgamma(), rnorm(), set.seed()) to Distributions and Random Numbers. The shelf data is small and invented: the models below have honest, mostly non-significant results, and the page says so wherever it matters. They are teaching examples of how to fit and read a model, not findings about libraries.
How R fits models
The chapter on statistical models in An Introduction to R (chapter 11) starts from one idea: a model is described by a formula, fitted by a function that takes the formula and a data frame, and the result is an object that other functions take apart. The manual’s own example is lm(y ~ x1 + x2, data = production), and it points out that the data argument makes the variables come from that data frame first, whether or not the frame is attached to the search path, which is the first reason not to use attach() (An Introduction to R: Linear models). The value of lm() is, technically, a list of class "lm", and generic functions that dispatch on that class display, extract and plot the results (An Introduction to R: Generic functions for extracting model information). That is the design this page follows: learn the fitting functions once, and learn the vocabulary of extractors once, because the same words (coef, residuals, predict, anova, logLik, AIC) work for most of them.
| Function | What it fits | Class of the result | Where in this page |
|---|---|---|---|
|
Linear models by least squares (QR decomposition) |
|
|
|
Generalised linear models by iteratively reweighted least squares |
|
|
|
Nonlinear least squares (Gauss-Newton, or the PORT library with bounds) |
|
|
|
General-purpose optimisation and root finding; no model object |
a plain list |
|
|
Maximum likelihood for a likelihood you write |
S4 class |
|
|
Analysis of variance, including multistratum designs with |
|
|
|
Local regression, smoothing splines, projection pursuit (all in |
their own classes |
|
|
Mixed models, robust regression, additive models, survival models (recommended packages shipped with R; named here, not run) |
their own classes |
Two further points of design are worth knowing before the code. First, the modelling functions evaluate the formula against the data frame with the same data-mask idea as with() and subset(), through model.frame() (see Evaluation and Non-Standard Evaluation): a variable not found in data is looked up in the formula’s environment, which is a source of confusing errors shown below. Second, the extractor generics are the interface and the list components are the implementation: coef(fit) and fit$coefficients return the same numbers for an lm today, but only the generic works for an nls fit, a model from another package, or a class that stores its coefficients under another name (see The model object as an S3 list). The R FAQ lists the differences from S that matter when reading old S-PLUS code: powers must be protected with I(), na.action is "na.omit" by default in R and unset in S, GLM family objects are implemented differently, terms objects are formulas with attributes, and unordered factors get treatment contrasts in R and Helmert contrasts in S (R FAQ 3.3.2: Models).
The data for the examples
A few examples use books alone (the number of pages on the publication year). The others describe how long a returned loan lasted, so it needs one row per loan with the book’s attributes next to it. That is a join of loans to books on the ISBN, and it needs a little cleaning first: the ISBNs in books.csv are hyphenated in some rows and the ones in loans.csv are bare digits, so both are normalised with normalise_isbn() from the scenario; the genre column holds nine spellings of six genres ("Fiction", "fiction", "Fiction "), so it is cleaned with tolower(trimws()), exactly as on the formulas page. The join is merge(loans, books, by = "isbn"), an inner join: a loan whose ISBN matches no book (loan L005 has a digit missing) drops out, and a book never loaned does not appear. The duration is the difference of two Date columns, a difftime, which as.numeric() turns into a plain number of days (see Dates and Times).
source("scenario.R")
# Clean the two columns that the join and the model use.
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
# One row per loan, with the book's pages, year and genre next to it.
loan_books <- merge(loans, books, by = "isbn")
c(loans = nrow(loans), books = nrow(books), merged = nrow(loan_books))
#> loans books merged
#> 60 16 59
# Loan duration in days: a difftime becomes a number.
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
class(loan_books$days)
#> [1] "numeric"
# Keep the returned loans with a sensible duration (this drops the NA
# dates, the unreturned loans and the loan returned before it began).
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
nrow(dur)
#> [1] 50
table(dur$genre)
#>
#> classic fiction mystery non-fiction poetry
#> 5 18 4 5 5
#> science fiction
#> 13
summary(dur$days)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 1.00 14.50 21.00 20.22 26.75 37.00
sum(complete.cases(dur[, c("days", "pages", "genre")]))
#> [1] 42Source: R documentation: merge
The numbers tell the story of the data. Sixty loans and sixteen books give 59 rows after the join (the loan with the short ISBN found no book). Of those, 50 are returned loans with a duration of zero or more days: seven loans are still out (the NA return dates), one has the impossible loan date 2025-02-30 (an NA date), and one was "returned" three days before it was loaned. Of the 50, only 42 have a known number of pages and a genre, because two books have NA pages. Every later script that models loan durations starts with these same six statements (the script above, with the comments left out), and filters further where it needs to.
Two warnings about this data apply to every model that follows. The sample is small (42 usable rows, and 13 distinct books among them), so every standard error is large and every test has little power; and the rows are not independent, because the same book appears in several loans and shares its pages and genre (the section on pitfalls shows what that does). The models fit and read correctly; the conclusions are weak, and the text says so.
Linear models with lm()
The formula interface
lm() takes a formula and a data frame. The call lm(pages ~ year, data = books) regresses the page count on the publication year: the response is on the left of the tilde, the predictors on the right, an intercept is added unless it is removed with -1 or +0, and the variables are found in data first. The value prints compactly: the call and the estimated coefficients. Everything else is hidden inside the object, and summary() returns a second object whose print method formats the full table (R documentation: lm).
source("scenario.R")
# The data argument says where the variables come from.
fit <- lm(pages ~ year, data = books)
fit
#>
#> Call:
#> lm(formula = pages ~ year, data = books)
#>
#> Coefficients:
#> (Intercept) year
#> 2067.0927 -0.8945
#>
class(fit)
#> [1] "lm"
typeof(fit)
#> [1] "list"
# summary() returns a new object whose print method does the formatting.
summary(fit)
#>
#> Call:
#> lm(formula = pages ~ year, data = books)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -185.617 -47.094 -1.516 56.433 147.117
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2067.0927 915.8707 2.257 0.0453 *
#> year -0.8945 0.4586 -1.950 0.0771 .
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Residual standard error: 89.9 on 11 degrees of freedom
#> (3 observations deleted due to missingness)
#> Multiple R-squared: 0.257, Adjusted R-squared: 0.1894
#> F-statistic: 3.804 on 1 and 11 DF, p-value: 0.07707
#>Source: R documentation: summary.lm
The summary() output has five parts, and each one answers a different question:
| Part | What it says |
|---|---|
|
The matched call, as |
|
The five-number summary of the residuals (minimum, quartiles, maximum) when there are more than five residual degrees of freedom; with fewer, the residuals themselves are printed. A residual is the observed value minus the fitted value. |
|
One row per estimated coefficient: the estimate, its standard error, the t value (estimate divided by standard error) and the two-sided p-value of the null hypothesis that the coefficient is zero, from a t distribution with the residual degrees of freedom. The stars are a legend ( |
|
The estimate of the error standard deviation, the square root of the residual sum of squares divided by the residual degrees of freedom; |
|
The share of the variance of the response explained by the model, the same share penalised for the number of predictors, and the F test of the model against the model with the intercept only. For a model without an intercept, R-squared is computed against zero rather than against the mean (R FAQ 7.41). |
Source: R documentation: summary.lm
Read the pages-on-year summary like this. Thirteen books have both a year and a page count (three were deleted for missing values, and the output says so), so there are 11 residual degrees of freedom: 13 observations minus two coefficients. The slope is about -0.89 pages per year with a standard error of 0.46, a t value of -1.95 and a p-value of 0.077: newer books in this tiny sample are a little shorter, and the evidence is weak. R-squared is 0.257, adjusted R-squared 0.189, and the F test of the whole model has the same p-value as the slope’s t test, because a model with one predictor has one test (F equals t squared: 3.804 is -1.950 squared).
The second model is the one most of the page uses: loan duration in days on the number of pages and the genre.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
# Loan duration on pages and genre; genre is character, so lm() makes a factor.
fit <- lm(days ~ pages + genre, data = dur)
fit
#>
#> Call:
#> lm(formula = days ~ pages + genre, data = dur)
#>
#> Coefficients:
#> (Intercept) pages genrefiction
#> 21.0543745 0.0009311 -3.7882244
#> genremystery genrenon-fiction genrepoetry
#> -0.5448319 0.2826825 -1.1407808
#> genrescience fiction
#> -0.8625698
#>
summary(fit)
#>
#> Call:
#> lm(formula = days ~ pages + genre, data = dur)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -16.451 -6.974 1.266 5.885 15.426
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 21.0543745 16.0118958 1.315 0.197
#> pages 0.0009311 0.0416789 0.022 0.982
#> genrefiction -3.7882244 8.2943031 -0.457 0.651
#> genremystery -0.5448319 7.7757046 -0.070 0.945
#> genrenon-fiction 0.2826825 11.2448822 0.025 0.980
#> genrepoetry -1.1407808 12.9878207 -0.088 0.931
#> genrescience fiction -0.8625698 4.8550186 -0.178 0.860
#>
#> Residual standard error: 9.226 on 35 degrees of freedom
#> (8 observations deleted due to missingness)
#> Multiple R-squared: 0.03251, Adjusted R-squared: -0.1333
#> F-statistic: 0.196 on 6 and 35 DF, p-value: 0.9758
#>Source: R documentation: lm
genre is a character column, so lm() turns it into a factor (the levels are sorted alphabetically, so "classic" comes first) and codes it with treatment contrasts: the intercept is the predicted duration for the reference level "classic" at zero pages, and each genre coefficient is the difference between that genre and the reference. The coefficient of pages is 0.0009 days per page with a p-value of 0.98, and every genre coefficient has a p-value above 0.6; the F test of the whole model has p = 0.976 and the adjusted R-squared is negative, which happens when a model explains less than a random predictor would by chance. Duration does not depend on pages or genre in this data, and the model says that correctly. A negative adjusted R-squared is not a bug; it is the adjustment working.
The extractor generics
The vocabulary below is the part to memorise. Each row is a generic function; the right-hand column says what it returns for an lm. The manual lists the same family (add1, alias, anova, coef, deviance, drop1, effects, family, formula, kappa, labels, plot, predict, print, proj, residuals, step, summary, vcov) in its chapter on extracting model information (An Introduction to R: Generic functions for extracting model information).
| Generic | Returns |
|---|---|
|
A named numeric vector of the estimated coefficients. |
|
The residuals, named by the row names of the data. |
|
The fitted values, named by row. |
|
Predictions for new rows of data, optionally with |
|
Confidence intervals for the coefficients (t based for an |
|
The estimated variance-covariance matrix of the coefficients. |
|
A sequential analysis of variance, or a comparison of nested models. |
|
The maximised log-likelihood (with its degrees of freedom), and the information criteria derived from it. |
|
The residual sum of squares (weighted if appropriate) for an |
|
Residual degrees of freedom, the number of observations used, the residual standard error. |
|
The formula, its terms object, the data frame the fit used and the design matrix. |
|
Leverage, Cook’s distance and the full set of influence diagnostics. |
|
Refit with a changed formula; test single additions or deletions; stepwise selection by AIC. |
|
Up to six diagnostic plots. |
Source: R documentation: coef
The next two scripts exercise them on the duration model. The first covers the numbers you read most often: the coefficients, residuals and fitted values, confidence intervals and the covariance matrix.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
fit <- lm(days ~ pages + genre, data = dur)
# The coefficients, as a named numeric vector.
coef(fit)
#> (Intercept) pages genrefiction
#> 21.0543744546 0.0009311033 -3.7882244396
#> genremystery genrenon-fiction genrepoetry
#> -0.5448318821 0.2826824677 -1.1407808409
#> genrescience fiction
#> -0.8625698375
round(coef(fit), 3)
#> (Intercept) pages genrefiction
#> 21.054 0.001 -3.788
#> genremystery genrenon-fiction genrepoetry
#> -0.545 0.283 -1.141
#> genrescience fiction
#> -0.863
# A coefficient matrix from the summary: estimate, standard error, t, p.
coef(summary(fit))["pages", ]
#> Estimate Std. Error t value Pr(>|t|)
#> 0.0009311033 0.0416789375 0.0223399001 0.9823037008
# Residuals and fitted values are named by the row names of the data.
head(residuals(fit), 4)
#> 1 2 3 4
#> -7.465406 1.534594 13.534594 7.534594
head(fitted(fit), 4)
#> 1 2 3 4
#> 17.46541 17.46541 17.46541 17.46541
head(resid(fit) + fitted(fit), 4)
#> 1 2 3 4
#> 10 19 31 25
head(dur$days, 4)
#> [1] 10 19 31 25
# Residuals sum to zero when there is an intercept.
all.equal(sum(residuals(fit)), 0)
#> [1] TRUE
# Confidence intervals for the coefficients (t-based).
round(confint(fit), 2)
#> 2.5 % 97.5 %
#> (Intercept) -11.45 53.56
#> pages -0.08 0.09
#> genrefiction -20.63 13.05
#> genremystery -16.33 15.24
#> genrenon-fiction -22.55 23.11
#> genrepoetry -27.51 25.23
#> genrescience fiction -10.72 8.99
round(confint(fit, "pages", level = 0.9), 4)
#> 5 % 95 %
#> pages -0.0695 0.0714
# The estimated covariance matrix of the coefficients.
round(vcov(fit)[1:3, 1:3], 3)
#> (Intercept) pages genrefiction
#> (Intercept) 256.381 -0.645 -121.067
#> pages -0.645 0.002 0.280
#> genrefiction -121.067 0.280 68.795
round(sqrt(diag(vcov(fit))), 4)
#> (Intercept) pages genrefiction
#> 16.0119 0.0417 8.2943
#> genremystery genrenon-fiction genrepoetry
#> 7.7757 11.2449 12.9878
#> genrescience fiction
#> 4.8550Source: R documentation: confint
The residuals and the fitted values add back up to the response (resid(fit) + fitted(fit) equals dur$days), and with an intercept in the model the residuals sum to zero. confint() gives 95 per cent intervals by default; for an lm they are the estimate plus and minus the t quantile times the standard error, and level = changes the width. Every interval for pages contains zero and the intervals for the genres are wide (the one for genrepoetry runs from -27.5 to 25.2 days): the data cannot distinguish any genre from another.
The second script covers the counts and the likelihood-based summaries. For a linear model deviance() is the residual sum of squares, and logLik() assumes normal errors with the variance estimated by maximum likelihood, which is why its degrees of freedom (8) are one more than the number of coefficients (7): the error variance counts. AIC() is minus twice the log-likelihood plus twice that number of parameters, and BIC() replaces the 2 by the logarithm of the number of observations (AIC(fit, k = log(nobs(fit))) is the BIC).
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
fit <- lm(days ~ pages + genre, data = dur)
# Counts and degrees of freedom.
nobs(fit)
#> [1] 42
df.residual(fit)
#> [1] 35
fit$rank
#> [1] 7
length(coef(fit))
#> [1] 7
# Residual sum of squares, and sigma.
deviance(fit)
#> [1] 2978.869
sum(residuals(fit)^2)
#> [1] 2978.869
sigma(fit)
#> [1] 9.225538
summary(fit)$sigma
#> [1] 9.225538
# Likelihood-based summaries.
logLik(fit)
#> 'log Lik.' -149.0896 (df=8)
AIC(fit)
#> [1] 314.1793
BIC(fit)
#> [1] 328.0806
AIC(fit, k = log(nobs(fit)))
#> [1] 328.0806
# The formula, the terms and the call.
formula(fit)
#> days ~ pages + genre
attr(terms(fit), "term.labels")
#> [1] "pages" "genre"
fit$call
#> lm(formula = days ~ pages + genre, data = dur)
fit$xlevels
#> $genre
#> [1] "classic" "fiction" "mystery" "non-fiction"
#> [5] "poetry" "science fiction"
#>Source: R documentation: AIC and BIC
The model frame and the design matrix
Behind every fit are two objects that formulas.adoc explains in detail: the model frame, a data frame holding exactly the variables the formula mentions after na.action has removed incomplete rows, and the model matrix, the numeric design matrix built from it (one column of ones for the intercept, the numeric predictors, and one 0/1 column per non-reference factor level). model.frame(fit) and model.matrix(fit) recover them from the fit, which is how you check what a model actually saw. The coefficients of a linear model are the least-squares solution of the normal equations for that matrix, which the last lines of the script verify (lm() itself uses a QR decomposition, which is numerically safer than forming the cross-product matrix; the manual describes it in the appendix on least-squares fitting, An Introduction to R: Least squares fitting and the QR decomposition).
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
fit <- lm(days ~ pages + genre, data = dur)
# The model frame: the variables the fit used, after na.action.
mf <- model.frame(fit)
dim(mf)
#> [1] 42 3
names(mf)
#> [1] "days" "pages" "genre"
head(mf, 3)
#> days pages genre
#> 1 10 214 fiction
#> 2 19 214 fiction
#> 3 31 214 fiction
# The design matrix: an intercept, a numeric column and treatment dummies.
mm <- model.matrix(fit)
dim(mm)
#> [1] 42 7
colnames(mm)
#> [1] "(Intercept)" "pages" "genrefiction"
#> [4] "genremystery" "genrenon-fiction" "genrepoetry"
#> [7] "genrescience fiction"
head(mm, 3)
#> (Intercept) pages genrefiction genremystery genrenon-fiction genrepoetry
#> 1 1 214 1 0 0 0
#> 2 1 214 1 0 0 0
#> 3 1 214 1 0 0 0
#> genrescience fiction
#> 1 0
#> 2 0
#> 3 0
# The coefficients are the least-squares solution of mm %*% b = y.
b <- solve(crossprod(mm), crossprod(mm, model.response(mf)))
all.equal(drop(b), coef(fit))
#> [1] TRUESource: R documentation: model.frame
The frame has 42 rows, not 50: the eight loans of books with unknown pages are gone, and names(mf) shows the three variables the formula used (days, pages, genre), not the whole dur.
Prediction
predict() evaluates the fitted model on new data. The data frame must contain every variable of the right-hand side under the same name, and a factor may only take levels that were in the fit. The interval argument adds either a confidence interval for the mean response (the uncertainty of the fitted line) or a prediction interval for a single new observation (that uncertainty plus the residual variance), and the prediction interval is always the wider one. se.fit = TRUE returns a list holding the fit, the standard errors, the residual degrees of freedom and the residual scale, and with no newdata the function returns the fitted values.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
fit <- lm(days ~ pages + genre, data = dur)
# New data needs the same variable names, and genre must be a known level.
new <- data.frame(pages = c(150, 300, 450), genre = "fiction")
predict(fit, newdata = new)
#> 1 2 3
#> 17.40582 17.54548 17.68515
# Confidence interval for the mean, prediction interval for a new loan.
predict(fit, newdata = new, interval = "confidence")
#> fit lwr upr
#> 1 17.40582 10.151590 24.66004
#> 2 17.54548 8.316660 26.77430
#> 3 17.68515 -3.288354 38.65865
predict(fit, newdata = new, interval = "prediction", level = 0.8)
#> fit lwr upr
#> 1 17.40582 4.4829530 30.32868
#> 2 17.54548 4.1113965 30.97957
#> 3 17.68515 -0.4069357 35.77723
# se.fit = TRUE returns a list: fit, se.fit, df and residual.scale.
p <- predict(fit, newdata = new, se.fit = TRUE)
names(p)
#> [1] "fit" "se.fit" "df" "residual.scale"
round(p$se.fit, 3)
#> 1 2 3
#> 3.573 4.546 10.331
p$df
#> [1] 35
# predict() with no newdata returns the fitted values.
all.equal(predict(fit), fitted(fit))
#> [1] TRUE
# A level the model has never seen is an error.
new2 <- data.frame(pages = 300, genre = "cookery")
tryCatch(predict(fit, newdata = new2), error = \(e) conditionMessage(e))
#> [1] "factor genre has new level cookery"
# A missing variable is an error too.
tryCatch(predict(fit, newdata = data.frame(pages = 300)),
error = \(e) conditionMessage(e)
)
#> [1] "object 'genre' not found"
# A model fitted with dur$pages in the formula has no variable called pages.
dollar_fit <- lm(dur$days ~ dur$pages)
tryCatch(
predict(dollar_fit, newdata = data.frame(pages = 300)),
warning = \(w) conditionMessage(w)
)
#> [1] "'newdata' had 1 row but variables found have 50 rows"Source: R documentation: predict.lm
Both errors at the end are the typical ones. factor genre has new level cookery says that a prediction cannot be made for a level the model never saw; and object 'genre' not found says the formula’s variable is missing from newdata. A related mistake is to write the data frame into the formula, lm(dur$days ~ dur$pages): the model then has a variable called dur$pages, not pages, and predict(fit, newdata = data.frame(pages = 300)) warns 'newdata' had 1 row but variables found have 50 rows and returns 50 predictions, one for each row of dur (NA for the 8 rows without a page count), instead of one prediction (the last lines of the script). Always use data =.
Comparing models with anova()
anova() on one lm is a sequential analysis of variance: each term’s sum of squares is the decrease in the residual sum of squares when that term is added to the terms before it, in the order they appear in the formula. The manual says the table is for a sequence of fitted models and that only for orthogonal experiments is the order of inclusion irrelevant (An Introduction to R: ANOVA tables). The R FAQ explains why the numbers change with the order of the terms, and what to use if you want the "Type III" sums of squares that other packages print: compare two models with anova(model1, model2), or call drop1(model), which shows the change from dropping each term from the full model (R FAQ 7.18: Why does the output from anova() depend on the order of factors in the model?).
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
# Nested models must be fitted to the same rows, so keep complete cases.
dur <- dur[complete.cases(dur[, c("days", "pages", "genre")]), ]
# anova() on one model: sequential (type I) sums of squares.
fit <- lm(days ~ pages + genre, data = dur)
anova(fit)
#> Analysis of Variance Table
#>
#> Response: days
#> Df Sum Sq Mean Sq F value Pr(>F)
#> pages 1 21.13 21.131 0.2483 0.6214
#> genre 5 78.98 15.795 0.1856 0.9662
#> Residuals 35 2978.87 85.111
# The order of the terms changes the table (R FAQ 7.18).
fit_rev <- lm(days ~ genre + pages, data = dur)
anova(fit_rev)
#> Analysis of Variance Table
#>
#> Response: days
#> Df Sum Sq Mean Sq F value Pr(>F)
#> genre 5 100.06 20.013 0.2351 0.9444
#> pages 1 0.04 0.042 0.0005 0.9823
#> Residuals 35 2978.87 85.111
# anova() on nested models compares them with an F test.
fit0 <- lm(days ~ 1, data = dur)
fit1 <- lm(days ~ pages, data = dur)
anova(fit0, fit1, fit)
#> Analysis of Variance Table
#>
#> Model 1: days ~ 1
#> Model 2: days ~ pages
#> Model 3: days ~ pages + genre
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 41 3079.0
#> 2 40 3057.8 1 21.131 0.2483 0.6214
#> 3 35 2978.9 5 78.976 0.1856 0.9662
# drop1() tests each term when it is dropped from the full model.
drop1(fit, test = "F")
#> Single term deletions
#>
#> Model:
#> days ~ pages + genre
#> Df Sum of Sq RSS AIC F value Pr(>F)
#> <none> 2978.9 192.99
#> pages 1 0.042 2978.9 190.99 0.0005 0.9823
#> genre 5 78.976 3057.8 184.09 0.1856 0.9662
# The models must use the same rows: this fit loses different observations.
fit_year <- lm(days ~ year, data = dur)
nobs(fit_year)
#> [1] 39
nobs(fit1)
#> [1] 42
tryCatch(anova(fit1, fit_year), error = \(e) conditionMessage(e))
#> [1] "models were not all fitted to the same size of dataset"Source: R documentation: anova.lm
The first two tables show the order dependence: pages before genre gives pages a sum of squares of 21.13, genre before pages gives genre 100.06 and pages 0.04, because pages and genre are related across books, so the design is not orthogonal. The three-model comparison tests each step: adding pages to the null model gives F = 0.248 (p = 0.62), adding genre to that gives F = 0.186 (p = 0.97). The last lines show a requirement of anova() on several models: they must be fitted to the same rows. lm(days ~ year) uses 39 rows because year has more missing values, so comparing it with the 42-row model is refused with models were not all fitted to the same size of dataset; the script keeps only complete cases of days, pages and genre first so that the nested models do compare.
Missing values and na.action
na.action decides what lm() does with rows that have an NA in any variable of the model. The default is whatever options("na.action") holds, "na.omit" in a fresh R session (the R FAQ lists it among the differences from S: R sets it, S does not), which drops the incomplete rows silently except for a note in summary(). The alternatives are na.fail, which stops with missing values in object, na.exclude, which drops the rows for the fit but remembers them, and na.pass, which passes the NA values on (the model-fitting code will then fail on them). The practical difference between na.omit and na.exclude is the length of residuals(), fitted() and predict(): with na.exclude, naresid() pads them with NA for the dropped rows, so the results line up with the original data frame and can be stored as columns.
source("scenario.R")
# The session-wide default is na.omit.
getOption("na.action")
#> [1] "na.omit"
# books$pages has two NAs and books$year one, so the two fits differ in n.
fit_omit <- lm(pages ~ year, data = books)
nobs(fit_omit)
#> [1] 13
length(residuals(fit_omit))
#> [1] 13
class(fit_omit$na.action)
#> [1] "omit"
fit_omit$na.action
#> 6 7 15
#> 6 7 15
#> attr(,"class")
#> [1] "omit"
# na.exclude keeps the same fit but pads residuals() and fitted() with NA.
fit_excl <- lm(pages ~ year, data = books, na.action = na.exclude)
all.equal(coef(fit_omit), coef(fit_excl))
#> [1] TRUE
length(residuals(fit_excl))
#> [1] 16
class(fit_excl$na.action)
#> [1] "exclude"
round(residuals(fit_excl), 1)
#> 1 2 3 4 5 6 7 8 9 10 11
#> 85.7 -13.4 -47.1 -61.3 37.2 NA NA 1.0 84.6 -1.5 -185.6
#> 12 13 14 15 16
#> 147.1 -85.2 56.4 NA -18.0
# The stored component is short; naresid() pads it, and the extractor
# generics (residuals(), fitted(), predict()) call it for na.exclude.
raw <- fit_excl$fitted.values
length(raw)
#> [1] 13
padded <- naresid(fit_excl$na.action, raw)
length(padded)
#> [1] 16
identical(padded, fitted(fit_excl))
#> [1] TRUE
head(unname(padded), 8)
#> [1] 301.3470 445.3617 261.0944 259.3054 263.7779 NA NA 278.9845
# Because the length matches the data, the residuals can be added back.
books$resid <- residuals(fit_excl)
books[is.na(books$resid), c("title", "year", "pages")]
#> title year pages
#> 6 A Map of Salt 2008 NA
#> 7 Garden of Small Fires NA 96
#> 15 Ten Winters 2022 NA
# na.fail refuses to proceed; na.action can also be set for the session.
tryCatch(lm(pages ~ year, data = books, na.action = na.fail),
error = \(e) conditionMessage(e)
)
#> [1] "missing values in object"fit$na.action is itself an object: the integer positions of the omitted rows with class "omit" (or "exclude"). The extractor generics call naresid() (for residuals and fitted values) and napredict() (for predictions) with it, which is why residuals(fit_excl) has 16 elements while fit_excl$fitted.values, the stored component, has 13. A trap follows from the default: two models with different predictors may be fitted to different rows, so their AIC() values and anova() comparisons are not comparable (the error shown in the previous section is the protective case; AIC() does not check). Check nobs() for every model you compare, or drop the incomplete rows once, before fitting.
Leverage, influence and Cook’s distance
A fit can be dominated by a few rows. Three ideas separate the cases. The leverage of an observation is its diagonal entry in the hat matrix (hatvalues()), a number between 0 and 1 that depends only on the predictors; the leverages add up to the number of coefficients, so their mean is p over n, and a common rule of thumb flags those above twice the mean. A standardised residual (rstandard()) divides each residual by its estimated standard deviation, and a studentised residual (rstudent()) uses the estimate of sigma from the fit without that observation. Cook’s distance (cooks.distance()) combines the two: how far the fitted values move when one observation is left out. influence.measures() collects the dfbetas, dffits, covariance ratios, Cook’s distances and hat values together with logical flags, and summary() of the result prints only the observations flagged on at least one of them (R documentation: influence.measures).
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
fit <- lm(days ~ pages + genre, data = dur)
# Leverage: the diagonal of the hat matrix; the average is p / n.
h <- hatvalues(fit)
round(head(h, 4), 3)
#> 1 2 3 4
#> 0.077 0.077 0.077 0.077
c(mean = mean(h), p_over_n = fit$rank / nobs(fit))
#> mean p_over_n
#> 0.1666667 0.1666667
names(which(h > 2 * fit$rank / nobs(fit)))
#> [1] "22" "23" "34" "35"
# Standardised and studentised residuals.
round(head(rstandard(fit), 4), 3)
#> 1 2 3 4
#> -0.842 0.173 1.527 0.850
round(head(rstudent(fit), 4), 3)
#> 1 2 3 4
#> -0.839 0.171 1.558 0.847
# Cook's distance: how much the fit moves when one observation is dropped.
cd <- cooks.distance(fit)
round(sort(cd, decreasing = TRUE)[1:3], 3)
#> 47 27 15
#> 0.179 0.135 0.111
# influence.measures() collects dfbetas, dffits, covratio, Cook and hat;
# summary() prints only the rows flagged as potentially influential.
im <- influence.measures(fit)
class(im)
#> [1] "infl"
dim(im$infmat)
#> [1] 42 11
colnames(im$infmat)[1:3]
#> [1] "dfb.1_" "dfb.pags" "dfb.gnrfc"
flagged <- apply(im$is.inf, 1, any)
sum(flagged)
#> [1] 5
# summary() prints only the rows that are flagged on at least one measure.
summary(im)
#> Potentially influential observations of
#> lm(formula = days ~ pages + genre, data = dur) :
#>
#> dfb.1_ dfb.pags dfb.gnrfc dfb.gnrm dfb.gnr- dfb.gnrp dfb.gnrsf dffit cov.r
#> 22 0.15 -0.12 -0.16 -0.14 -0.13 -0.15 -0.11 0.17 1.92_*
#> 23 0.23 -0.19 -0.25 -0.22 -0.21 -0.23 -0.17 0.28 1.89_*
#> 34 0.00 0.00 0.00 0.00 0.49 0.00 0.00 0.84 2.12_*
#> 35 0.00 0.00 0.00 0.00 -0.49 0.00 0.00 -0.84 2.12_*
#> 52 0.01 -0.02 0.00 0.00 0.00 -0.01 0.03 -0.04 1.69_*
#> cook.d hat
#> 22 0.00 0.37
#> 23 0.01 0.37
#> 34 0.10 0.50
#> 35 0.10 0.50
#> 52 0.00 0.28Leverage in this data is set by how many loans fall in each genre cell, because genre has five dummy columns for 42 rows. Once the rows with missing pages are dropped, non-fiction has only two loans left, and a fitted value that depends on two observations gives each of them a leverage of 0.50; classic has five loans, and rows 22 and 23, which sit apart from the other three in pages, have leverage 0.37. Those four rows (22, 23, 34 and 35) are the ones above twice the mean (2 times 7 / 42, about 0.33), and they carry the flags for dffits and the covariance ratio. Cook’s distances are all small (the largest is 0.18), so no single loan drives the fit. The same quantities exist for a glm, computed from the final weighted least-squares fit; the R 4.6.0 NEWS says several influence measures for "glm" objects now use Pearson instead of deviance residuals, and that weighted.residuals() now returns the weighted working residuals for them.
Diagnostic plots
plot() on an lm produces up to six diagnostic plots, chosen by the which argument, and by default draws four (which = c(1, 2, 3, 5)): the manual’s chapter says "four plots", which is the default. They are, in order of the which number:
which |
Plot | What to look for |
|---|---|---|
1 |
Residuals against fitted values |
A flat band around zero with no curve. A curve says the mean model is wrong (add a transformation, a polynomial or a missing predictor); a funnel says the variance grows with the mean. |
2 |
Normal Q-Q plot of the standardised residuals |
Points near the line. Heavy tails or an S shape say the errors are not normal. For a |
3 |
Scale-location: the square root of the absolute standardised residuals against the fitted values |
A flat trend. A rising trend is non-constant variance. |
4 |
Cook’s distance for each observation |
Bars standing out from the rest; the |
5 |
Residuals against leverage |
Points in the upper or lower right corner, where both leverage and residual are large. Since R 4.6.0 the plot no longer draws a smoother. |
6 |
Cook’s distance against leverage divided by one minus leverage |
Points beyond the contour lines of constant Cook’s distance, which fan out from the origin. |
Source: R documentation: plot.lm
A script has no screen, so the plots are drawn to a device that is thrown away: svg() to a file in tempfile(), or pdf(NULL), which draws and keeps nothing. par(mfrow = c(2, 2)) puts four plots on one page; plot() returns NULL invisibly. The numbers each plot is drawn from are cheap to compute, and the script prints them so that the page does not depend on a picture. The plots are not reproduced here.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
fit <- lm(days ~ pages + genre, data = dur)
# The plots are drawn on the current device. A script has no screen, so draw
# to a file in the session's temporary directory (not reproduced here).
plot_file <- tempfile(fileext = ".svg")
svg(plot_file, width = 8, height = 8)
par(mfrow = c(2, 2))
plot(fit)
invisible(dev.off())
file.exists(plot_file)
#> [1] TRUE
file.size(plot_file) > 1000
#> [1] TRUE
# which = selects plots: 1 residuals vs fitted, 2 normal Q-Q, 3 scale-location,
# 4 Cook's distance, 5 residuals vs leverage, 6 Cook's distance vs leverage.
pdf(NULL)
par(mfrow = c(2, 3))
res <- plot(fit, which = 1:6, id.n = 3)
invisible(dev.off())
res
#> NULL
# The default is which = c(1, 2, 3, 5); this is the same, spelled out.
formals(stats:::plot.lm)$which
#> c(1, 2, 3, 5)
# What each plot shows, computed by hand (no drawing needed).
# 1: residuals against fitted values.
head(cbind(fitted = fitted(fit), resid = resid(fit)), 3)
#> fitted resid
#> 1 17.46541 -7.465406
#> 2 17.46541 1.534594
#> 3 17.46541 13.534594
# 2: standardised residuals against normal quantiles.
q <- qqnorm(rstandard(fit), plot.it = FALSE)
head(cbind(theoretical = q$x, sample = q$y), 3)
#> theoretical sample
#> 1 -0.60133175 -0.8424142
#> 2 0.02984524 0.1731672
#> 3 1.38299413 1.5272758
# 3: square root of absolute standardised residuals against fitted values.
head(sqrt(abs(rstandard(fit))), 3)
#> 1 2 3
#> 0.9178312 0.4161337 1.2358300
# 4: Cook's distance; the three largest are the labelled bars.
names(sort(cooks.distance(fit), decreasing = TRUE))[1:3]
#> [1] "47" "27" "15"
# 5: standardised residuals against leverage.
head(cbind(leverage = hatvalues(fit), std_resid = rstandard(fit)), 3)
#> leverage std_resid
#> 1 0.07727525 -0.8424142
#> 2 0.07727525 0.1731672
#> 3 0.07727525 1.5272758
# 6: Cook's distance against leverage / (1 - leverage).
h <- hatvalues(fit)
head(cbind(x = h / (1 - h), cooks = cooks.distance(fit)), 3)
#> x cooks
#> 1 0.0837468 0.008490270
#> 2 0.0837468 0.000358758
#> 3 0.0837468 0.027906482Source: R documentation: plot.lm
The NEWS records the history of the plots: R 4.0.0 changed plot.lm() for glm fits (the smoother’s iter.smooth now defaults to 0, and the residuals plotted are Pearson residuals, because deviance residuals do not in general have zero mean), and made Cook’s distance use the weighted residuals; R 4.3.0 replaced the normal Q-Q plot of a glm by the half-normal plot; R 4.6.0 removed the smoother from the residuals-versus-leverage plot (panel.raw = panel.smooth brings it back).
Weights, offsets and subsets
Three arguments change which observations count and how. weights multiplies each observation’s contribution to the sum of squares: use precision weights (the inverse of each observation’s variance), or frequency weights for rows that stand for several identical cases. Here every book gets the same total weight however many times it was loaned, by weighting each loan with one over the number of loans of its book, which reduces the dominance of often-loaned books. offset() is a term whose coefficient is fixed at one; it is how a known baseline enters, and it matters most in glm() for rates. subset selects rows with a logical or index expression evaluated in the data, before the model frame is built. A zero weight removes an observation from the fit (and from the degrees of freedom) while it stays in the model frame.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
# Weights: each book gets the same total weight, however often it was loaned.
copies <- table(dur$isbn)
dur$n_loans <- as.vector(copies[dur$isbn])
fit_w <- lm(days ~ pages, data = dur, weights = 1 / n_loans)
fit_u <- lm(days ~ pages, data = dur)
rbind(unweighted = coef(fit_u), weighted = coef(fit_w))
#> (Intercept) pages
#> unweighted 17.86927 0.006786498
#> weighted 18.14086 0.006196129
length(weights(fit_w))
#> [1] 42
head(weighted.residuals(fit_w), 3)
#> 1 2 3
#> -3.8648197 -0.1905851 4.7083944
# subset: the rows are chosen before the model frame is built.
fit_sf <- lm(days ~ pages, data = dur, subset = genre == "science fiction")
nobs(fit_sf)
#> [1] 13
fit_sf$call
#> lm(formula = days ~ pages, data = dur, subset = genre == "science fiction")
# offset: a term with a fixed coefficient of one, here 21 days as a baseline.
fit_o <- lm(days ~ pages + offset(rep(21, nrow(dur))), data = dur)
coef(fit_o)
#> (Intercept) pages
#> -3.130729022 0.006786498
fit_o2 <- lm(I(days - 21) ~ pages, data = dur)
all.equal(coef(fit_o), coef(fit_o2))
#> [1] TRUE
# A zero weight drops an observation from the fit, but not from the frame.
w <- ifelse(dur$genre == "poetry", 0, 1)
fit_0 <- lm(days ~ pages, data = dur, weights = w)
c(frame_rows = nrow(model.frame(fit_0)), df_residual = df.residual(fit_0))
#> frame_rows df_residual
#> 42 35Source: R documentation: lm
The unweighted slope of pages is 0.0068 and the weighted one 0.0062: weighting changes the estimates only slightly here. fit_sf used 13 loans of science fiction and keeps the subset argument in its call, so update() can reuse it. The offset model is the same as fitting days - 21 as the response, which the last check confirms. With zero weights for the 5 poetry loans, the model frame keeps all 42 rows but the residual degrees of freedom are 35: 42 rows, 5 given zero weight, 2 coefficients. R-devel (the upcoming R 4.7.0, not run here) adds an argument wtol to lm() and glm() for treating very small weights as zero, and fixes lm(rock)-style calls where the first argument is a data frame; both are listed under "Changes in R-devel" in the NEWS.
Updating a model, and choosing terms
update(fit, new.formula) refits the model with the stored call and a changed formula, in which . stands for the corresponding part of the old formula: . ~ . + year adds a term, log(.) ~ . transforms the response, . ~ . - genre drops one. Other arguments can be changed too (update(fit, subset = ...)). The manual explains the same facility, adds that a data= argument in the original call is passed on to update() and its allies, and names add1(), drop1() and step() as the functions for exploring incremental sequences of models (An Introduction to R: Updating fitted models). update() re-evaluates the stored call, so the data it names must still exist, with the same name, where update() is called.
add1() shows what adding each term in scope would do, drop1() what removing each term from the current model would do; both report the change in the residual sum of squares and in AIC, and test = "F" adds an F test. step() does this repeatedly: it adds or drops terms to minimise AIC and returns the final model with the path in its anova component. Be sceptical of stepwise selection: the p-values and intervals of the final model ignore the search that produced it, and with n = 42 and weak signal the search is fitting noise. It is shown here because it is in the manual, and the null model wins.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
dur <- dur[complete.cases(dur[, c("days", "pages", "genre")]), ]
fit <- lm(days ~ pages + genre, data = dur)
# update() re-fits the stored call with a changed formula: "." is the old part.
fit2 <- update(fit, . ~ . + year)
formula(fit2)
#> days ~ pages + genre + year
fit2$call
#> lm(formula = days ~ pages + genre + year, data = dur)
# Change the response, drop a term, or change another argument.
formula(update(fit, log(.) ~ .))
#> log(days) ~ pages + genre
formula(update(fit, . ~ . - genre))
#> days ~ pages
update(fit, subset = genre != "poetry")$call
#> lm(formula = days ~ pages + genre, data = dur, subset = genre !=
#> "poetry")
nobs(update(fit, subset = genre != "poetry"))
#> [1] 37
# add1() and drop1() show what each single change would do (AIC and F tests).
# year has missing values, so add1() warns that it used fewer rows.
tryCatch(add1(fit, scope = ~ . + year), warning = \(w) conditionMessage(w))
#> [1] "using the 39/42 rows from a combined fit"
dur_y <- dur[!is.na(dur$year), ]
fit_y <- lm(days ~ pages + genre, data = dur_y)
add1(fit_y, scope = ~ . + year, test = "F")
#> Single term additions
#>
#> Model:
#> days ~ pages + genre
#> Df Sum of Sq RSS AIC F value Pr(>F)
#> <none> 2865.4 181.58
#> year 1 85.283 2780.1 182.40 0.951 0.337
drop1(fit_y, test = "F")
#> Single term deletions
#>
#> Model:
#> days ~ pages + genre
#> Df Sum of Sq RSS AIC F value Pr(>F)
#> <none> 2865.4 181.58
#> pages 1 0.037 2865.4 179.58 0.0004 0.9840
#> genre 5 71.622 2937.0 172.54 0.1600 0.9753
# step() adds and drops terms by AIC; trace = 0 keeps the search quiet.
small <- lm(days ~ 1, data = dur)
best <- step(fit, scope = list(lower = small, upper = fit), trace = 0)
formula(best)
#> days ~ 1
best$anova
#> Step Df Deviance Resid. Df Resid. Dev AIC
#> 1 NA NA 35 2978.869 192.9884
#> 2 - genre 5 78.97640 40 3057.845 184.0874
#> 3 - pages 1 21.13073 41 3078.976 182.3767Source: R documentation: add1 and drop1
step() has its own reference page: R documentation: step.
add1() printed the warning using the 39/42 rows from a combined fit when year was offered, because year has missing values that the base fit did not have: the table would compare models fitted to different rows, so R refits on the common rows and says so. The rest of the script works on the 39 rows where year is known. step() removed genre (the AIC fell from 192.99 to 184.09) and then pages (to 182.38), ending with days ~ 1, the model with only a mean. That is the honest summary of this data.
Factor predictors, contrasts, interactions and polynomials
A factor in a formula becomes a set of columns built by contrasts. The default for an unordered factor is contr.treatment (compare each level with the first), the default for an ordered factor is contr.poly, and getOption("contrasts") shows the pair; the R FAQ notes that S uses Helmert contrasts for unordered factors, a deliberate difference. The coefficient meanings change with the coding while the fitted values do not. relevel() changes the reference level, C(f, contr.sum) or contrasts(f) <- ... change the coding, and a factor with fewer than two levels cannot be coded at all (An Introduction to R: Contrasts; the contrast matrices and C() are explained on Formulas).
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
dur <- dur[complete.cases(dur[, c("days", "pages", "genre")]), ]
dur$genre <- factor(dur$genre)
# The default coding for an unordered factor is treatment contrasts.
getOption("contrasts")
#> unordered ordered
#> "contr.treatment" "contr.poly"
contrasts(dur$genre)[1:3, ]
#> fiction mystery non-fiction poetry science fiction
#> classic 0 0 0 0 0
#> fiction 1 0 0 0 0
#> mystery 0 1 0 0 0
# The reference level is the first level in alphabetical order: "classic".
levels(dur$genre)
#> [1] "classic" "fiction" "mystery" "non-fiction"
#> [5] "poetry" "science fiction"
fit <- lm(days ~ genre, data = dur)
round(coef(fit), 2)
#> (Intercept) genrefiction genremystery
#> 21.40 -3.94 -0.65
#> genrenon-fiction genrepoetry genrescience fiction
#> 0.10 -1.40 -0.86
# relevel() changes the reference: intercept = mean of fiction.
dur$genre <- relevel(dur$genre, ref = "fiction")
fit_f <- lm(days ~ genre, data = dur)
round(coef(fit_f), 2)
#> (Intercept) genreclassic genremystery
#> 17.46 3.94 3.29
#> genrenon-fiction genrepoetry genrescience fiction
#> 4.04 2.54 3.08
# Group means are the same whatever the coding.
tapply(dur$days, dur$genre, mean) |> round(2)
#> fiction classic mystery non-fiction poetry
#> 17.46 21.40 20.75 21.50 20.00
#> science fiction
#> 20.54
all.equal(unname(fitted(fit)), unname(fitted(fit_f)))
#> [1] TRUE
# Sum contrasts: coefficients are deviations from the average of the means.
fit_s <- lm(days ~ C(genre, contr.sum), data = dur)
round(coef(fit_s), 2)
#> (Intercept) C(genre, contr.sum)1 C(genre, contr.sum)2
#> 20.27 -2.81 1.13
#> C(genre, contr.sum)3 C(genre, contr.sum)4 C(genre, contr.sum)5
#> 0.47 1.22 -0.28
# A factor with a single level cannot be used.
one <- dur[dur$genre == "fiction", ]
one$genre <- factor(one$genre)
tryCatch(lm(days ~ genre, data = one), error = \(e) conditionMessage(e))
#> [1] "contrasts can be applied only to factors with 2 or more levels"Source: R documentation: relevel
With treatment contrasts and "classic" as the reference, the intercept (21.40) is the mean duration of classic loans; after relevel(ref = "fiction") it is 17.46, the mean of fiction loans, and every coefficient is a difference from that. The group means printed by tapply() are the same in both fits, and so are the fitted values. With sum contrasts the intercept is the unweighted average of the six group means (20.27) and the coefficients are deviations from it, with the last level’s deviation implied (the coefficients sum to zero). The error contrasts can be applied only to factors with 2 or more levels is what a factor with one level produces, usually after subsetting a data frame without droplevels().
Interactions and polynomial terms use the formula operators from the formulas page; the models themselves are ordinary. a * b means a + b + a:b, so an interaction without its main effects (a:b alone) forces a common intercept, which is rarely what is meant. In a formula ^ means crossing, so a quadratic term is written I(pages^2) (R FAQ 3.3.2 gives the same warning), or with poly(pages, 2), which builds orthogonal polynomials: the fitted values are identical to those from I(), the coefficients differ, and predict() re-uses the stored basis so that new data is coded consistently. poly(pages, 2, raw = TRUE) gives the same coefficients as I().
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
dur <- dur[complete.cases(dur[, c("days", "pages", "genre")]), ]
# I() protects arithmetic: pages^2 inside I() is a square, not a crossing.
fit_sq <- lm(days ~ pages + I(pages^2), data = dur)
round(coef(fit_sq), 6)
#> (Intercept) pages I(pages^2)
#> 21.192553 -0.022325 0.000054
# poly() builds orthogonal polynomials: same fit, different coefficients.
fit_poly <- lm(days ~ poly(pages, 2), data = dur)
round(coef(fit_poly), 3)
#> (Intercept) poly(pages, 2)1 poly(pages, 2)2
#> 19.690 4.597 3.705
all.equal(fitted(fit_sq), fitted(fit_poly))
#> [1] TRUE
# poly(raw = TRUE) gives the same coefficients as I().
fit_raw <- lm(days ~ poly(pages, 2, raw = TRUE), data = dur)
all.equal(unname(coef(fit_raw)), unname(coef(fit_sq)))
#> [1] TRUE
# predict() re-uses the basis of poly(), so new data is coded the same way.
predict(fit_poly, newdata = data.frame(pages = c(100, 300)))
#> 1 2
#> 19.50127 19.36620
# Interaction: the slope of pages may differ between fiction and the rest.
dur$is_fiction <- dur$genre == "fiction"
fit_int <- lm(days ~ pages * is_fiction, data = dur)
round(coef(fit_int), 4)
#> (Intercept) pages is_fictionTRUE
#> 20.4357 0.0009 -24.8699
#> pages:is_fictionTRUE
#> 0.1035
colnames(model.matrix(fit_int))
#> [1] "(Intercept)" "pages" "is_fictionTRUE"
#> [4] "pages:is_fictionTRUE"
# y ~ a + b + a:b is the same model as y ~ a * b.
fit_int2 <- lm(days ~ pages + is_fiction + pages:is_fiction, data = dur)
all.equal(coef(fit_int), coef(fit_int2))
#> [1] TRUE
# Marginality: an interaction without its main effects forces the lines
# through a common intercept; R fits it, but the model is rarely meant.
coef(lm(days ~ pages:is_fiction, data = dur))
#> (Intercept) pages:is_fictionFALSE pages:is_fictionTRUE
#> 20.000674141 0.002140103 -0.011477717
# Is the interaction needed?
anova(lm(days ~ pages + is_fiction, data = dur), fit_int)
#> Analysis of Variance Table
#>
#> Model 1: days ~ pages + is_fiction
#> Model 2: days ~ pages * is_fiction
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 39 2984.3
#> 2 38 2951.7 1 32.678 0.4207 0.5205Source: R documentation: poly
The coefficient of pages is 0.0009 for loans that are not fiction and 0.1035 larger for the fiction ones; the F test of the interaction has p = 0.52. A quadratic in pages has coefficients -0.0223 and 0.000054: tiny numbers because pages ranges over hundreds, which is a reason to prefer poly() (or to rescale the predictor) when the raw coefficients are unreadable.
Errors and warnings from lm()
The messages below are the ones that real data produces most often. Each one is caught with tryCatch() so that the script exits cleanly (see Conditions and Error Handling).
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
# A name that is neither a column of data nor a variable in scope.
tryCatch(
lm(pages ~ no_such_column, data = books),
error = \(e) conditionMessage(e)
)
#> [1] "object 'no_such_column' not found"
# Variables of different lengths in the calling environment.
x <- 1:10
y <- 1:9
tryCatch(lm(y ~ x), error = \(e) conditionMessage(e))
#> [1] "variable lengths differ (found for 'x')"
# A response that is all NA after na.action.
tryCatch(lm(c(NA_real_, NA_real_) ~ c(1, 2)), error = \(e) conditionMessage(e))
#> [1] "0 (non-NA) cases"
# Collinear columns: the redundant coefficient is NA, not an error.
books$pages2 <- 2 * books$pages
fit <- lm(year ~ pages + pages2, data = books)
coef(fit)
#> (Intercept) pages pages2
#> 2077.0839354 -0.2872619 NA
alias(fit)$Complete
#> (Intercept) pages
#> pages2 0 2
summary(fit)$coefficients
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2077.0839354 43.7972181 47.425020 4.489586e-14
#> pages -0.2872619 0.1472853 -1.950377 7.706806e-02
# predict() on new data that breaks the collinearity: the 4.3.0+ warning.
new <- data.frame(pages = 200, pages2 = 500)
tryCatch(predict(fit, newdata = new), warning = \(w) conditionMessage(w))
#> [1] "prediction from rank-deficient fit; attr(*, \"non-estim\") has doubtful cases"
predict(fit, newdata = new, rankdeficient = "NA")
#> 1
#> NA
predict(fit, newdata = data.frame(pages = 200, pages2 = 400))
#> 1
#> 2019.632
# A character response is coerced to numeric: a warning, then an error.
tryCatch(lm(title ~ pages, data = books), warning = \(w) conditionMessage(w))
#> [1] "NAs introduced by coercion"
suppressWarnings(tryCatch(lm(title ~ pages, data = books),
error = \(e) conditionMessage(e)
))
#> [1] "NA/NaN/Inf in 'y'"Source: R documentation: lm
-
object 'no_such_column' not found: a name in the formula is neither a column ofdatanor a variable in the formula’s environment. Check the spelling, and check thatdata =was given. -
variable lengths differ (found for 'x'): the variables came from the workspace and have different lengths (here 10 and 9); keep the variables as columns of one data frame so that this cannot happen. -
0 (non-NA) cases: afterna.action, nothing is left to fit. -
A coefficient reported as
NAis not an error:lm()hassingular.ok = TRUEby default, so a column that is a linear combination of earlier ones (pages2is twicepages) is dropped and its coefficient isNA, with no warning.alias(fit)shows the dependency, andsummary()printsCoefficients: (1 not defined because of singularities). The test is numerical, with a tolerance of1e-7(thetolargument of the defaultmethod = "qr"). -
prediction from rank-deficient fit; attr(*, "non-estim") has doubtful cases: a prediction is requested for new data for which the dropped coefficient matters. The wording and the behaviour changed in R 4.3.0. Until thenpredict.lm()warned for every prediction from a rank-deficient fit; since R 4.3.0 there is arankdeficientargument whose default,"warnif", warns only when some rows ofnewdataare not estimable (the NEWS describes the feature as still somewhat experimental),"NA"returnsNAfor those rows, and the last call shows that an estimable row (one that keeps the collinearity,pages2equal to twicepages) predicts without any warning. -
NAs introduced by coercionfollowed byNA/NaN/Inf in 'y': a character column used as the response is coerced to numeric, which producesNA, and the fit stops. Convert the column deliberately, or model a count or a factor withglm().
The R FAQ’s list of differences from S also warns about lm(y ~ x^3): in a formula ^ is the crossing operator, so the cube must be protected with I(x^3) (R FAQ 3.3.2: Models).
Generalised linear models
A linear model assumes a normal response with a constant variance and a mean that is a linear function of the predictors. A generalised linear model (GLM) keeps the linear predictor and relaxes the other two assumptions. The manual states the structure in three assumptions (An Introduction to R: Generalized linear models):
-
The predictors act on the response only through a single linear predictor,
eta = beta1 x1 + ... + betap xp. -
The response has a distribution from the exponential family (normal, binomial, Poisson, gamma, inverse Gaussian) whose variance is a function of its mean and possibly a scale parameter.
-
The mean is a smooth invertible function of the linear predictor,
mu = m(eta), and the inverse function is the link.
The combination of a distribution and a link is a family. The chapter lists the families and the links that are available for each (An Introduction to R: Families); the table below lists the default link of each family first, which is not the order of the manual’s table, and adds cauchit, which ?family lists for the binomial family.
| Family | Links | Typical use |
|---|---|---|
|
logit, probit, cloglog, cauchit, log |
Yes or no outcomes, and counts of successes out of trials |
|
identity, log, inverse |
The same model as |
|
inverse, identity, log |
Positive, right-skewed continuous responses: durations, costs |
|
|
Positive responses with a variance that grows faster than the Gamma’s |
|
log, identity, sqrt |
Counts |
|
chosen with |
Quasi-likelihood: only the mean-variance relation is specified, and the dispersion is estimated |
Source: An Introduction to R: Families
glm(formula, family = family, data = data) fits the model. The family argument is a family object or a function that generates one (binomial and binomial() are equivalent; binomial(link = "probit") chooses a non-default link). Under the hood glm() calls glm.fit(), which uses iteratively reweighted least squares (IRLS): starting from fitted values derived from the data, it repeatedly forms a working response and weights from the current linear predictor, solves a weighted least-squares problem, and stops when the deviance changes by less than epsilon (1e-8) or after maxit (25) iterations; glm.control() sets both. Each pass is an ordinary weighted least-squares fit (a working response z built from the linear predictor and the residuals, with weights mu * (1 - mu) for the logistic model), which is why a glm object inherits from lm and contains a QR decomposition.
Logistic regression for overdue loans
A loan is overdue when it was kept longer than loan_days, 21 days. The response overdue is logical (a logical, a 0/1 numeric, or a factor whose first level means failure are all accepted by family = binomial) and the link is the logit, the log of the odds. The coefficients are therefore changes in log odds of an overdue loan per unit of the predictor, with the same treatment-contrast coding of genre as before. The printed object shows the call, the coefficients, the null and residual deviances and the AIC, and summary() adds the coefficient table, now with z values (normal approximation) rather than t values, because the dispersion is fixed at one.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
# A loan is overdue when it was kept longer than loan_days (21) days.
dur$overdue <- dur$days > loan_days
table(dur$overdue)
#>
#> FALSE TRUE
#> 26 24
# family = binomial: logit link, and a logical (or 0/1, or factor) response.
gfit <- glm(overdue ~ pages + genre, family = binomial, data = dur)
gfit
#>
#> Call: glm(formula = overdue ~ pages + genre, family = binomial, data = dur)
#>
#> Coefficients:
#> (Intercept) pages genrefiction
#> 3.957005 -0.009405 -2.456090
#> genremystery genrenon-fiction genrepoetry
#> -1.529257 -2.311059 -2.678581
#> genrescience fiction
#> -0.936101
#>
#> Degrees of Freedom: 41 Total (i.e. Null); 35 Residual
#> (8 observations deleted due to missingness)
#> Null Deviance: 57.84
#> Residual Deviance: 55.46 AIC: 69.46
summary(gfit)
#>
#> Call:
#> glm(formula = overdue ~ pages + genre, family = binomial, data = dur)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 3.957005 3.916816 1.010 0.312
#> pages -0.009405 0.009855 -0.954 0.340
#> genrefiction -2.456090 2.064504 -1.190 0.234
#> genremystery -1.529257 1.874970 -0.816 0.415
#> genrenon-fiction -2.311059 2.686954 -0.860 0.390
#> genrepoetry -2.678581 3.172557 -0.844 0.399
#> genrescience fiction -0.936101 1.132083 -0.827 0.408
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> Null deviance: 57.843 on 41 degrees of freedom
#> Residual deviance: 55.456 on 35 degrees of freedom
#> (8 observations deleted due to missingness)
#> AIC: 69.456
#>
#> Number of Fisher Scoring iterations: 4
#>Source: R documentation: glm
Of the 50 returned loans, 24 were overdue. The summary says what the linear-model one said: nothing is significant. The null deviance (57.84 on 41 degrees of freedom) is the deviance of the intercept-only model, the residual deviance (55.46 on 35) that of the fitted one, and their difference, 2.39 on 6 degrees of freedom, is a chi-squared statistic for the whole model: far from significant. The AIC is 69.46, higher than that of simpler models (the next section compares them), and the fit needed four Fisher scoring iterations (the IRLS iterations).
Odds ratios, intervals and predictions
exp(coef(fit)) converts the log odds into odds ratios: the multiplicative change in the odds of an overdue loan. An odds ratio of 0.991 for pages means that each extra page multiplies the odds of an overdue loan by 0.991 (0.39 for 100 more pages), with an interval, 0.97 to 1.01, that contains 1. confint(fit) for a glm is the profile likelihood interval, which can be asymmetric and is more reliable than the Wald interval from confint.default() (the estimate plus or minus 1.96 standard errors); since R 4.4.0 it is a method in stats, where it was previously a stub that loaded MASS, and it prints a message Waiting for profiling to be done... while it works. The script collects the message so that it does not reach the console.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
dur$overdue <- dur$days > loan_days
gfit <- glm(overdue ~ pages + genre, family = binomial, data = dur)
# Coefficients are log odds ratios; exp() gives odds ratios.
round(exp(coef(gfit)), 3)
#> (Intercept) pages genrefiction
#> 52.300 0.991 0.086
#> genremystery genrenon-fiction genrepoetry
#> 0.217 0.099 0.069
#> genrescience fiction
#> 0.392
# The profile-likelihood interval (stats::confint.glm since R 4.4.0) prints a
# message while it profiles; collect it instead of letting it reach the console.
msgs <- character()
ci <- withCallingHandlers(confint(gfit), message = \(m) {
msgs <<- c(msgs, conditionMessage(m))
invokeRestart("muffleMessage")
})
msgs
#> [1] "Waiting for profiling to be done...\n"
round(exp(ci), 3)
#> 2.5 % 97.5 %
#> (Intercept) 0.041 462823.767
#> pages 0.969 1.009
#> genrefiction 0.001 3.511
#> genremystery 0.003 6.977
#> genrenon-fiction 0.000 15.740
#> genrepoetry 0.000 22.845
#> genrescience fiction 0.034 3.513
# The Wald interval from the coefficient table is quicker and symmetric.
round(exp(confint.default(gfit)), 3)
#> 2.5 % 97.5 %
#> (Intercept) 0.024 112851.048
#> pages 0.972 1.010
#> genrefiction 0.001 4.905
#> genremystery 0.005 8.547
#> genrenon-fiction 0.001 19.207
#> genrepoetry 0.000 34.450
#> genrescience fiction 0.043 3.607
# predict(): type = "link" (the default) is the log odds, "response" the
# probability.
new <- data.frame(pages = c(150, 300, 450), genre = "fiction")
predict(gfit, newdata = new)
#> 1 2 3
#> 0.09010392 -1.32070740 -2.73151872
predict(gfit, newdata = new, type = "response")
#> 1 2 3
#> 0.52251075 0.21070062 0.06113893
plogis(predict(gfit, newdata = new))
#> 1 2 3
#> 0.52251075 0.21070062 0.06113893
round(head(fitted(gfit), 4), 3)
#> 1 2 3 4
#> 0.375 0.375 0.375 0.375
# Predicted classes at the 0.5 threshold against the observed ones.
observed <- model.response(model.frame(gfit))
table(observed, predicted = fitted(gfit) > 0.5)
#> predicted
#> observed FALSE TRUE
#> FALSE 20 3
#> TRUE 12 7Source: R documentation: confint
Reading the output: the odds ratio of 52.3 for the intercept is the odds of an overdue classic loan with zero pages, an extrapolation, and its interval runs from 0.04 to over 400,000 because the data say almost nothing about it. The Wald and profile intervals agree in sign and differ in the tails (for genrefiction 0.001 to 4.9 against 0.001 to 3.5). predict() has a type argument that matters for GLMs: type = "link" (the default) is on the linear-predictor scale, the log odds; type = "response" is the predicted probability, the inverse link applied (plogis() for the logit); type = "terms" gives the contribution of each term. A fiction book of 150 pages has a predicted probability of 0.52 of being overdue, and one of 450 pages 0.06. The cross-tabulation at a 0.5 threshold classifies 27 of 42 loans correctly, which is slightly better than guessing the majority class (23 of 42 loans were on time), and it is evaluated on the data the model was fitted to.
Deviance, anova and residuals
The deviance of a GLM is twice the difference between the log-likelihood of a saturated model and that of the fitted one; it plays the part the residual sum of squares plays for lm(). anova(fit) on a glm is a sequential analysis of deviance. The R 4.4.0 NEWS changed its default: it now computes test statistics and p-values, a chi-squared test when the dispersion is fixed and an F test when it is estimated, and test = FALSE suppresses them (so the explicit test = "Chisq" in older code is no longer needed, but it is harmless). anova(g0, g1, g2, test = "Chisq") compares nested models with the likelihood-ratio test, and the drop in deviance by hand against pchisq() gives the same number.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
dur$overdue <- dur$days > loan_days
dur <- dur[complete.cases(dur[, c("pages", "genre")]), ]
gfit <- glm(overdue ~ pages + genre, family = binomial, data = dur)
# anova() on a glm is sequential too; test = "Chisq" adds the p-values
# (the default since R 4.4.0 for a glm).
anova(gfit, test = "Chisq")
#> Analysis of Deviance Table
#>
#> Model: binomial, link: logit
#>
#> Response: overdue
#>
#> Terms added sequentially (first to last)
#>
#>
#> Df Deviance Resid. Df Resid. Dev Pr(>Chi)
#> NULL 41 57.843
#> pages 1 0.38841 40 57.454 0.5331
#> genre 5 1.99880 35 55.456 0.8493
# Giving test a false value suppresses the p-values.
anova(gfit, test = FALSE)
#> Analysis of Deviance Table
#>
#> Model: binomial, link: logit
#>
#> Response: overdue
#>
#> Terms added sequentially (first to last)
#>
#>
#> Df Deviance Resid. Df Resid. Dev
#> NULL 41 57.843
#> pages 1 0.38841 40 57.454
#> genre 5 1.99880 35 55.456
# Nested models: the likelihood-ratio (deviance) test.
g0 <- glm(overdue ~ 1, family = binomial, data = dur)
g1 <- glm(overdue ~ pages, family = binomial, data = dur)
anova(g0, g1, gfit, test = "Chisq")
#> Analysis of Deviance Table
#>
#> Model 1: overdue ~ 1
#> Model 2: overdue ~ pages
#> Model 3: overdue ~ pages + genre
#> Resid. Df Resid. Dev Df Deviance Pr(>Chi)
#> 1 41 57.843
#> 2 40 57.454 1 0.38841 0.5331
#> 3 35 55.456 5 1.99880 0.8493
# The same test by hand: the drop in deviance against a chi-squared.
dev_drop <- deviance(g1) - deviance(gfit)
c(drop = dev_drop, df = g1$df.residual - gfit$df.residual)
#> drop df
#> 1.998799 5.000000
pchisq(dev_drop, df = 5, lower.tail = FALSE)
#> [1] 0.8493112
# AIC and BIC, as for lm.
AIC(g0, g1, gfit)
#> df AIC
#> g0 1 59.84283
#> g1 2 61.45442
#> gfit 7 69.45562
BIC(g0, g1, gfit)
#> df BIC
#> g0 1 61.58050
#> g1 2 64.92976
#> gfit 7 81.61931
# Residuals come in several types.
r_dev <- residuals(gfit, type = "deviance")
r_prs <- residuals(gfit, type = "pearson")
r_res <- residuals(gfit, type = "response")
r_wrk <- residuals(gfit, type = "working")
round(rbind(
deviance = r_dev[1:4], pearson = r_prs[1:4],
response = r_res[1:4], working = r_wrk[1:4]
), 3)
#> 1 2 3 4
#> deviance -0.969 -0.969 1.401 1.401
#> pearson -0.774 -0.774 1.292 1.292
#> response -0.375 -0.375 0.625 0.625
#> working -1.599 -1.599 2.668 2.668
# The squared deviance residuals add up to the deviance; the squared Pearson
# residuals to Pearson's chi-squared statistic, the basis of the dispersion.
all.equal(sum(r_dev^2), deviance(gfit))
#> [1] TRUE
c(pearson = sum(r_prs^2), dispersion = sum(r_prs^2) / df.residual(gfit))
#> pearson dispersion
#> 41.802869 1.194368
summary(gfit)$dispersion
#> [1] 1Source: R documentation: anova.glm
residuals(fit, type =) has four kinds for a GLM. The deviance residuals (the default) are the signed square roots of each observation’s contribution to the deviance, so their squares add up to the deviance. The Pearson residuals divide the raw difference by the standard deviation implied by the variance function, and their squares add up to Pearson’s chi-squared statistic. Response residuals are observed minus fitted on the response scale, and working residuals are those of the last IRLS step. The Pearson statistic divided by the residual degrees of freedom (1.19 here) estimates the dispersion: a binomial or Poisson model should give about 1, and a value much larger than 1 is overdispersion, where quasibinomial or quasipoisson (or a mixed model) is the usual remedy.
Warnings: separation and non-convergence
Two warnings appear whenever a logistic regression meets a pattern the data cannot support, and they are the two a user meets first. glm.fit: fitted probabilities numerically 0 or 1 occurred means that some fitted probabilities are indistinguishable from 0 or 1 in floating point. If it appears together with glm.fit: algorithm did not converge, the usual cause is separation: a predictor (or a combination of predictors) classifies the observations perfectly, so the likelihood has no finite maximum and IRLS pushes the coefficients to larger and larger values until it runs out of iterations. The estimates and standard errors are then meaningless (the Wald z statistic and p-value are useless), and the remedies are to remove or coarsen the offending predictor, to use a penalised method (Firth’s bias-reduced logistic regression is available in add-on packages such as brglm2 or logistf; named only, not used here), or to collect more data. The overdue flag is defined from the duration, so the duration separates it perfectly. The second warning alone, algorithm did not converge, means that the limit of maxit iterations was reached; raising maxit through glm.control() may help, and if it does not the model is probably badly specified.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
dur$overdue <- dur$days > loan_days
# Warnings do not stop glm(), so collect them while the model is fitted.
collect <- function(expr) {
msgs <- character()
value <- withCallingHandlers(expr, warning = \(w) {
msgs <<- c(msgs, conditionMessage(w))
invokeRestart("muffleWarning")
})
list(value = value, warnings = msgs)
}
# Complete separation: overdue is defined from days, so days predicts it
# perfectly.
sep <- collect(glm(overdue ~ days, family = binomial, data = dur))
sep$warnings
#> [1] "glm.fit: algorithm did not converge"
#> [2] "glm.fit: fitted probabilities numerically 0 or 1 occurred"
range(round(fitted(sep$value), 6))
#> [1] 0 1
# The estimate is not finite, only large, and its standard error explodes:
# the Wald z statistic is useless here.
coef(summary(sep$value))[, "Std. Error"] > 1000
#> (Intercept) days
#> TRUE TRUE
sep$value$converged
#> [1] FALSE
# Too few iterations: "algorithm did not converge".
few <- collect(glm(overdue ~ pages + genre,
family = binomial, data = dur,
control = glm.control(maxit = 2)
))
few$warnings
#> [1] "glm.fit: algorithm did not converge"
c(converged = few$value$converged, iter = few$value$iter)
#> converged iter
#> 0 2Source: R documentation: glm.control
The script collects the warnings with withCallingHandlers() and invokeRestart("muffleWarning"), so it can print them in order; glm() returns normally after a warning, and fit$converged is the programmatic check. Both texts are produced by glm.fit(), which is why they begin with glm.fit:. In the separation case the fitted probabilities run from 0 to 1 (to six decimals), both standard errors are larger than 1000, and fit$converged is FALSE.
Counts with Poisson and quasi-Poisson
Loans per month is a count. family = poisson models the log of the expected count as linear in the predictors, so exp(coef) is a rate ratio. The script first builds the counts: table() leaves out months with no loans, which would bias the model, so the months are given as the levels of a factor. The time index t is the month number.
source("scenario.R")
# Loans per month: a count, so family = poisson (log link by default).
# table() drops months with no loans, so give it every month as a level.
month_starts <- seq(as.Date("2025-01-01"), by = "month", length.out = 19)
months <- format(month_starts[-19], "%Y-%m")
counts <- table(factor(format(loans$loan_date, "%Y-%m"), levels = months))
monthly <- data.frame(month = months, n = as.vector(counts))
monthly$t <- seq_len(nrow(monthly))
head(monthly, 4)
#> month n t
#> 1 2025-01 4 1
#> 2 2025-02 2 2
#> 3 2025-03 3 3
#> 4 2025-04 3 4
sum(monthly$n)
#> [1] 59
summary(monthly$n)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.000 3.000 3.000 3.278 4.000 7.000
pfit <- glm(n ~ t, family = poisson, data = monthly)
summary(pfit)
#>
#> Call:
#> glm(formula = n ~ t, family = poisson, data = monthly)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 1.29028 0.26531 4.863 1.15e-06 ***
#> t -0.01103 0.02512 -0.439 0.661
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> (Dispersion parameter for poisson family taken to be 1)
#>
#> Null deviance: 17.040 on 17 degrees of freedom
#> Residual deviance: 16.847 on 16 degrees of freedom
#> AIC: 72.687
#>
#> Number of Fisher Scoring iterations: 5
#>
# exp(coefficient) is a rate ratio: the monthly growth in loans.
exp(coef(pfit))
#> (Intercept) t
#> 3.6338209 0.9890337
round(exp(confint.default(pfit)), 3)
#> 2.5 % 97.5 %
#> (Intercept) 2.160 6.112
#> t 0.942 1.039
# Dispersion: Pearson chi-squared per residual degree of freedom.
sum(residuals(pfit, type = "pearson")^2) / df.residual(pfit)
#> [1] 0.9423944
# quasipoisson estimates the dispersion; the estimates stay the same and
# the standard errors are scaled by the square root of the dispersion.
qfit <- glm(n ~ t, family = quasipoisson, data = monthly)
all.equal(coef(pfit), coef(qfit))
#> [1] TRUE
summary(qfit)$dispersion
#> [1] 0.9423944
round(rbind(
poisson = coef(summary(pfit))[, "Std. Error"],
quasi = coef(summary(qfit))[, "Std. Error"]
), 4)
#> (Intercept) t
#> poisson 0.2653 0.0251
#> quasi 0.2576 0.0244
# A quasi-likelihood is not a likelihood, so there is no AIC.
AIC(qfit)
#> [1] NA
# An offset turns counts into rates: loans per member-day of exposure
# (here a made-up exposure, the number of days in each month).
monthly$days <- as.numeric(diff(month_starts))
ofit <- glm(n ~ t + offset(log(days)), family = poisson, data = monthly)
round(coef(ofit), 4)
#> (Intercept) t
#> -2.1202 -0.0112Source: R documentation: glm
There are 59 dated loans over 18 months, 3.3 per month on average, and one month with no loans at all (the Min. of 0), which is the month a plain table() would have dropped. The slope is -0.011 on the log scale, a rate ratio of 0.989 per month with a Wald confidence interval from 0.94 to 1.04: no trend. The Poisson model assumes that the variance equals the mean, and the Pearson statistic per residual degree of freedom estimates how far from that the data are: 0.94. family = quasipoisson fits the same coefficients (the all.equal() check is true) with the standard errors scaled by the square root of the estimated dispersion; here the dispersion is slightly below 1, so the standard errors shrink a little (0.0244 instead of 0.0251 for the slope), and one would normally keep the Poisson ones rather than claim more precision than the model assumes. A quasi-likelihood is not a likelihood, so AIC(qfit) is NA. An offset(log(exposure)) turns counts into rates: with the number of days in each month as the exposure (made up for the illustration), the slope is the log change of loans per day.
Positive durations with the Gamma family
Loan durations are positive and right-skewed. A Gamma family with a log link models the mean duration as the exponential of a linear predictor, with a variance proportional to the square of the mean (a constant coefficient of variation). The default link of Gamma() is the inverse, its canonical link, which can predict negative means; the log link is the usual choice. The dispersion is estimated rather than fixed, so the summary reports t values and the line Dispersion parameter for Gamma family taken to be 0.1992879; a zero or negative response is rejected with non-positive values not allowed for the 'Gamma' family.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
# Durations are positive and right-skewed: a Gamma family with a log link.
dur <- dur[complete.cases(dur[, c("pages", "genre")]), ]
min(dur$days)
#> [1] 1
gam_fit <- glm(days ~ pages, family = Gamma(link = "log"), data = dur)
summary(gam_fit)
#>
#> Call:
#> glm(formula = days ~ pages, family = Gamma(link = "log"), data = dur)
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2.8894961 0.1898057 15.223 <2e-16 ***
#> pages 0.0003354 0.0006591 0.509 0.614
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> (Dispersion parameter for Gamma family taken to be 0.1992879)
#>
#> Null deviance: 13.428 on 41 degrees of freedom
#> Residual deviance: 13.375 on 40 degrees of freedom
#> AIC: 316.26
#>
#> Number of Fisher Scoring iterations: 5
#>
# With the log link, exp(coef) multiplies the mean duration.
round(exp(coef(gam_fit)), 4)
#> (Intercept) pages
#> 17.9842 1.0003
# The dispersion is estimated (the Gamma family has a free scale parameter).
summary(gam_fit)$dispersion
#> [1] 0.1992879
# The default link of the Gamma family is the inverse (the canonical link).
Gamma()$link
#> [1] "inverse"
glm(days ~ pages, family = Gamma, data = dur)$family$link
#> [1] "inverse"
# Compare with the Gaussian fit on the same rows: same mean model, other scale.
gau <- glm(days ~ pages, family = gaussian, data = dur)
c(gaussian = AIC(gau), gamma_log = AIC(gam_fit))
#> gaussian gamma_log
#> 305.2783 316.2565
# A zero or negative response is rejected by the Gamma family.
bad <- dur
bad$days[1] <- 0
tryCatch(glm(days ~ pages, family = Gamma(link = "log"), data = bad),
error = \(e) conditionMessage(e)
)
#> [1] "non-positive values not allowed for the 'Gamma' family"Source: R documentation: family
exp(coef) multiplies the mean: the intercept is a mean of 18 days at zero pages, and each extra page multiplies the mean duration by 1.0003. The AIC of the Gaussian fit (305.3) is lower than that of the Gamma log-link fit (316.3) on the same rows: AIC compares models for the same response, and here the roughly uniform durations (1 to 37 days) suit the normal model better than the Gamma.
Grouped binomial data
When the unit of analysis is a group, the binomial response is a two-column matrix of successes and failures, or a proportion with the number of trials as weights; the manual’s example of the blind men of Kalythos uses the matrix form (An Introduction to R: The glm() function). The three formulations give the same coefficients. Grouping also makes goodness-of-fit testing possible, because the residual deviance of a model with fewer parameters than groups can be compared with a chi-squared distribution (the approximation needs expected counts that are not tiny, and these counts are, so treat the p-value as a demonstration).
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
dur$overdue <- dur$days > loan_days
# Grouped binomial data: successes and failures per genre.
tab <- table(dur$genre, dur$overdue)
counts <- data.frame(
genre = rownames(tab), late = tab[, "TRUE"],
on_time = tab[, "FALSE"], row.names = NULL
)
counts$n <- counts$late + counts$on_time
counts
#> genre late on_time n
#> 1 classic 3 2 5
#> 2 fiction 8 10 18
#> 3 mystery 2 2 4
#> 4 non-fiction 3 2 5
#> 5 poetry 3 2 5
#> 6 science fiction 5 8 13
# The response may be a two-column matrix of successes and failures.
grp <- glm(cbind(late, on_time) ~ genre, family = binomial, data = counts)
round(coef(grp), 3)
#> (Intercept) genrefiction genremystery
#> 0.405 -0.629 -0.405
#> genrenon-fiction genrepoetry genrescience fiction
#> 0.000 0.000 -0.875
# The same model from a proportion and the number of trials as weights.
counts$prop <- counts$late / counts$n
grp2 <- glm(prop ~ genre, family = binomial, weights = n, data = counts)
all.equal(coef(grp), coef(grp2))
#> [1] TRUE
# And from the 0/1 rows: the same estimates (the deviance differs, as the
# grouped data has fewer, larger units).
rows <- glm(overdue ~ genre, family = binomial, data = dur)
all.equal(coef(grp), coef(rows), tolerance = 1e-6)
#> [1] TRUE
round(c(grouped = deviance(grp), rows = deviance(rows)), 3)
#> grouped rows
#> 0.000 67.789
# Six parameters for six groups leaves no residual degrees of freedom.
df.residual(grp)
#> [1] 0
# Grouped data allow a goodness-of-fit test of the intercept-only model: its
# residual deviance is compared with a chi-squared on the residual df.
null <- glm(cbind(late, on_time) ~ 1, family = binomial, data = counts)
c(deviance = deviance(null), df = df.residual(null))
#> deviance df
#> 1.445313 5.000000
pchisq(deviance(null), df.residual(null), lower.tail = FALSE)
#> [1] 0.9192931
# quasibinomial keeps the mean-variance relation and estimates a dispersion.
quasi_null <- glm(cbind(late, on_time) ~ 1,
family = quasibinomial, data = counts
)
round(summary(quasi_null)$dispersion, 3)
#> [1] 0.287
round(sum(residuals(null, type = "pearson")^2) / df.residual(null), 3)
#> [1] 0.287Source: R documentation: glm
The six-parameter model on six groups is saturated (zero residual degrees of freedom, zero deviance), and the model with only an intercept has a residual deviance of 1.45 on 5 degrees of freedom, a goodness-of-fit p-value of 0.92: the groups do not differ in their overdue rate. The deviance of the same model on the 50 individual rows (67.8) is not zero or comparable, because the deviance of ungrouped binary data cannot be used as a goodness-of-fit statistic.
Family objects
A family object is a list of functions and expressions that together define the model and the fitting: the link function and its inverse, the variance function, the deviance residuals and the AIC. Reading them makes the algebra of the table above concrete. fam$linkfun maps a mean to the linear predictor and fam$linkinv maps back, fam$mu.eta is the derivative of the inverse link that IRLS needs, fam$variance is V(mu), and fam$dev.resids the squared deviance residuals. make.link() builds the link part alone. A string that is not a known link is refused with ‘nonsense’ link not recognised.
# A family generator is a function; calling it gives a family object.
fam <- binomial()
class(fam)
#> [1] "family"
fam$family
#> [1] "binomial"
fam$link
#> [1] "logit"
names(fam)
#> [1] "family" "link" "linkfun" "linkinv" "variance"
#> [6] "dev.resids" "aic" "mu.eta" "initialize" "validmu"
#> [11] "valideta" "simulate" "dispersion"
# A different link is chosen with link = (a string or a name).
probit <- binomial(link = "probit")
probit$link
#> [1] "probit"
tryCatch(binomial(link = "nonsense"), error = \(e) conditionMessage(e))
#> [1] "‘nonsense’ link not recognised"
# linkfun() maps a mean to the linear predictor, linkinv() maps it back.
p <- c(0.1, 0.5, 0.9)
fam$linkfun(p)
#> [1] -2.197225 0.000000 2.197225
all.equal(fam$linkfun(p), qlogis(p))
#> [1] TRUE
fam$linkinv(fam$linkfun(p))
#> [1] 0.1 0.5 0.9
probit$linkfun(p)
#> [1] -1.281552 0.000000 1.281552
all.equal(probit$linkfun(p), qnorm(p))
#> [1] TRUE
# mu.eta() is the derivative of linkinv(), used by the fitting algorithm.
eta <- c(-1, 0, 1)
fam$mu.eta(eta)
#> [1] 0.1966119 0.2500000 0.1966119
all.equal(fam$mu.eta(eta), dlogis(eta))
#> [1] TRUE
# variance() is the variance function V(mu): mu (1 - mu) for the binomial.
fam$variance(p)
#> [1] 0.09 0.25 0.09
poisson()$variance(c(1, 4, 9))
#> [1] 1 4 9
Gamma()$variance(c(1, 2, 3))
#> [1] 1 4 9
gaussian()$variance(c(1, 2, 3))
#> [1] 1 1 1
# dev.resids() gives the squared deviance residuals: 2 * log(1 / 0.8) here.
fam$dev.resids(y = 1, mu = 0.8, wt = 1)
#> [1] 0.4462871
-2 * log(0.8)
#> [1] 0.4462871
# The links each family offers (R-intro 11.6.1), tried by calling the generator.
vapply(
c("logit", "probit", "cloglog", "cauchit", "log"),
\(l) binomial(link = l)$link, character(1)
)
#> logit probit cloglog cauchit log
#> "logit" "probit" "cloglog" "cauchit" "log"
poisson(link = "sqrt")$link
#> [1] "sqrt"
quasi(link = "log", variance = "mu")$family
#> [1] "quasi"
# family() extracts the family from a fitted model.
fit <- glm(carb ~ wt, family = poisson, data = mtcars)
family(fit)
#>
#> Family: poisson
#> Link function: log
#>Source: R documentation: family
For the binomial family the variance is mu * (1 - mu), for the Poisson it equals mu, for the Gamma it is mu squared and for the gaussian family it is constant. The family() generic extracts the family from a fitted model, which is how predict() finds the inverse link.
The glm object
A glm fit is an lm fit with more components, and its class vector says so: class(gfit) is c("glm", "lm"), so any generic with no glm method falls back to the lm method (hatvalues() is an example), and inherits(gfit, "lm") is true. The extra components include the family, the linear predictors, the IRLS weights and iteration count, the null deviance and degrees of freedom, and the response y. summary(gfit) is a summary.glm, with its own print method.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
dur$overdue <- dur$days > loan_days
gfit <- glm(overdue ~ pages + genre, family = binomial, data = dur)
lfit <- lm(days ~ pages + genre, data = dur)
# A glm is an lm in the S3 sense: class c("glm", "lm").
class(gfit)
#> [1] "glm" "lm"
inherits(gfit, "lm")
#> [1] TRUE
inherits(lfit, "glm")
#> [1] FALSE
inherits(gfit, c("lm", "glm"), which = TRUE)
#> [1] 2 1
# The components of the list.
names(gfit)
#> [1] "coefficients" "residuals" "fitted.values"
#> [4] "effects" "R" "rank"
#> [7] "qr" "family" "linear.predictors"
#> [10] "deviance" "aic" "null.deviance"
#> [13] "iter" "weights" "prior.weights"
#> [16] "df.residual" "df.null" "y"
#> [19] "converged" "boundary" "model"
#> [22] "na.action" "call" "formula"
#> [25] "terms" "data" "offset"
#> [28] "control" "method" "contrasts"
#> [31] "xlevels"
# The parts of the IRLS fit: iterations, weights, working response.
gfit$iter
#> [1] 4
gfit$converged
#> [1] TRUE
round(head(gfit$weights, 3), 4)
#> 1 2 3
#> 0.2343 0.2343 0.2343
round(head(gfit$linear.predictors, 3), 3)
#> 1 2 3
#> -0.512 -0.512 -0.512
head(gfit$y, 3)
#> 1 2 3
#> 0 0 1
gfit$df.null
#> [1] 41
gfit$family$family
#> [1] "binomial"
# Methods found for "glm" first, then for "lm": summary() is the glm method...
class(summary(gfit))
#> [1] "summary.glm"
# ...while hatvalues() has no glm method of its own and uses the lm one.
round(head(hatvalues(gfit), 3), 4)
#> 1 2 3
#> 0.077 0.077 0.077
is.null(getS3method("hatvalues", "glm", optional = TRUE))
#> [1] TRUESource: R documentation: glm
Nonlinear least squares with nls()
Some models are linear in their predictors and some are not. The manual separates the two cases: certain nonlinear models can be fitted by a GLM, but in the majority of cases the nonlinear curve-fitting problem is one of nonlinear optimisation, where there is no guarantee that the procedure converges on satisfactory estimates and every method needs initial guesses for the parameters (An Introduction to R: Nonlinear least squares and maximum likelihood models). nls() minimises the sum of squared residuals of a formula whose right-hand side is an arbitrary R expression in the data and the parameters, by Gauss-Newton iteration from the start you supply. Parameters are the names in start; everything else in the formula is data.
The loan durations give a natural curve to fit. The share of returned loans that were back after at most d days is the empirical cumulative distribution function of the durations, ecdf(), which rises from near zero to one and looks like a logistic curve. The logistic curve Asym / (1 + expxmid - d) / scal has an asymptote Asym, a midpoint xmid (the duration at which half of the asymptote is reached) and a scale scal. The points of a cumulative curve are not independent, so the standard errors below are optimistic; the example is about fitting a curve, not about inference.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
# Share of returned loans that were back after at most d days: the empirical
# cumulative distribution of the loan durations, at d = 1, 2, ..., 40.
cdf <- ecdf(dur$days)
curve_data <- data.frame(d = 1:40, share = cdf(1:40))
head(curve_data, 5)
#> d share
#> 1 1 0.02
#> 2 2 0.02
#> 3 3 0.02
#> 4 4 0.06
#> 5 5 0.06
tail(curve_data, 3)
#> d share
#> 38 38 1
#> 39 39 1
#> 40 40 1
# A logistic curve: share = Asym / (1 + exp((xmid - d) / scal)).
# nls() needs a starting guess for every parameter.
nfit <- nls(share ~ Asym / (1 + exp((xmid - d) / scal)),
data = curve_data,
start = list(Asym = 1, xmid = 20, scal = 5)
)
nfit
#> Nonlinear regression model
#> model: share ~ Asym/(1 + exp((xmid - d)/scal))
#> data: curve_data
#> Asym xmid scal
#> 1.06 20.93 5.77
#> residual sum-of-squares: 0.03216
#>
#> Number of iterations to convergence: 4
#> Achieved convergence tolerance: 6.411e-06
summary(nfit)
#>
#> Formula: share ~ Asym/(1 + exp((xmid - d)/scal))
#>
#> Parameters:
#> Estimate Std. Error t value Pr(>|t|)
#> Asym 1.05959 0.01761 60.17 <2e-16 ***
#> xmid 20.93045 0.31684 66.06 <2e-16 ***
#> scal 5.77012 0.22845 25.26 <2e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Residual standard error: 0.02948 on 37 degrees of freedom
#>
#> Number of iterations to convergence: 4
#> Achieved convergence tolerance: 6.411e-06
#>Source: R documentation: nls
The printed object lists the model, the data and the estimates, with the residual sum of squares and the iteration report: four Gauss-Newton iterations to a convergence tolerance of 6.4e-06, below the default tol of 1e-05. The fit says that half of the loans are back by day 20.9 (xmid), that the curve rises over a scale of about 5.8 days, and that the asymptote is estimated at 1.06, slightly above the true ceiling of one. summary() adds standard errors, t values and the residual standard error with its degrees of freedom (40 points minus 3 parameters).
The extractor generics work as for lm(): coef, vcov, residuals, fitted, predict, deviance, df.residual, logLik, AIC and sigma. R 4.6.0 added a hatvalues() method for nls fits. confint() for an nls fit is a profile-likelihood interval (R 4.4.0 moved the glm and nls methods from MASS into stats, so the MASS namespace is no longer loaded) and prints Waiting for profiling to be done...; a model with fewer parameters, here with the asymptote fixed at one, is compared with anova(small, big). Bounds on the parameters need algorithm = "port", which uses the PORT library, as nlminb() does.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
cdf <- ecdf(dur$days)
curve_data <- data.frame(d = 1:40, share = cdf(1:40))
nfit <- nls(share ~ Asym / (1 + exp((xmid - d) / scal)),
data = curve_data,
start = list(Asym = 1, xmid = 20, scal = 5)
)
# The extractor generics work as for lm().
round(coef(nfit), 3)
#> Asym xmid scal
#> 1.06 20.93 5.77
round(vcov(nfit), 5)
#> Asym xmid scal
#> Asym 0.00031 0.00477 0.00305
#> xmid 0.00477 0.10039 0.04692
#> scal 0.00305 0.04692 0.05219
df.residual(nfit)
#> [1] 37
deviance(nfit)
#> [1] 0.03216431
sigma(nfit)
#> [1] 0.02948399
logLik(nfit)
#> 'log Lik.' 85.75801 (df=4)
AIC(nfit)
#> [1] -163.516
round(head(fitted(nfit), 3), 4)
#> [1] 0.0325 0.0384 0.0454
round(head(residuals(nfit), 3), 4)
#> [1] -0.0125 -0.0184 -0.0254
predict(nfit, newdata = data.frame(d = c(14, 21, 28))) |> round(3)
#> [1] 0.245 0.533 0.819
# confint() profiles the likelihood: it prints a message while it works.
msgs <- character()
ci <- withCallingHandlers(confint(nfit), message = \(m) {
msgs <<- c(msgs, conditionMessage(m))
invokeRestart("muffleMessage")
})
msgs
#> [1] "Waiting for profiling to be done...\n"
round(ci, 3)
#> 2.5% 97.5%
#> Asym 1.027 1.096
#> xmid 20.338 21.573
#> scal 5.340 6.237
# Fix the asymptote at 1 by dropping the parameter: a nested model, so F test.
nfit1 <- nls(share ~ 1 / (1 + exp((xmid - d) / scal)),
data = curve_data,
start = list(xmid = 20, scal = 5)
)
round(coef(nfit1), 3)
#> xmid scal
#> 20.055 5.224
anova(nfit1, nfit)
#> Analysis of Variance Table
#>
#> Model 1: share ~ 1/(1 + exp((xmid - d)/scal))
#> Model 2: share ~ Asym/(1 + exp((xmid - d)/scal))
#> Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
#> 1 38 0.045203
#> 2 37 0.032164 1 0.013039 14.999 0.0004231 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# Bounds need algorithm = "port": here Asym may not exceed 1.
bounded <- nls(share ~ Asym / (1 + exp((xmid - d) / scal)),
data = curve_data,
start = list(Asym = 0.9, xmid = 20, scal = 5), algorithm = "port",
upper = c(Asym = 1, xmid = 100, scal = 100)
)
round(coef(bounded), 3)
#> Asym xmid scal
#> 1.000 20.055 5.224
bounded$convInfo$stopMessage
#> [1] "both X-convergence and relative convergence (5)"
# The convergence settings of a fit and the defaults.
nfit$convInfo$isConv
#> [1] TRUE
nfit$convInfo$finIter
#> [1] 4
str(nls.control())
#> List of 7
#> $ maxiter : num 50
#> $ tol : num 1e-05
#> $ minFactor : num 0.000977
#> $ printEval : logi FALSE
#> $ warnOnly : logi FALSE
#> $ scaleOffset: num 0
#> $ nDcentral : logi FALSESource: R documentation: predict.nls
The F test of the fixed-asymptote model against the free one gives F = 15.0 and p = 0.0004: the data prefer an asymptote slightly above one, which is a symptom of the curve being a poor description of a distribution that is closer to uniform. The bounded fit, with Asym <= 1, converges at the bound (Asym = 1) with the same xmid and scal as the model that fixed it. str(nls.control()) lists the convergence settings: maxiter = 50, tol = 1e-05, minFactor = 1/1024 (the smallest step-halving factor), printEval, warnOnly, scaleOffset and nDcentral. The last two arrived in R 4.1.0: scaleOffset adds a constant to the denominator of the relative-offset convergence criterion, for fits that are expected to be exact, and nDcentral = TRUE uses central differences for the numerical derivatives.
Self-starting models
Finding start values is the hard part of nls(). A self-starting model is a function with an attached routine that computes good starting values from the data, so start can be omitted. SSlogis(input, Asym, xmid, scal) is the logistic curve above; the others in stats are SSasymp, SSasympOff, SSasympOrig, SSbiexp, SSfol, SSfpl, SSgompertz, SSmicmen and SSweibull. The manual’s example fits the Michaelis-Menten model with nls(y ~ SSmicmen(x, Vm, K), df). getInitial() shows the start values that the model computed, and selfStart() builds a self-starting model from a function and an initial routine.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
cdf <- ecdf(dur$days)
curve_data <- data.frame(d = 1:40, share = cdf(1:40))
# SSlogis() is a self-starting model: its initial function finds the start.
sfit <- nls(share ~ SSlogis(d, Asym, xmid, scal), data = curve_data)
round(coef(sfit), 3)
#> Asym xmid scal
#> 1.06 20.93 5.77
sfit$convInfo$finIter
#> [1] 0
# It is the same curve as the hand-written formula, same parametrisation.
hand <- nls(share ~ Asym / (1 + exp((xmid - d) / scal)),
data = curve_data,
start = list(Asym = 1, xmid = 20, scal = 5)
)
all.equal(coef(sfit), coef(hand), tolerance = 1e-5)
#> [1] TRUE
# getInitial() shows the start that the model computed.
round(getInitial(share ~ SSlogis(d, Asym, xmid, scal), data = curve_data), 3)
#> Asym xmid scal
#> 1.06 20.93 5.77
# SSlogis is of class "selfStart": a function with an initial function.
class(SSlogis)
#> [1] "selfStart"
attr(SSlogis, "pnames")
#> [1] "Asym" "xmid" "scal"
round(SSlogis(c(10, 21, 30), Asym = 1, xmid = 21, scal = 5.8), 4)
#> [1] 0.1305 0.5000 0.8252
# The self-starting models that stats provides.
ss <- ls("package:stats", pattern = "^SS")
ss[vapply(ss, \(f) inherits(get(f, "package:stats"), "selfStart"), logical(1))]
#> [1] "SSasymp" "SSasympOff" "SSasympOrig" "SSbiexp" "SSfol"
#> [6] "SSfpl" "SSgompertz" "SSlogis" "SSmicmen" "SSweibull"Source: R documentation: SSlogis
SSlogis reaches the same estimates as the hand-written formula, and needs no Gauss-Newton iterations of its own (finIter is 0): the routine that computes the start values has already found the least-squares solution. R 4.1.0 changed the interface a little: nls() now passes trace and control to getInitial(), so the initial function of a self-starting model must have ... in its argument list.
Why nls() fails, and the messages it gives
nls() fails in a small number of recognisable ways, and they are worth knowing by their texts. singular gradient matrix at initial parameter estimates means the matrix of derivatives of the model with respect to the parameters has no full rank at the start values: either the model has a redundant parameter (case 1 below: the two slopes b and c cannot be told apart) or the start is so far from the data that the model is flat in some direction. singular gradient (without the rest) is the same problem met after some iterations. step factor ... reduced below 'minFactor' of ... means that Gauss-Newton proposed a step that did not reduce the residual sum of squares, and halving it down to minFactor still did not; it points to a bad start or a model that does not fit. number of iterations exceeded maximum of 50 is the iteration limit; it also happens for artificial data with zero residuals, where the convergence criterion (a ratio of the size of the residual component to the size of the rest) compares two vanishing numbers and never becomes small enough, and the help page warns against using nls() on such artificial data; scaleOffset (R 4.1.0) is the supported workaround. Missing value or an infinity produced when evaluating the model is a model that is not finite at the current parameters (case 4: dividing by scal = 0). Changing control = nls.control(warnOnly = TRUE) turns the convergence failures into warnings and returns the last iterate, which you can inspect but should not trust.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
cdf <- ecdf(dur$days)
curve_data <- data.frame(d = 1:40, share = cdf(1:40))
logistic <- share ~ Asym / (1 + exp((xmid - d) / scal))
msg <- \(expr) tryCatch(expr, error = \(e) conditionMessage(e))
# 1. A parameter that is not identifiable: a + b * d + c * d has two slopes.
msg(nls(share ~ a + b * d + c * d,
data = curve_data,
start = list(a = 0, b = 0.02, c = 0.01)
))
#> [1] "singular gradient matrix at initial parameter estimates"
# 2. A start so far from the data that the gradient is singular on the way.
msg(nls(logistic,
data = curve_data, start = list(Asym = 1, xmid = 40, scal = 10)
))
#> [1] "singular gradient"
# 3. A start from which the Gauss-Newton step keeps overshooting.
msg(nls(logistic,
data = curve_data, start = list(Asym = 0.3, xmid = 73, scal = 9)
))
#> [1] "step factor 0.000488281 reduced below 'minFactor' of 0.000976562"
# 4. A start where the model is not finite (scal = 0 divides by zero).
msg(nls(logistic,
data = curve_data, start = list(Asym = 1, xmid = 20, scal = 0)
))
#> [1] "Missing value or an infinity produced when evaluating the model"
# 5. Artificial data with zero residuals: the convergence test never triggers.
exact <- data.frame(d = 1:10, y = 2 * exp(0.3 * 1:10))
msg(nls(y ~ a * exp(b * d), data = exact, start = list(a = 1, b = 0.2)))
#> [1] "number of iterations exceeded maximum of 50"
# scaleOffset adds a constant to the denominator of the convergence criterion.
zero_ok <- nls(y ~ a * exp(b * d),
data = exact, start = list(a = 1, b = 0.2),
control = nls.control(scaleOffset = 1)
)
round(coef(zero_ok), 4)
#> a b
#> 2.0 0.3
# 6. Too few iterations: maxiter is 50 by default.
msg(nls(logistic,
data = curve_data, start = list(Asym = 1, xmid = 20, scal = 5),
control = nls.control(maxiter = 2)
))
#> [1] "number of iterations exceeded maximum of 2"
# warnOnly = TRUE turns the error into a warning and returns the last iterate.
tryCatch(
nls(logistic,
data = curve_data,
start = list(Asym = 1, xmid = 20, scal = 5),
control = nls.control(maxiter = 2, warnOnly = TRUE)
),
warning = \(w) conditionMessage(w)
)
#> [1] "number of iterations exceeded maximum of 2"Source: R documentation: nls.control
In practice: plot the data and the curve at the start values before fitting (the manual’s advice), rescale the predictors so the parameters are of similar magnitude, use a self-starting model if one exists, try algorithm = "port" with bounds when parameters must stay positive, and compare with a linear transformation if one exists. Community packages such as minpack.lm offer a Levenberg-Marquardt alternative that is more forgiving of poor starts; it is named here and not used. R-devel (the upcoming R 4.7.0, not run here) lets the named lower and upper bounds of nls() be given for only the parameters you want to constrain and in any order.
Optimisation and maximum likelihood
nls() is least squares specialised to a model formula. The general tools take any R function of the parameters. The manual names optim(), nlm() and nlminb() as the nonlinear optimisation routines, and shows two uses: minimising a sum of squared errors, and maximising a log-likelihood, the latter being a method that applies "even if the errors are not normal" (An Introduction to R: Maximum likelihood). The function to minimise takes the parameter vector as its first argument; to maximise, return the negative, or set control = list(fnscale = -1) for optim().
optim() and the negative log-likelihood
The loan durations (50 values, mean 20.2 days, standard deviation 8.6) give a likelihood to work with. Assume for a moment that they are normal: the negative log-likelihood is minus the sum of dnorm(x, mean, sd, log = TRUE). Working with log = TRUE avoids the underflow of multiplying many small densities, and the standard deviation is parametrised on the log scale so that the optimiser cannot propose a negative value. optim() offers six methods: "Nelder-Mead" (the default, derivative-free and robust, but slow), "BFGS" (a quasi-Newton method using gradients, by finite differences if you give none), "CG" (conjugate gradients), "L-BFGS-B" (BFGS with box constraints, the only method that accepts lower and upper), "SANN" (simulated annealing) and "Brent" (one dimension). The result is a list with the best parameters par, the value of the function there, the counts of function and gradient evaluations, a convergence code (0 is success) and a message. With hessian = TRUE the Hessian at the solution is added.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
days <- dur$days
c(n = length(days), mean = mean(days), sd = sd(days))
#> n mean sd
#> 50.000000 20.220000 8.633797
# The negative log-likelihood of a normal sample; sigma is optimised on the log
# scale so that it stays positive.
nll_normal <- function(par, x) {
-sum(dnorm(x, mean = par[1], sd = exp(par[2]), log = TRUE))
}
# Nelder-Mead (the default): derivative-free, needs only the function.
nm <- optim(c(mu = 10, log_sigma = 1), nll_normal, x = days)
names(nm)
#> [1] "par" "value" "counts" "convergence" "message"
round(nm$par, 3)
#> mu log_sigma
#> 20.216 2.146
nm$convergence
#> [1] 0
nm$value
#> [1] 178.2261
nm$counts
#> function gradient
#> 59 NA
# Closed form: the mean, and the standard deviation with divisor n.
c(mu = mean(days), log_sigma = log(sqrt(mean((days - mean(days))^2))))
#> mu log_sigma
#> 20.220000 2.145583
# BFGS uses gradients (here by finite differences) and can return the Hessian.
bfgs <- optim(c(mu = 10, log_sigma = 1), nll_normal,
x = days,
method = "BFGS", hessian = TRUE
)
round(bfgs$par, 4)
#> mu log_sigma
#> 20.2200 2.1456
bfgs$counts
#> function gradient
#> 61 17
round(bfgs$hessian, 3)
#> mu log_sigma
#> mu 0.684 0
#> log_sigma 0.000 100
# The inverse of the Hessian of the negative log-likelihood is the
# variance-covariance matrix of the estimates; the square roots of its
# diagonal are the standard errors.
se <- sqrt(diag(solve(bfgs$hessian)))
round(se, 4)
#> mu log_sigma
#> 1.2087 0.1000
# Check: the standard error of the mean is sd / sqrt(n) (with divisor n here).
sqrt(mean((days - mean(days))^2) / length(days))
#> [1] 1.208732Source: R documentation: optim
The estimates agree with the closed forms (the mean, and the standard deviation with divisor n instead of n - 1) to within 0.005 for Nelder-Mead and to four decimals for BFGS. The Hessian of the negative log-likelihood is the observed information; its inverse is the estimated covariance matrix of the maximum-likelihood estimates, and the square roots of its diagonal are the standard errors: 1.21 for the mean, which equals the standard deviation over the square root of n (the check at the end), and 0.10 for the log standard deviation.
Bounds, maximising and gradients
Real parameters often have a natural range. For a Gamma model of the durations, shape and rate must be positive, so "L-BFGS-B" with lower bounds is the right method. Without bounds the unconstrained methods may step to a negative value, where dgamma() returns NaN with the warning NaNs produced; the script catches it. The maximum-likelihood fit is not the same as the method of moments (shape = mean^2 / variance), as the comparison shows: 3.57 against 5.48, because the durations are close to uniform between 1 and 37 days, not Gamma shaped. An analytic gradient given as gr makes the search quicker and more accurate; here the score of the Gamma log-likelihood is written with digamma(). control takes maxit (the iteration limit), reltol (the relative tolerance, about 1.5e-08), parscale (the typical size of each parameter) and fnscale; a convergence of 1 says that the limit of iterations was reached, 10 degeneracy of the Nelder-Mead simplex, and 51 or 52 a warning or an error from L-BFGS-B.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
days <- dur$days
# A gamma model of the durations: both parameters must be positive, so use
# L-BFGS-B, the only optim() method that accepts lower and upper bounds.
nll_gamma <- function(par, x) {
-sum(dgamma(x, shape = par[1], rate = par[2], log = TRUE))
}
bounded <- optim(c(shape = 1, rate = 0.1), nll_gamma,
x = days,
method = "L-BFGS-B", lower = c(1e-4, 1e-4), upper = c(Inf, Inf)
)
round(bounded$par, 4)
#> shape rate
#> 3.5673 0.1764
bounded$convergence
#> [1] 0
bounded$message
#> [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
# Method of moments for comparison: shape = mean^2 / var, rate = mean / var.
m <- mean(days)
v <- var(days)
round(c(shape = m^2 / v, rate = m / v), 4)
#> shape rate
#> 5.4848 0.2713
# Without bounds, Nelder-Mead can wander into negative parameters: dgamma()
# then returns NaN, which optim() reports as a warning.
caught <- tryCatch(
optim(c(shape = 1, rate = 0.1), nll_gamma, x = days),
warning = \(w) conditionMessage(w)
)
caught
#> [1] "NaNs produced"
# To maximise, either negate the function or set control$fnscale to -1.
loglik_gamma <- function(par, x) {
sum(dgamma(x, shape = par[1], rate = par[2], log = TRUE))
}
maxi <- optim(c(shape = 1, rate = 0.1), loglik_gamma,
x = days,
method = "L-BFGS-B", lower = c(1e-4, 1e-4), control = list(fnscale = -1)
)
all.equal(maxi$par, bounded$par, tolerance = 1e-4)
#> [1] TRUE
all.equal(maxi$value, -bounded$value, tolerance = 1e-6)
#> [1] TRUE
# A gradient (gr =) makes the search faster and more accurate: here the
# analytic score.
score_nll <- function(par, x) {
n <- length(x)
c(
-(n * (log(par[2]) - digamma(par[1])) + sum(log(x))),
-(n * par[1] / par[2] - sum(x))
)
}
with_grad <- optim(c(shape = 2, rate = 0.1), nll_gamma,
gr = score_nll, x = days,
method = "L-BFGS-B", lower = c(1e-4, 1e-4)
)
round(with_grad$par, 4)
#> shape rate
#> 3.5671 0.1764
with_grad$counts
#> function gradient
#> 12 12
# control: maxit stops early, and convergence = 1 says so (0 means success).
early <- optim(c(shape = 3, rate = 0.15), nll_gamma,
x = days,
method = "L-BFGS-B", lower = c(1e-4, 1e-4), control = list(maxit = 2)
)
early$convergence
#> [1] 1
early$message
#> [1] "NEW_X"
early$counts
#> function gradient
#> 6 6Source: R documentation: optim
One-dimensional problems: optimize(), uniroot(), nlm() and nlminb()
optimize(f, interval) (alias optimise) minimises a function of one variable on an interval by golden-section search combined with parabolic interpolation, and maximum = TRUE maximises. Its default tolerance is about 1.2e-04, so the result of the exponential rate (0.04947 against the exact 1 / mean of 0.04946) is accurate to four digits unless tol is tightened. uniroot(f, interval) finds a zero of a function of one variable (Brent’s method) and needs the function to change sign across the interval; otherwise it stops with f() values at end points not of opposite sign. Solving the score equation of the exponential model gives the same rate. nlm() is a Newton-type minimiser that returns the estimate, the gradient, optionally the Hessian, and a code (1 means the gradient is close to zero); nlminb() uses the PORT library, accepts bounds, and returns the parameters, the objective, a convergence code and a message. Its relative-convergence message is of the same kind as the one nls(algorithm = "port") reports. The R-devel NEWS (in its section for the scheduled R 4.6.2, not run here) lists a fix of a long-standing bug in nlminb() with constraints.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
days <- dur$days
# optimize() minimises a function of one variable on an interval.
nll_exp <- function(rate) -sum(dexp(days, rate = rate, log = TRUE))
op <- optimize(nll_exp, interval = c(0.001, 1))
names(op)
#> [1] "minimum" "objective"
round(unlist(op), 5)
#> minimum objective
#> 0.04947 200.33361
1 / mean(days)
#> [1] 0.04945598
# The default tolerance is about 1.2e-4; tighten it for more digits.
round(optimize(nll_exp, interval = c(0.001, 1), tol = 1e-10)$minimum, 6)
#> [1] 0.049456
# maximum = TRUE maximises instead.
mx <- optimize(\(rate) -nll_exp(rate), interval = c(0.001, 1), maximum = TRUE)
round(mx$maximum, 5)
#> [1] 0.04947
# uniroot() finds a zero of a function of one variable: where the score is 0.
score <- \(rate) length(days) / rate - sum(days)
root <- uniroot(score, interval = c(0.001, 1), tol = 1e-10)
round(root$root, 6)
#> [1] 0.049456
names(root)
#> [1] "root" "f.root" "iter" "init.it" "estim.prec"
tryCatch(uniroot(score, interval = c(0.1, 1)), error = \(e) conditionMessage(e))
#> [1] "f() values at end points not of opposite sign"
# nlm(): a Newton-type minimiser; its Hessian is optional.
nll_normal <- function(par) -sum(dnorm(days, par[1], exp(par[2]), log = TRUE))
nl <- nlm(nll_normal, p = c(15, 2), hessian = TRUE)
names(nl)
#> [1] "minimum" "estimate" "gradient" "hessian" "code"
#> [6] "iterations"
round(nl$estimate, 4)
#> [1] 20.2200 2.1456
nl$code
#> [1] 1
round(sqrt(diag(solve(nl$hessian))), 4)
#> [1] 1.2087 0.1000
# nlminb(): the PORT library; it accepts bounds, like nls(algorithm = "port").
nb <- nlminb(c(10, 1), nll_normal, lower = c(-Inf, -5), upper = c(Inf, 5))
names(nb)
#> [1] "par" "objective" "convergence" "iterations" "evaluations"
#> [6] "message"
round(nb$par, 4)
#> [1] 20.2200 2.1456
nb$convergence
#> [1] 0
nb$message
#> [1] "relative convergence (4)"Source: R documentation: optimize
The other functions of the script have their own reference pages: R documentation: uniroot, R documentation: nlm and R documentation: nlminb.
All four methods find the same estimates of the normal model, 20.22 and 2.1456 for the log standard deviation, and the Hessian from nlm() reproduces the standard errors from optim(). A start that sends exp(par[2]) to infinity makes nlm() warn Inf replaced by maximum positive value; the start values in the script avoid it, which is a reminder that the function you give an optimiser must be finite wherever the optimiser may go.
Maximum likelihood with stats4::mle()
Writing optim() calls, extracting the Hessian and inverting it is the same boilerplate for every likelihood. The stats4 package (part of R) wraps it in mle(minuslogl, start, method, fixed, nobs, lower, upper, ...): you give the negative log-likelihood as a function of named parameters and a list of starting values, and get an object of S4 class "mle" with coef(), vcov(), logLik(), AIC(), BIC() (given nobs), summary(), profile() and confint() methods. method defaults to "BFGS", or to "L-BFGS-B" when bounds are given. The function’s parameters are its formal arguments, and the data are found in the calling environment. Because the object is S4, the printing uses show() and the slots are @ fields (see S4 Classes and Methods).
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
library(stats4)
days <- dur$days
# mle() takes the negative log-likelihood as a function of the parameters by
# name, and a list of starting values. The data is found in the global
# environment.
nll <- function(shape = 2, rate = 0.1) {
-sum(dgamma(days, shape = shape, rate = rate, log = TRUE))
}
fit <- mle(nll,
start = list(shape = 2, rate = 0.1), method = "L-BFGS-B",
lower = c(1e-4, 1e-4), nobs = length(days)
)
class(fit)
#> [1] "mle"
#> attr(,"package")
#> [1] "stats4"
isS4(fit)
#> [1] TRUE
fit
#>
#> Call:
#> mle(minuslogl = nll, start = list(shape = 2, rate = 0.1), method = "L-BFGS-B",
#> nobs = length(days), lower = c(1e-04, 1e-04))
#>
#> Coefficients:
#> shape rate
#> 3.5672911 0.1764259
summary(fit)
#> Maximum likelihood estimation
#>
#> Call:
#> mle(minuslogl = nll, start = list(shape = 2, rate = 0.1), method = "L-BFGS-B",
#> nobs = length(days), lower = c(1e-04, 1e-04))
#>
#> Coefficients:
#> Estimate Std. Error
#> shape 3.5672911 0.68263235
#> rate 0.1764259 0.03625128
#>
#> -2 log L: 368.9553
# The accessors are S4 methods of the generics from stats.
round(coef(fit), 4)
#> shape rate
#> 3.5673 0.1764
round(vcov(fit), 5)
#> shape rate
#> shape 0.46599 0.02304
#> rate 0.02304 0.00131
round(sqrt(diag(vcov(fit))), 4)
#> shape rate
#> 0.6826 0.0363
logLik(fit)
#> 'log Lik.' -184.4777 (df=2)
AIC(fit)
#> [1] 372.9553
BIC(fit)
#> [1] 376.7794
nobs(fit)
#> [1] 50Source: R documentation: mle
The Gamma model of the durations has shape 3.57 (standard error 0.68) and rate 0.176 (0.036), and -2 log L of 368.96, which logLik() reports as -184.48 with 2 degrees of freedom. The estimates equal those of optim() in the previous section, as they must: mle() calls it. The R 4.6.0 NEWS says that AIC() and BIC() now check whether the methods namespace (not stats4) is loaded, still use stats4::logLik and stats4::nobs, so packages need not load stats4 explicitly; and that profiling "mle" objects, which could fail because the profile started from the full maximum-likelihood estimate, now starts from the profile estimate of the previous step.
confint() on an "mle" object is a profile likelihood interval: it re-optimises the other parameters at a grid of values of each parameter and finds where the likelihood ratio crosses the chi-squared cutoff. It prints Profiling... on the console. Unlike the Wald interval from vcov(), it is not symmetric, which is the better description of a parameter, like the shape of a Gamma, that has a skewed likelihood. fixed = holds a parameter at a value; since R 4.0.0 an argument of the likelihood may be a vector (the coefficients are then named p1 and p2), with bounds and starting values given for the whole vector.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
library(stats4)
days <- dur$days
nll <- function(shape = 2, rate = 0.1) {
-sum(dgamma(days, shape = shape, rate = rate, log = TRUE))
}
fit <- mle(nll,
start = list(shape = 2, rate = 0.1), method = "L-BFGS-B",
lower = c(1e-4, 1e-4)
)
# confint() on an mle object profiles the likelihood; it prints "Profiling..."
# on the console and returns a matrix.
ci <- confint(fit)
#> Profiling...
round(ci, 3)
#> 2.5 % 97.5 %
#> shape 2.397 5.085
#> rate 0.114 0.257
# Wald intervals from the covariance matrix are symmetric around the estimate.
est <- coef(fit)
se <- sqrt(diag(vcov(fit)))
z <- qnorm(0.975)
round(cbind(lower = est - z * se, upper = est + z * se), 3)
#> lower upper
#> shape 2.229 4.905
#> rate 0.105 0.247
# Fixing a parameter: the rate is fixed, so only the shape is estimated.
fit_fixed <- mle(nll,
start = list(shape = 2), fixed = list(rate = 0.1),
method = "BFGS"
)
round(coef(fit_fixed), 4)
#> shape rate
#> 2.2231 0.1000
# Since R 4.0.0 an argument of the likelihood may be a vector (the coefficients
# are then named p1, p2, ...); bounds are given for the whole vector.
nll_vec <- function(p = c(2, 0.1)) {
-sum(dgamma(days, shape = p[1], rate = p[2], log = TRUE))
}
fit_vec <- mle(nll_vec,
start = list(p = c(2, 0.1)), method = "L-BFGS-B",
lower = c(1e-4, 1e-4)
)
round(coef(fit_vec), 4)
#> p1 p2
#> 3.5673 0.1764The profile interval for the shape, 2.40 to 5.09, reaches 1.17 below the estimate of 3.57 and 1.52 above it, where the Wald interval (2.23 to 4.91) is symmetric, 1.34 each way; the rate interval, 0.114 to 0.257, is skewed in the same direction. With rate fixed at 0.1 the shape is estimated at 2.22.
Some non-standard models
The manual closes its chapter on models with a short list of facilities beyond lm() and glm() (An Introduction to R: Some non-standard models): mixed models (lme() and nlme() in the recommended nlme package, which make heavy use of formulas), local approximating regressions (loess(), in stats), robust regression (lqs() and rlm() in the recommended MASS package), additive models (mgcv, also recommended), and tree-based models. Survival models are in the recommended survival package. They share the interface of this page, a formula, a data argument and the same generics, so the habits carry over: the recommended packages ship with R and are installed with it, and are named here without being run.
The model object as an S3 list
A fitted lm is a list. class(fit) is "lm", is.list(fit) is true, and names(fit) lists its components: coefficients, residuals, effects, rank, fitted.values, assign, qr, df.residual, contrasts, xlevels, call, terms and model. The qr component is the QR decomposition of the design matrix, terms the terms object (a formula with attributes, one of which is the environment), and model the model frame. The script prints only classes, dimensions and names for terms and model, and uses str(fit, max.level = 1, give.attr = FALSE), which lists the components without descending: calling str() on fit$terms would print <environment: ...> with a memory address that changes every run.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
dur <- dur[complete.cases(dur[, c("days", "pages", "genre")]), ]
fit <- lm(days ~ pages + genre, data = dur)
# A fitted model is a plain list with a class attribute.
class(fit)
#> [1] "lm"
is.list(fit)
#> [1] TRUE
length(fit)
#> [1] 13
names(fit)
#> [1] "coefficients" "residuals" "effects" "rank"
#> [5] "fitted.values" "assign" "qr" "df.residual"
#> [9] "contrasts" "xlevels" "call" "terms"
#> [13] "model"
# Component by component. terms and model carry an environment, so only
# their classes and sizes are printed.
fit$rank
#> [1] 7
fit$df.residual
#> [1] 35
fit$call
#> lm(formula = days ~ pages + genre, data = dur)
class(fit$qr)
#> [1] "qr"
names(fit$qr)
#> [1] "qr" "qraux" "pivot" "tol" "rank"
dim(fit$qr$qr)
#> [1] 42 7
class(fit$terms)
#> [1] "terms" "formula"
class(fit$model)
#> [1] "data.frame"
dim(fit$model)
#> [1] 42 3
fit$contrasts
#> $genre
#> [1] "contr.treatment"
#>
# str() with max.level = 1 lists the components without descending into them.
str(fit, max.level = 1, give.attr = FALSE)
#> List of 13
#> $ coefficients : Named num [1:7] 21.054374 0.000931 -3.788224 -0.544832 0.282682 ...
#> $ residuals : Named num [1:42] -7.47 1.53 13.53 7.53 5.53 ...
#> $ effects : Named num [1:42] -127.609 4.597 8.574 -0.263 -1.599 ...
#> $ rank : int 7
#> $ fitted.values: Named num [1:42] 17.5 17.5 17.5 17.5 17.5 ...
#> $ assign : int [1:7] 0 1 2 2 2 2 2
#> $ qr :List of 5
#> $ df.residual : int 35
#> $ contrasts :List of 1
#> $ xlevels :List of 1
#> $ call : language lm(formula = days ~ pages + genre, data = dur)
#> $ terms :Classes 'terms', 'formula' language days ~ pages + genre
#> $ model :'data.frame': 42 obs. of 3 variables:Source: R documentation: lm
The list is the representation; the generics are the interface. coef(fit) returns fit$coefficients for an lm, but an nls fit has no coefficients component at all (its parameters live in the m component, a closure), coef() still works on it, and $ on a list matches names partially, so fit$coef returns the coefficients too, which is a convenience at the console and a hazard in a package. Removing the class with unclass() exposes the plain list and turns every generic back into an error (no applicable method for 'predict' applied to an object of class "list"); restoring the class restores the behaviour. summary(fit) returns a different object of class "summary.lm", a list with call, terms, residuals, coefficients (a matrix), aliased, sigma, df, r.squared, adj.r.squared, fstatistic and cov.unscaled, and its own print method: what you saw as the "summary" is just that method running. To pull a number out of a model, take it from the summary object by name (summary(fit)$r.squared) rather than parsing the printed text.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
dur <- dur[complete.cases(dur[, c("days", "pages", "genre")]), ]
fit <- lm(days ~ pages + genre, data = dur)
# coef(fit) is the generic; fit$coefficients is the component it returns for lm.
identical(coef(fit), fit$coefficients)
#> [1] TRUE
# $ matches names partially, so fit$coef also works, but only by accident.
identical(fit$coef, fit$coefficients)
#> [1] TRUE
# A generic dispatches on the class, so it keeps working for other model
# classes, where the component may have another name or be computed.
nfit <- nls(days ~ a + b * pages, data = dur, start = list(a = 10, b = 0.01))
names(nfit)
#> [1] "m" "convInfo" "data" "call" "dataClasses"
#> [6] "control"
is.null(nfit$coefficients)
#> [1] TRUE
round(coef(nfit), 4)
#> a b
#> 17.8693 0.0068
# Removing the class exposes the list; the generics no longer apply.
raw <- unclass(fit)
class(raw)
#> [1] "list"
is.null(raw$coefficients)
#> [1] FALSE
tryCatch(predict(raw), error = \(e) conditionMessage(e))
#> [1] "no applicable method for 'predict' applied to an object of class \"list\""
inherits(raw, "lm")
#> [1] FALSE
# Setting the class back restores the behaviour.
class(raw) <- "lm"
all.equal(fitted(raw), fitted(fit))
#> [1] TRUE
# summary() returns another S3 object, with its own print method.
s <- summary(fit)
class(s)
#> [1] "summary.lm"
names(s)
#> [1] "call" "terms" "residuals" "coefficients"
#> [5] "aliased" "sigma" "df" "r.squared"
#> [9] "adj.r.squared" "fstatistic" "cov.unscaled"
round(s$r.squared, 4)
#> [1] 0.0325
round(s$adj.r.squared, 4)
#> [1] -0.1333
round(s$fstatistic, 4)
#> value numdf dendf
#> 0.196 6.000 35.000
round(s$sigma, 4)
#> [1] 9.2255
s$df
#> [1] 7 35 7
class(s$coefficients)
#> [1] "matrix" "array"
round(s$coefficients["pages", ], 4)
#> Estimate Std. Error t value Pr(>|t|)
#> 0.0009 0.0417 0.0223 0.9823
# Printing it calls the print method for class "summary.lm", found in stats.
is.function(getS3method("print", "summary.lm"))
#> [1] TRUE
is.function(getS3method("print", "lm"))
#> [1] TRUESource: R documentation: getS3method
methods(class = "lm") lists the S3 methods registered for the class, plus the S4 methods that stats and methods register for it (coerce, initialize, show, slotsFromS3): 40 entries in the checked installation. It shows which generics have an lm method of their own. Notably coef has none (the default method reads the coefficients component), anova, confint, predict, residuals, summary, plot and vcov do, and a glm has its own methods and inherits the rest from lm through the class vector c("glm", "lm"). getS3method(generic, class) returns the method itself, and optional = TRUE returns NULL instead of an error when it does not exist, which is how the script shows that hatvalues has no glm method.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
dur <- dur[complete.cases(dur[, c("days", "pages", "genre")]), ]
fit <- lm(days ~ pages + genre, data = dur)
gfit <- glm(I(days > loan_days) ~ pages + genre, family = binomial, data = dur)
# methods() lists the S3 methods registered or visible for a class.
m <- methods(class = "lm")
class(m)
#> [1] "MethodsFunction"
m
#> [1] add1 alias anova case.names coerce
#> [6] confint cooks.distance deviance dfbeta dfbetas
#> [11] dffits drop1 dummy.coef effects extractAIC
#> [16] family formula hatvalues influence initialize
#> [21] kappa labels logLik model.frame model.matrix
#> [26] nobs plot predict print proj
#> [31] qr residuals rstandard rstudent show
#> [36] simulate slotsFromS3 summary variable.names vcov
#> see '?methods' for accessing help and source code
# Which of the generics used in this page have an lm method of their own?
wanted <- c(
"anova", "coef", "confint", "predict", "residuals", "summary", "plot",
"vcov"
)
setNames(wanted %in% format(m), wanted)
#> anova coef confint predict residuals summary plot vcov
#> TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE
# A glm has its own methods and inherits those of lm.
head(format(methods(class = "glm")), 8)
#> [1] "add1" "anova" "coerce" "confint"
#> [5] "cooks.distance" "deviance" "dfbetas" "dffits"
class(gfit)
#> [1] "glm" "lm"
# Which method will run? getS3method() finds it.
environmentName(environment(getS3method("summary", "lm")))
#> [1] "stats"
is.function(getS3method("residuals", "glm"))
#> [1] TRUE
is.null(getS3method("hatvalues", "glm", optional = TRUE))
#> [1] TRUESource: R documentation: methods
A wrapper class of your own
You can build on a model object in the two ways described on S3 Classes and Generics. A wrapper holds the fit in a component and defines methods for the generics it wants to support; a generic without a method falls through to its default, which for residuals() reads a residuals component that the wrapper does not have and returns NULL, with no warning. A subclass prepends a class to the fit’s own class vector, keeps all lm behaviour, and adds a method that calls NextMethod() to print the extra information first. The wrapper in this script prints a one-line summary; the subclass prints a header and then the ordinary lm print.
source("scenario.R")
# A small wrapper class around an lm fit: a constructor, a print method and
# methods for the generics that should keep working.
new_shelf_fit <- function(fit, label) {
stopifnot(inherits(fit, "lm"), is.character(label), length(label) == 1L)
structure(list(fit = fit, label = label), class = "shelf_fit")
}
print.shelf_fit <- function(x, ...) {
cat("<shelf_fit> ", x$label, "\n", sep = "")
cat(" formula:", deparse1(formula(x$fit)), "\n")
r2 <- round(summary(x$fit)$r.squared, 3)
cat(" n =", nobs(x$fit), " R-squared =", r2, "\n")
invisible(x)
}
coef.shelf_fit <- function(object, ...) coef(object$fit)
predict.shelf_fit <- function(object, ...) predict(object$fit, ...)
books$genre <- tolower(trimws(books$genre))
sf <- new_shelf_fit(lm(pages ~ year, data = books), "pages on year")
sf
#> <shelf_fit> pages on year
#> formula: pages ~ year
#> n = 13 R-squared = 0.257
coef(sf)
#> (Intercept) year
#> 2067.0926829 -0.8945014
predict(sf, newdata = data.frame(year = c(1990, 2020))) |> round(1)
#> 1 2
#> 287.0 260.2
# Without a method, a generic falls through to its default method, which looks
# for a $residuals component in the wrapper and finds none.
residuals(sf)
#> NULL
length(residuals(sf$fit))
#> [1] 13
# The other design: a subclass that keeps the lm behaviour and adds to it.
new_shelf_lm <- function(fit, label) {
stopifnot(inherits(fit, "lm"))
fit$label <- label
class(fit) <- c("shelf_lm", class(fit))
fit
}
print.shelf_lm <- function(x, ...) {
cat("<shelf_lm> ", x$label, "\n", sep = "")
NextMethod()
}
sl <- new_shelf_lm(lm(pages ~ year, data = books), "pages on year")
class(sl)
#> [1] "shelf_lm" "lm"
sl
#> <shelf_lm> pages on year
#>
#> Call:
#> lm(formula = pages ~ year, data = books)
#>
#> Coefficients:
#> (Intercept) year
#> 2067.0927 -0.8945
#>
round(head(residuals(sl), 3), 2)
#> 1 2 3
#> 85.65 -13.36 -47.09
inherits(sl, "lm")
#> [1] TRUEA tidy interface with broom
The broom package turns model objects into data frames, which is convenient when many models are compared or the results are plotted. tidy() returns one row per coefficient (term, estimate, standard error, statistic, p-value, and a confidence interval with conf.int = TRUE), glance() one row for the whole model (R-squared, sigma, F statistic, log-likelihood, AIC, BIC, deviance, degrees of freedom, nobs), and augment() the data with the fitted values, residuals, hat values, Cook’s distance and standardised residuals added as columns that start with a dot. tidy() and glance() have methods for glm, nls, htest and many other classes, so the output has the same shape across models; tidy(gfit, exponentiate = TRUE) returns odds ratios. The columns print as a tibble, a data frame class from the tidyverse (see The tidyverse and dplyr).
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
library(broom)
dur <- dur[complete.cases(dur[, c("days", "pages", "genre")]), ]
fit <- lm(days ~ pages + genre, data = dur)
# tidy(): one row per coefficient, as a tibble.
tidy(fit)
#> # A tibble: 7 × 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 21.1 16.0 1.31 0.197
#> 2 pages 0.000931 0.0417 0.0223 0.982
#> 3 genrefiction -3.79 8.29 -0.457 0.651
#> 4 genremystery -0.545 7.78 -0.0701 0.945
#> 5 genrenon-fiction 0.283 11.2 0.0251 0.980
#> 6 genrepoetry -1.14 13.0 -0.0878 0.931
#> 7 genrescience fiction -0.863 4.86 -0.178 0.860
tidy(fit, conf.int = TRUE)[1:2, ]
#> # A tibble: 2 × 7
#> term estimate std.error statistic p.value conf.low conf.high
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 21.1 16.0 1.31 0.197 -11.5 53.6
#> 2 pages 0.000931 0.0417 0.0223 0.982 -0.0837 0.0855
# glance(): one row for the whole model.
glance(fit)
#> # A tibble: 1 × 12
#> r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 0.0325 -0.133 9.23 0.196 0.976 6 -149. 314. 328.
#> # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>
# augment(): the data plus fitted values, residuals and influence measures.
augment(fit)[1:3, ]
#> # A tibble: 3 × 10
#> .rownames days pages genre .fitted .resid .hat .sigma .cooksd .std.resid
#> <chr> <dbl> <int> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 10 214 fiction 17.5 -7.47 0.0773 9.26 0.00849 -0.842
#> 2 2 19 214 fiction 17.5 1.53 0.0773 9.36 0.000359 0.173
#> 3 3 31 214 fiction 17.5 13.5 0.0773 9.04 0.0279 1.53
names(augment(fit))
#> [1] ".rownames" "days" "pages" "genre" ".fitted"
#> [6] ".resid" ".hat" ".sigma" ".cooksd" ".std.resid"
# The same functions work for a glm.
gfit <- glm(I(days > loan_days) ~ pages + genre, family = binomial, data = dur)
tidy(gfit, exponentiate = TRUE)[1:3, ]
#> # A tibble: 3 × 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 52.3 3.92 1.01 0.312
#> 2 pages 0.991 0.00986 -0.954 0.340
#> 3 genrefiction 0.0858 2.06 -1.19 0.234
glance(gfit)
#> # A tibble: 1 × 8
#> null.deviance df.null logLik AIC BIC deviance df.residual nobs
#> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <int>
#> 1 57.8 41 -27.7 69.5 81.6 55.5 35 42The broom reference has separate pages for glance() (broom package documentation: glance.lm) and augment() (broom package documentation: augment.lm). The base-R equivalents are coef(summary(fit)), c(summary(fit)$r.squared, AIC(fit)) and cbind(model.frame(fit), fitted(fit), residuals(fit)); broom’s value is the uniform shape, not new numbers.
Randomness: set.seed(), simulate() and the bootstrap
Anything that draws random numbers depends on the state of the random number generator, so set.seed() goes before it, in the same script (see Distributions and Random Numbers). simulate(fit, nsim) draws new responses from the fitted model (for an lm, the fitted values plus normal noise with the estimated residual variance) as a data frame with one column per simulation, named sim_1, sim_2 and so on. A nonparametric bootstrap resamples rows of the data with replacement and refits the model, and the spread of the refitted coefficients estimates their sampling variability without assuming normal errors; a parametric bootstrap refits on responses simulated from the model. With 13 books the bootstrap interval for the slope of pages on year, from -3.6 to 3.1, is much wider than the t interval of -1.9 to 0.1: the bootstrap exposes how little 13 points say, and a single resample can contain the same few books many times.
source("scenario.R")
books <- books[complete.cases(books[, c("year", "pages")]), ]
fit <- lm(pages ~ year, data = books)
# simulate() draws new responses from the fitted model; it uses the random
# number generator, so set the seed first.
set.seed(282)
sim <- simulate(fit, nsim = 3)
dim(sim)
#> [1] 13 3
names(sim)
#> [1] "sim_1" "sim_2" "sim_3"
round(head(sim, 3))
#> sim_1 sim_2 sim_3
#> 1 272 308 430
#> 2 405 257 380
#> 3 291 444 282
# The same seed reproduces the draw.
set.seed(282)
identical(simulate(fit, nsim = 3), sim)
#> [1] TRUE
# A nonparametric bootstrap of the slope: resample the rows, refit, collect.
set.seed(282)
slopes <- replicate(500, {
idx <- sample(nrow(books), replace = TRUE)
coef(lm(pages ~ year, data = books[idx, ]))[["year"]]
})
round(quantile(slopes, c(0.025, 0.5, 0.975)), 3)
#> 2.5% 50% 97.5%
#> -3.581 -0.931 3.092
# Compare with the t-based interval.
round(confint(fit)["year", ], 3)
#> 2.5 % 97.5 %
#> -1.904 0.115
# A parametric bootstrap from the model: simulate, refit, collect.
set.seed(282)
sim500 <- simulate(fit, nsim = 500)
slopes2 <- vapply(sim500, \(y) coef(lm(y ~ books$year))[[2L]], numeric(1))
round(sd(slopes2), 3)
#> [1] 0.459
round(summary(fit)$coefficients["year", "Std. Error"], 3)
#> [1] 0.459Source: R documentation: simulate
The standard deviation of the parametric bootstrap slopes (0.459) reproduces the standard error that summary() reports (0.459), because the simulated data come from exactly the model whose standard error is being checked. The same results come out on every run because the seed is set.
Pitfalls
The mistakes below are the ones that produce wrong conclusions without producing an error.
-
Independence. The 42 usable loans come from 13 books, and the rows for one book share
pagesandgenre. The standard errors of a book-level predictor computed from loan-level rows are too small, because the effective sample size for that predictor is the number of books. Model one row per book (here, the mean duration per book) or use a mixed model (nlme::lme(),lme4::lmer(); the R FAQ explains whylmer()shows no p-values, R FAQ 7.35: Why are p-values not displayed when using lmer()?). -
No intercept.
y ~ x - 1andy ~ x + 0drop the intercept, and R-squared is then computed against zero: a model that explains nothing can still report 75 per cent (R FAQ 7.41: Why does summary() report strange results for the R^2 estimate when I fit a linear model with no intercept?). -
A numeric code is not a factor. A genre coded 1 to 6 and used as a number is a single slope; wrap it in
factor()to get one coefficient per level. -
In-sample fit. The residual standard error is computed on the rows the model was fitted to and is optimistic about new data. For a linear model the leave-one-out prediction error needs no refits: it is the residual divided by one minus the hat value, which the script checks against actual refits.
-
Different rows in different models.
AIC()of two models fitted to different rows are not comparable; comparenobs(). -
Coefficient tests are not variable tests. A factor with six levels has five coefficients; use
anova()ordrop1()to test the factor as a whole. -
Stepwise selection gives optimistic p-values and intervals, as noted above.
-
Extrapolation. A prediction outside the range of the predictors is a statement about the model, not the data (the intercept of the pages-and-genre model is the duration at zero pages).
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
books$isbn <- normalise_isbn(books$isbn)
loan_books <- merge(loans, books, by = "isbn")
loan_books$days <- as.numeric(loan_books$return_date - loan_books$loan_date)
dur <- loan_books[!is.na(loan_books$days) & loan_books$days >= 0, ]
dur <- dur[complete.cases(dur[, c("days", "pages", "genre")]), ]
# 1. The loans are not independent: a book contributes several rows.
table(table(dur$isbn))
#>
#> 1 2 3 4 6
#> 2 3 4 1 3
length(unique(dur$isbn))
#> [1] 13
# Rows of the same book share pages and genre, so the effective sample size for
# those variables is the number of books. Model one row per book instead.
per_book <- aggregate(days ~ isbn + pages + genre, data = dur, FUN = mean)
nrow(per_book)
#> [1] 13
book_fit <- lm(days ~ pages, data = per_book)
round(coef(summary(book_fit)), 4)
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 18.1409 3.4729 5.2236 0.0003
#> pages 0.0062 0.0123 0.5034 0.6246
loan_fit <- lm(days ~ pages, data = dur)
round(coef(summary(loan_fit)), 4)
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 17.8693 3.7175 4.8068 0.000
#> pages 0.0068 0.0129 0.5258 0.602
# 2. No intercept: R-squared is computed against zero, not against the mean
# (R FAQ 7.41), so it is not comparable with the usual one.
with_int <- lm(days ~ pages, data = dur)
no_int <- lm(days ~ pages - 1, data = dur)
round(c(
with_intercept = summary(with_int)$r.squared,
no_intercept = summary(no_int)$r.squared
), 3)
#> with_intercept no_intercept
#> 0.007 0.751
y <- dur$days
round(1 - sum(residuals(no_int)^2) / sum(y^2), 3)
#> [1] 0.751
round(1 - sum(residuals(no_int)^2) / sum((y - mean(y))^2), 3)
#> [1] -0.567
# 3. A numeric code is not a factor: genre as 1..6 forces a straight line.
dur$genre_code <- as.integer(factor(dur$genre))
length(coef(lm(days ~ genre_code, data = dur)))
#> [1] 2
length(coef(lm(days ~ factor(genre_code), data = dur)))
#> [1] 6
# 4. Fitting and testing on the same rows overstates the fit: compare the
# residual standard error with a leave-one-out prediction error.
loo <- vapply(seq_len(nrow(dur)), \(i) {
m <- lm(days ~ pages, data = dur[-i, ])
dur$days[i] - predict(m, newdata = dur[i, ])
}, numeric(1))
round(c(sigma = sigma(with_int), loo_rmse = sqrt(mean(loo^2))), 2)
#> sigma loo_rmse
#> 8.74 8.96
all.equal(unname(loo), unname(residuals(with_int) / (1 - hatvalues(with_int))))
#> [1] TRUEVersion notes
The modelling functions are old and stable; most NEWS entries since R 4.0.0 are fixes and small additions, listed here from the NEWS files of the corresponding releases (R NEWS). The R 4.0.0 entries are in the same file; entries older than R 4.0.0 are in the R 3.x NEWS.
| Version | Change |
|---|---|
R 4.0.0 |
|
R 4.1.0 |
|
R 4.3.0 |
|
R 4.3.3 |
|
R 4.4.0 |
The |
R 4.5.0 |
New |
R 4.5.2 |
The help pages of |
R 4.6.0 |
The default |
The R-devel NEWS (the upcoming 4.7.0, never run here) lists, among others, a wtol argument for lm() and glm() that treats very small weights as zero; lm(rock) and glm(rock) style calls where the first argument is a data frame, and update() on such a fit; correct AIC() and logLik() for a glm with zero weights; model.frame.glm() using the same logic as model.frame.lm(); a fix for drop1.glm() on models fitted with y = FALSE; anova(<glm>, test = tt) partially matching the string tt; named lower and upper bounds of nls() in any order and for a subset of the parameters; and summary(<aov>) warning on an almost perfect fit, as summary(<lm>) does (R-devel NEWS: Changes in R-devel). The same file’s section for the scheduled R 4.6.2 lists a fix for a long-standing bug in nlminb() with constraints, and an update of influence.measures() to match dffits(), dfbetas() and covratio().
Where older tutorials differ
| Older tutorial writes | Current practice | Why | Since |
|---|---|---|---|
|
|
|
always |
|
|
the |
R 4.4.0 |
|
|
chi-squared for fixed dispersion, F for estimated dispersion are the defaults; |
R 4.4.0 |
|
character columns, turned into factors by the model function (levels sorted alphabetically, stored in |
|
R 4.0.0 |
The residuals-versus-leverage plot has a smoother; the |
no smoother in the leverage plot; a half-normal Q-Q plot of absolute standardised deviance residuals for a |
the NEWS says the smoother never made much sense, and deviance residuals of a |
R 4.6.0, R 4.3.0 |
|
warns only when |
the old warning could not tell estimable from non-estimable predictions |
R 4.3.0 |
|
|
|
always |
Hand-picked starting values for every |
a self-starting model ( |
the start is computed from the data; bounds keep parameters in range |
always |
Type III sums of squares from another package, or the default |
|
R’s default |
always |
|
|
the exponential of the mean of the logarithm is the geometric mean, not the mean |
always |
|
|
a global option changes every later model in the session, including those in other packages |
always |
|
|
|
always |
|
a model chosen from subject knowledge, cross-validation or penalised regression |
the p-values after stepwise selection ignore the search |
always |
Source: R FAQ 3.3.2: Models
Two community habits are worth a sentence each. Tutorials written for other languages often fit a model and then parse the printed summary(); in R, take the number from the summary object (summary(fit)$r.squared) or use broom::glance(). And the "base R, tidyverse and tidymodels" split in modelling tutorials is about workflow (the tidymodels packages wrap fitting engines, including lm() and glm(), behind one interface), not about the fitting functions on this page, which they call.
Where to go next
-
Formulas for the formula language, terms, model frames, model matrices and the contrast matrices behind every factor coefficient.
-
S3 Classes and Generics for the dispatch that makes
summary(fit)findsummary.lm(), and for writing your own methods. -
S4 Classes and Methods for the class system
stats4::mle()uses. -
Descriptive Statistics and Tests for summaries,
t.test(),chisq.test()andaov()withTukeyHSD(). -
Distributions and Random Numbers for the
d,p,qandrfunctions the likelihoods use,set.seed(), andRNGkind(). -
Base Graphics and Devices for drawing the diagnostic plots to a file, and ggplot2 and the Grammar of Graphics for plotting fitted values and intervals.
-
Data Frames and Factors for the data the models consume, and Missing and Special Values for
NAhandling. -
Conditions and Error Handling for
tryCatch()andwithCallingHandlers(), used here to catch the model warnings. -
The tidyverse and dplyr for the data-frame side of tibbles, which broom returns.
References
-
An Introduction to R: 11 Statistical models in R; 11.2 Linear models, 11.3 Generic functions for extracting model information, 11.4.1 ANOVA tables, 11.5 Updating fitted models, 11.6 Generalized linear models, 11.6.1 Families, 11.6.2 The glm() function, 11.7 Nonlinear least squares and maximum likelihood models, 11.7.2 Maximum likelihood, 11.8 Some non-standard models, 11.1.1 Contrasts, Least squares fitting and the QR decomposition
-
R FAQ 3.3.2: Models
-
R FAQ 7.18: Why does the output from anova() depend on the order of factors in the model?
-
R FAQ 7.35: Why are p-values not displayed when using lmer()?
-
R FAQ 7.41: Why does summary() report strange results for the R^2 estimate when I fit a linear model with no intercept?
-
R NEWS: Changes in R 4.0.0 to 4.6.1; R-devel NEWS