Formulas
|
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. Everything on this page is base R (the base, stats, graphics and utils packages that ship with R), and every example runs; the lintr 3.4.0 default linters were used only to check the scripts. Two add-on packages are named in a sentence each and not run: purrr 1.2.2 and ggplot2 4.0.3. A formula is the small piece of R syntax that looks like y ~ x. It is how R programmers say "y depends on x" to a function, and it is the way almost every statistical modelling function in R is told what to fit. It is also an ordinary R object with a class, a length, an environment and a place in the language, and base R reuses it well outside modelling: aggregate(pages ~ genre, ...) summarises by group, xtabs(~ genre, ...) counts, boxplot(pages ~ genre, ...) splits. This page explains what the tilde operator builds, how the formula language expands the operators +, -, :, *, /, ^, %in%, . and I(), how terms(), model.frame() and model.matrix() turn a formula and a data frame into a design matrix, how factors are coded by contrasts, how to change a formula with update() and build one with reformulate(), and how the other formula-using functions of base R read it.
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; a # nolint: comment marks a line where the idiom on show is what a default linter dislikes, so the linter does not hide the lesson. 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"). The books data frame has the columns pages, year and genre that most formulas below use. Its genre column holds nine spellings of six genres ("Fiction", "fiction", "Fiction " and so on), so the scripts that model by genre first run books$genre <- tolower(trimws(books$genre)); cleaning data is the subject of Strings and Regular Expressions and Data Frames. Where a message contains quotation marks, they are the typographic quotes R prints in a UTF-8 locale. A formula created inside a function prints an extra <environment: ...> line holding a memory address that changes on every run; the examples avoid printing that line and use environmentName() or deparse1() instead, and say so where it matters.
This page is the third of the three metaprogramming pages. A formula is a call, so everything in Language Objects and Quoting applies to it, and the way model functions evaluate formulas against a data frame is the data-mask idea of Evaluation and Non-Standard Evaluation. The fitting of models themselves (lm(), glm(), residual diagnostics, anova()) belongs to Statistical Models; this page stops where the model is fitted and only uses lm() to show what a formula does.
What a formula is
The tilde operator
The tilde ~ separates the left-hand side and the right-hand side of a model formula. Its help page is short: the left-hand side is optional, and a formula has mode call whose components are ~, the left-hand side (if present) and the right-hand side, in that order, so a one-sided formula has two components. The operator does not evaluate either side. ~ is a primitive of the kind called a "special" (the same kind as quote() and if): R hands it the unevaluated expressions on both sides, and it returns the call itself, with two additions that make it a formula, a class attribute and an environment. That is why a formula can mention variables that do not exist (no_such_response ~ no_such_predictor is a perfectly good object), why y ~ x costs nothing to create, and why it is the model-fitting function, not the tilde, that later decides what the names mean.
source("scenario.R")
# `~` quotes both sides and returns a formula.
f <- pages ~ genre
f
#> pages ~ genre
class(f)
#> [1] "formula"
typeof(f)
#> [1] "language"
length(f)
#> [1] 3
f[[1]]
#> `~`
f[[2]]
#> pages
f[[3]]
#> genre
# A one-sided formula has no left-hand side: length 2, and [[2]] is the rhs.
g <- ~genre
length(g)
#> [1] 2
g[[2]]
#> genre
# Nothing is evaluated, so the names need not exist.
no_such_response ~ no_such_predictor
#> no_such_response ~ no_such_predictor
exists("no_such_response")
#> [1] FALSE
# `~` is a primitive "special": its arguments are not evaluated.
typeof(`~`)
#> [1] "special"
is.primitive(`~`)
#> [1] TRUESource: R documentation: ~
The tilde has a very low precedence: it binds more loosely than | and || and more tightly than the assignments and ?, so pages ~ year + genre is read as pages ~ (year + genre) without parentheses, and a formula can contain whole logical expressions on either side. The table of operators in Operators and Precedence lists it, and the manual lists it as the operator used for model formulae, either unary or binary (The R Language Definition: Operators). The unary form is the one-sided formula: ~ genre has no response, which is what xtabs() and cor.test() expect, and what model.matrix() needs when there is nothing to predict.
A formula is a call until it is evaluated
Quoting a formula gives a call of class "call", and evaluating that call is what produces the formula. Everything R can do to a call it can do to a formula: index it with [[ (f[[2]] is the left-hand side), replace parts, convert it to a list, ask for all the names in it with all.vars(). This is how the model-fitting functions take a formula apart, and how a program can build or edit one (see Building and inspecting formulas below).
# Quoting `y ~ x` gives a call; evaluating the call makes the formula.
cl <- quote(pages ~ year)
class(cl)
#> [1] "call"
is.call(cl)
#> [1] TRUE
f <- eval(cl)
class(f)
#> [1] "formula"
identical(f, pages ~ year)
#> [1] TRUE
# A formula is a call with two attributes: the class and the environment.
names(attributes(f))
#> [1] "class" ".Environment"
attr(f, "class")
#> [1] "formula"
attr(f, ".Environment")
#> <environment: R_GlobalEnv>
# Because it is a call, it can be indexed and modified like one.
f[[3]] <- quote(year + genre)
f
#> pages ~ year + genre
f[[2]] <- quote(log(pages))
f
#> log(pages) ~ year + genre
as.list(f)
#> [[1]]
#> `~`
#>
#> [[2]]
#> log(pages)
#>
#> [[3]]
#> year + genre
#>
deparse1(f)
#> [1] "log(pages) ~ year + genre"
all.vars(f)
#> [1] "pages" "year" "genre"Source: R documentation: ~
The two attributes are the whole difference between the call quote(pages ~ year) and the formula pages ~ year. print() of a formula (the print.formula() method) shows the environment only when it is not the global environment, through its showEnv argument, whose default is !identical(e, .GlobalEnv); that is the reason a formula built in a function prints an extra line. The same formula can be turned back into text with deparse1(), which is the safest way to show or compare formulas in output.
The environment of a formula
The formula keeps the environment in which it was created, and that is its second purpose. The help page of formula() says it plainly: a formula object has an associated environment, and model.frame() uses that environment, not the parent environment of the call, to evaluate the variables that are not found in the supplied data argument. So when lm(pages ~ year, data = books, subset = year > cutoff) runs, pages and year are looked up in books first, and cutoff, which is not a column, is looked up in the environment of the formula, which is the frame of the function that wrote the formula. The R language definition makes a related point about the call stack: some functions, in particular those that use and manipulate model formulas, need to simulate dynamic scope by accessing the call stack directly (The R Language Definition: The call stack). The environments themselves, and parent.frame(), are the subject of Scoping, Environments and Closures.
source("scenario.R")
# A formula made at the console remembers the global environment.
f <- pages ~ year
environment(f)
#> <environment: R_GlobalEnv>
environmentName(environment(f))
#> [1] "R_GlobalEnv"
# A formula made inside a function remembers that call's frame.
make_formula <- function() {
cutoff <- 2000 # nolint: object_usage_linter.
pages ~ year
}
g <- make_formula()
environmentName(environment(g))
#> [1] ""
identical(environment(g), globalenv())
#> [1] FALSE
ls(environment(g))
#> [1] "cutoff"
get("cutoff", envir = environment(g))
#> [1] 2000
# print() shows the environment when it is not the global one. That line holds
# a memory address, which changes on every run, so the page omits it.
print(g, showEnv = FALSE)
#> pages ~ year
# Variables that are not columns of `data` are looked up in that environment.
fit_since <- function(data) {
cutoff <- 2000
lm(pages ~ year, data = data,
subset = year > cutoff) # nolint: object_usage_linter.
}
coef(fit_since(books))
#> (Intercept) year
#> -5747.76495 2.97732
# The same call fails at the console: `cutoff` does not exist there.
cutoff_missing <- tryCatch(
lm(pages ~ year, data = books, subset = year > cutoff),
error = \(e) conditionMessage(e)
)
cutoff_missing
#> [1] "object 'cutoff' not found"Source: R documentation: formula
Two practical consequences follow. A formula that is stored and used later (in a list, in a package function, saved with saveRDS()) carries a whole environment with it, and everything that environment references stays alive and is serialized with it; a formula built at the console carries only the global environment, which is saved as a reference to the global environment and not as a copy, but a formula built inside a function carries a copy of the function’s frame. And a formula that works in one place can fail in another: the last call above works inside fit_since() and fails at the console, because the two formulas see different environments. When the data and every variable the formula needs are columns of data, none of this matters; it starts to matter as soon as a formula refers to a variable that is not a column, or is passed between functions.
formula() and as.formula()
A string, a fitted model or a data frame can be turned into a formula with formula(), which is a generic function, or as.formula(), which is almost identical and additionally keeps the attributes of an object that already is a formula. The help page describes the methods: for a fitted model the default method looks for a "formula" component, then a "terms" component, then the formula argument of the call, and finally a "formula" attribute; for a data frame the first column becomes the response and the remaining columns the right-hand side, or, with one column, a one-sided formula. DF2formula() is the helper behind the data frame method. The argument env of as.formula() (default parent.frame()) sets the environment of a formula built from a string; an existing formula keeps its own.
source("scenario.R")
# formula() and as.formula() turn a string into a formula.
f <- as.formula("pages ~ year + genre")
f
#> pages ~ year + genre
class(f)
#> [1] "formula"
environmentName(environment(f))
#> [1] "R_GlobalEnv"
# The default `env` is parent.frame(): inside a function, the function's frame.
build <- function(text) as.formula(text)
environmentName(environment(build("pages ~ year")))
#> [1] ""
identical(environment(build("pages ~ year")), globalenv())
#> [1] FALSE
# env makes the environment explicit.
shared <- new.env()
h <- as.formula("pages ~ year", env = shared)
identical(environment(h), shared)
#> [1] TRUE
# A formula passes through unchanged, so as.formula() is safe on either form.
identical(as.formula(f), f)
#> [1] TRUE
identical(formula(f), f)
#> [1] TRUE
# formula() also reads the formula out of a fitted model, and a data frame
# becomes first column ~ the rest.
fit <- lm(pages ~ year, data = books)
formula(fit)
#> pages ~ year
formula(books[c("pages", "year", "genre")])
#> pages ~ year + genre
formula(data.frame(pages = 1))
#> ~pages
# A length-one character vector is the supported input. Longer ones are
# deprecated since R 4.0.0.
res <- tryCatch(
formula(c("pages ~", "year")),
warning = \(w) conditionMessage(w)
)
res
#> [1] "Using formula(x) is deprecated when x is a character vector of length > 1.\n Consider formula(paste(x, collapse = \" \")) instead."Source: R documentation: formula
Two version details are worth knowing (full list in Version notes). Since R 3.6.0 formula(model.frame(frml, ...)) returns the original frml in every case, and the earlier behaviour is available as DF2formula(). And since R 4.0.0 a character vector of length greater than one is deprecated as input: the old behaviour parsed the elements as separate lines of code and used the first complete expression, which silently truncated bad input. The last example above shows the real warning text, which also suggests the replacement, formula(paste(x, collapse = " ")). The help page adds that building the formula by evaluating a call, eval(call("~", quote(foo + bar))), is an order of magnitude more efficient than formula(c("~", "foo + bar")).
The formula language
The right-hand side of a model formula is written in a small language of its own, in which the usual arithmetic operators mean something different. The formula operators are similar in effect to the Wilkinson and Rogers notation of programs such as Glim and Genstat, with one inevitable change: the operator written . there becomes : in R, because the period is a valid character in a name (An Introduction to R: Defining statistical models; formulae). The notation is summarised there (based on Chambers and Hastie, 1992), and the help page ?formula repeats it.
| Written | Meaning | Expands to |
|---|---|---|
|
|
the response is |
|
include |
both sets of terms |
|
include |
|
|
the interaction (tensor product) of |
one term |
|
crossing |
|
|
|
|
|
nesting shorthand |
|
|
all terms in the bracket, with interactions up to order 2 |
|
|
insulate |
one term, whose column is the computed value |
|
on the right of a model fit: all columns of |
see the sections below |
|
remove the intercept |
|
|
a term whose coefficient is fixed at one |
no column in the design matrix |
The best way to learn the language is to look at what terms() makes of an expression, because the term labels are exactly the columns that will be built. The script below applies the operators to the abstract names a, b and c (nothing is evaluated, so they need not exist) and prints the term labels, and then the intercept for the six ways of writing it.
# What each operator expands to: the term labels and the intercept of terms().
expand <- function(f) {
tt <- terms(f)
paste(attr(tt, "term.labels"), collapse = ", ")
}
forms <- list(
y ~ a + b,
y ~ a:b,
y ~ a * b,
y ~ (a + b + c)^2,
y ~ (a + b + c)^2 - a:b,
y ~ a + b %in% a,
y ~ a / b,
y ~ a + b - a,
y ~ a * b - a,
y ~ a + a
)
expansions <- data.frame(
formula = vapply(forms, deparse1, character(1L)),
terms = vapply(forms, expand, character(1L))
)
print(expansions, right = FALSE)
#> formula terms
#> 1 y ~ a + b a, b
#> 2 y ~ a:b a:b
#> 3 y ~ a * b a, b, a:b
#> 4 y ~ (a + b + c)^2 a, b, c, a:b, a:c, b:c
#> 5 y ~ (a + b + c)^2 - a:b a, b, c, a:c, b:c
#> 6 y ~ a + b %in% a a, a:b
#> 7 y ~ a/b a, a:b
#> 8 y ~ a + b - a b
#> 9 y ~ a * b - a b, a:b
#> 10 y ~ a + a a
# The intercept is the attribute "intercept": 1 when present, 0 when removed.
with_without <- list(
y ~ x, y ~ x - 1, y ~ x + 0, y ~ 0 + x, y ~ 1, y ~ x - 1 + 1
)
data.frame(
formula = vapply(with_without, deparse1, character(1L)),
intercept = vapply(
with_without, \(f) attr(terms(f), "intercept"), numeric(1L)
)
)
#> formula intercept
#> 1 y ~ x 1
#> 2 y ~ x - 1 0
#> 3 y ~ x + 0 0
#> 4 y ~ 0 + x 0
#> 5 y ~ 1 1
#> 6 y ~ x - 1 + 1 1Source: R documentation: terms.formula
A few lines of the table deserve a remark. Row 8, y ~ a + b - a, shows that - removes a term, not a variable, and row 9, y ~ a * b - a, shows the trap in the language: removing the main effect a while keeping the interaction a:b produces a model that violates marginality. R does what it is told, and An Introduction to R warns that fitting a model with an interaction but without the corresponding main effects will in general lead to surprising results and is for experts only. Row 10 shows that a repeated term is counted once. The intercept rows agree with the R FAQ on the differences from S: in R y ~ x + 0 is an alternative to y ~ x - 1, and a model with no parameters at all is written y ~ 0 (R FAQ 3.3.2: Models). The last row, y ~ x - 1 + 1, shows that later terms win: the 1 puts the intercept back.
Arithmetic inside a formula: I() and function calls
The symbols +, -, *, ^ and / therefore have two lives. In the formula language they combine terms; inside the parentheses of a function call they are plain arithmetic, so log(pages), sqrt(year) or poly(year, 2) are calls whose value becomes a column, and the label of that term is the text of the call. To use an operator arithmetically at the top level of a term the expression is wrapped in I(), the "as is" function (R documentation: I); it adds the class "AsIs" to its argument, and model code treats the result as the plain value. The R FAQ states the classic difference from S: where S regressed on a power with lm(y ~ x^3), R needs the power insulated, lm(y ~ I(x^3)).
source("scenario.R")
# Inside a formula, ^ means "crossing to degree 2", not squaring.
attr(terms(pages ~ year^2), "term.labels")
#> [1] "year"
# I() ("as is") protects arithmetic.
attr(terms(pages ~ year + I(year^2)), "term.labels")
#> [1] "year" "I(year^2)"
# The columns of the design matrix show the difference.
colnames(model.matrix(pages ~ year^2, data = books))
#> [1] "(Intercept)" "year"
colnames(model.matrix(pages ~ year + I(year^2), data = books))
#> [1] "(Intercept)" "year" "I(year^2)"
colnames(model.matrix(pages ~ I(year - 2000), data = books))
#> [1] "(Intercept)" "I(year - 2000)"
# Functions of a variable are terms of their own. poly() does not accept NA,
# and it runs before any na.action can drop the row (year is NA in one book).
poly_error <- tryCatch(
model.matrix(pages ~ poly(year, 2), data = books),
error = \(e) conditionMessage(e)
)
poly_error
#> [1] "missing values are not allowed in 'poly'"
dated <- books[!is.na(books$year), ]
colnames(model.matrix(pages ~ log(year) + poly(year, 2), data = dated))
#> [1] "(Intercept)" "log(year)" "poly(year, 2)1" "poly(year, 2)2"
# a + b means two terms; I(a + b) means one term, the sum.
d <- data.frame(y = 1:4, a = c(1, 2, 3, 4), b = c(2, 1, 4, 3))
colnames(model.matrix(y ~ a + b, d))
#> [1] "(Intercept)" "a" "b"
colnames(model.matrix(y ~ I(a + b), d))
#> [1] "(Intercept)" "I(a + b)"
model.matrix(y ~ I(a + b), d)[, 2]
#> 1 2 3 4
#> 3 3 7 7
# AsIs objects keep the class "AsIs" and print like the value.
class(I(books$year[1:3]))
#> [1] "AsIs"
I(books$year[1:3])
#> [1] 1974 1813 2019
# Variable names that are not syntactic are quoted with backticks.
odd <- data.frame(y = 1:3, `page count` = c(10, 20, 35), check.names = FALSE)
coef(lm(y ~ `page count`, odd))
#> (Intercept) `page count`
#> 0.28947368 0.07894737Source: R documentation: I
Note the first lines: pages ~ year^2 does not square the year, because ^ means "crossing to the given degree", and crossing a single variable with itself gives the variable, so the model has the single column year. The poly() error is a real subtlety of how formulas are evaluated: the variables in the formula (including calls such as poly(year, 2)) are evaluated first, in model.frame(), and rows with missing values are dropped afterwards by the na.action, so a function that rejects NA fails before the action can drop the row. The help page of model.frame() adds a related warning: because the variables are evaluated before the rows are dropped by subset, data-dependent bases such as orthogonal polynomials (poly()) or splines (bs() of the splines package) are computed on the full data set rather than on the subset, and the same holds for summary statistics such as mean(). A variable name that is not syntactic is written between backticks, as in the last example; the help page of formula() cautions that not all code that uses formulas will accept such names.
The dot, the intercept and the other conventions
The symbol . has two special meanings, and the help page of formula() keeps them apart. In a formula given to a model-fitting function together with a data argument it stands for all the columns of data that are not otherwise in the formula, and terms(formula, data = data) expands it. In the context of update.formula() only, it stands for whatever was previously in that part of the formula (see Updating formulas and fitted models). Without data, terms() cannot expand it, and says so.
The intercept is always part of the model unless it is removed with -1, + 0 or 0 +; a model without an intercept changes how summary.lm() computes R squared, which is the subject of R FAQ 7.41 and is shown in Pitfalls.
Anything else in a formula is a call to an ordinary function: the vertical bar in y ~ a | b has no meaning to terms() (lattice, which ships with R, gives it one), and s(a) is a term labelled s(a) unless it is declared a "special" in terms(). The specials argument marks calls of the named functions in the "specials" attribute, which is how aov() finds Error() terms and how the survival package (3.8-6, which ships with R) finds strata() and cluster(). The help page cautions that specials are nothing more than syntax, so y ~ stats::offset(w) + x does not fix the coefficient of w: the call is to stats::offset, not to offset.
# `|` is just a function call to terms(); packages give it a meaning.
attr(terms(y ~ a | b), "term.labels")
#> [1] "a | b"
# So are calls to any function, unless they are declared as specials.
attr(terms(y ~ s(a) + b), "term.labels")
#> [1] "s(a)" "b"
attr(terms(y ~ s(a) + b, specials = "s"), "specials")
#> $s
#> [1] 2
#>
# Backtick-quoted names that are not syntactic.
attr(terms(y ~ `page count` + `year`), "term.labels")
#> [1] "`page count`" "year"
# A formula with no variables at all is allowed.
attr(terms(~1), "term.labels")
#> character(0)
attr(terms(y ~ 1), "intercept")
#> [1] 1
length(attr(terms(~0), "term.labels"))
#> [1] 0Source: R documentation: terms.formula
Terms objects
terms() is a generic function whose formula method does the symbolic work: it expands the operators, orders the terms, replaces . and records what the model needs. Its value is a terms object, which is the formula itself with attributes; the help page notes that, unlike in S, there is no formula attribute: the object is itself a formula. Every model fit stores one (fit$terms), and model.frame() and model.matrix() read it.
source("scenario.R")
tt <- terms(pages ~ genre * year)
class(tt)
#> [1] "terms" "formula"
tt[[2]]
#> pages
names(attributes(tt))
#> [1] "variables" "factors" "term.labels" "order" "intercept"
#> [6] "response" "class" ".Environment"
# The pieces you are most likely to use.
attr(tt, "term.labels")
#> [1] "genre" "year" "genre:year"
attr(tt, "order")
#> [1] 1 1 2
attr(tt, "intercept")
#> [1] 1
attr(tt, "response")
#> [1] 1
attr(tt, "variables")
#> list(pages, genre, year)
attr(tt, "factors")
#> genre year genre:year
#> pages 0 0 0
#> genre 1 0 1
#> year 0 1 1
# The environment attribute is the formula's environment.
environmentName(attr(tt, ".Environment"))
#> [1] "R_GlobalEnv"
# formula() turns a terms object back into a plain formula, and simplify = TRUE
# writes the expansion out.
formula(tt)
#> pages ~ genre * year
terms(pages ~ genre * year, simplify = TRUE)[[3]]
#> genre + year + genre:yearSource: R documentation: terms.object
The attributes, in the order the help page of terms.object lists them:
| Attribute | Content |
|---|---|
|
a call to |
|
an integer matrix of variables by terms: |
|
the labels of the terms, after re-ordering, without offsets; non-syntactic names are quoted with backticks so that the formula can be rebuilt from them |
|
the order of interaction of each term ( |
|
|
|
the index in |
|
present when the model has |
|
present when |
|
optional; set by model fitting to help |
|
the environment of the formula |
Source: R documentation: terms.object
The argument data is used to expand the dot, keep.order keeps the terms in the order written instead of main effects first, simplify = TRUE returns the formula expanded and simplified, specials marks special calls, and allowDotAsName = TRUE treats . as an ordinary name for non-standard uses of formulas. The arguments abb and neg.out are not implemented in R and have been formally deprecated since R 4.4.0.
source("scenario.R")
# `.` means "every column of data that is not already used".
numbers <- books[c("pages", "year", "genre")]
tt <- terms(pages ~ ., data = numbers)
formula(tt)
#> pages ~ year + genre
attr(tt, "term.labels")
#> [1] "year" "genre"
# With the left-hand side removed from the pool, only the other columns count.
formula(terms(log(pages) ~ ., data = numbers))
#> log(pages) ~ year + genre
formula(terms(pages ~ . - genre, data = numbers))
#> pages ~ (year + genre) - genre
formula(terms(pages ~ .^2, data = numbers))
#> pages ~ (year + genre)^2
# Without data, `.` cannot be expanded.
no_data <- tryCatch(
terms(pages ~ .),
error = \(e) conditionMessage(e)
)
no_data
#> [1] "'.' in formula and no 'data' argument"
# Term order: main effects first, then interactions, unless keep.order = TRUE.
attr(terms(y ~ a:b + a + b), "term.labels")
#> [1] "a" "b" "a:b"
attr(terms(y ~ a:b + a + b, keep.order = TRUE), "term.labels")
#> [1] "a:b" "a" "b"
# Marking special function calls.
tt <- terms(y ~ a + offset(w) + strata(g), specials = c("strata", "offset"))
attr(tt, "specials")
#> $strata
#> [1] 4
#>
#> $offset
#> [1] 3
#>
attr(tt, "offset")
#> [1] 3Source: R documentation: terms.formula
The line pages ~ (year + genre) - genre is how the dot expansion looks: the columns of data that are not the response are put between parentheses, and the operators then act on the group, which is also why pages ~ .^2 gives all two-way interactions of the other columns. Marking a function as special is only syntax, and R 4.6.0 made terms(<formula>, specials = ) treat non-syntactic specials more gracefully; the help page of terms.formula adds that since R 4.6.0 the entry "<pkg>::<name>" in specials marks calls of the backquoted function name, not calls written <pkg>::<name>(...); the reverse was true before R 4.6.0, but that was undocumented and accidental.
The factors matrix shows how R decides between contrast coding and dummy coding for a factor, which matters in the design matrix below. The entry 2 appears when a variable occurs in an interaction whose margin is missing from the model:
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
three <- books[books$genre %in% c("classic", "fiction", "poetry"), ]
three$era <- factor(ifelse(three$year >= 2000, "new", "old"))
# In "factors", 1 means "coded by contrasts" and 2 means "coded by dummy
# variables for every level", which happens when the margin is missing.
attr(terms(~ genre + genre:era), "factors")
#> genre genre:era
#> genre 1 2
#> era 0 1
attr(terms(~ genre:era), "factors")
#> genre:era
#> genre 2
#> era 2
# The same difference in the design matrix.
colnames(model.matrix(~ genre + genre:era, data = three))
#> [1] "(Intercept)" "genrefiction" "genrepoetry"
#> [4] "genreclassic:eraold" "genrefiction:eraold" "genrepoetry:eraold"
colnames(model.matrix(~ genre:era, data = three))
#> [1] "(Intercept)" "genreclassic:eranew" "genrefiction:eranew"
#> [4] "genrepoetry:eranew" "genreclassic:eraold" "genrefiction:eraold"
#> [7] "genrepoetry:eraold"Source: R documentation: terms.object
The rule is about margins. In genre:era the factor era is coded by contrasts (1) when the term left after removing era, which is genre, is in the model, and the factor genre is coded by dummy variables (2) when the term left after removing genre, which is era, is not. In the first model the interaction therefore has one column for each level of genre, times the single contrast column of era; in the second model neither margin is present, so both factors are coded by dummy variables and the interaction has one column for each of the six combinations. R tracks this automatically, which is why terms(), and not the formula text, is the reliable description of the model.
Changing a terms object
delete.response() returns a terms object with no response, which is what predict() uses on the right-hand side to build the design matrix of new data. drop.terms() removes terms by position in labels(termobj) and optionally keeps the response, and the [ method of terms selects terms and always keeps the response. R 4.3.3 fixed several bugs in drop.terms() (dropx of length zero, offsets, no invented response), and R 4.4.0 made it honour keep.response = FALSE, its documented default; R 4.5.0 made drop.terms(y ~ w, 1) work.
tt <- terms(y ~ z + x + w)
tt
#> y ~ z + x + w
#> attr(,"variables")
#> list(y, z, x, w)
#> attr(,"factors")
#> z x w
#> y 0 0 0
#> z 1 0 0
#> x 0 1 0
#> w 0 0 1
#> attr(,"term.labels")
#> [1] "z" "x" "w"
#> attr(,"order")
#> [1] 1 1 1
#> attr(,"intercept")
#> [1] 1
#> attr(,"response")
#> [1] 1
#> attr(,".Environment")
#> <environment: R_GlobalEnv>
# delete.response() keeps the right-hand side; it is what predict() uses.
nr <- delete.response(tt)
attr(nr, "response")
#> [1] 0
formula(nr)
#> ~z + x + w
# drop.terms() removes terms by position. keep.response = FALSE is the
# documented default, and R 4.4.0 made the function honour it.
formula(drop.terms(tt, 2:3, keep.response = TRUE))
#> y ~ z
formula(drop.terms(tt, 2:3))
#> ~z
# [.terms keeps the response, always.
formula(tt[-1])
#> y ~ x + w
formula(tt[2:3])
#> y ~ x + wFrom a formula to a model frame
model.frame() evaluates the variables of the formula and returns a data frame with one column per variable or expression, in the order of the variables attribute, plus the terms object as the "terms" attribute. The procedure is given in the help page: all the variables in formula, subset and the extra arguments are looked up first in data and then in the environment of the formula; then the subset expression is evaluated and used as a row index; then the na.action function is applied (and may add attributes to the frame); and finally the levels of factors are adjusted according to drop.unused.levels and xlev. The na.action is taken, in order, from an attribute of data, from options("na.action"), and then na.fail; the factory-fresh option is na.omit, so rows with a missing value in any variable of the model disappear silently, and the number of observations can be smaller than the number of rows. The "na.action" attribute of the result records which rows went (this is what summary.lm() reports as observations deleted due to missingness).
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
mf <- model.frame(pages ~ year + genre, data = books)
class(mf)
#> [1] "data.frame"
dim(mf)
#> [1] 13 3
names(mf)
#> [1] "pages" "year" "genre"
head(mf, 4)
#> pages year genre
#> 1 387 1974 science fiction
#> 2 432 1813 classic
#> 3 214 2019 fiction
#> 4 198 2021 fiction
# 16 books, 13 rows: the NA in pages (two books) and in year (one book)
# are gone.
nrow(books)
#> [1] 16
na.action(mf)
#> 6 7 15
#> 6 7 15
#> attr(,"class")
#> [1] "omit"
class(na.action(mf))
#> [1] "omit"
rownames(books)[unclass(na.action(mf))]
#> [1] "6" "7" "15"
# na.action = na.pass keeps them.
dim(model.frame(pages ~ year + genre, data = books, na.action = na.pass))
#> [1] 16 3
getOption("na.action")
#> [1] "na.omit"
# subset is evaluated in the data (then in the formula's environment).
dim(model.frame(pages ~ year, data = books, subset = year > 2015))
#> [1] 7 2
# Transformations are computed once, and the column is named after the call.
names(model.frame(log(pages) ~ I(year - 2000) + genre, data = books))
#> [1] "log(pages)" "I(year - 2000)" "genre"
# The response and other extras have accessors.
y <- model.response(mf)
head(y)
#> 1 2 3 4 5 8
#> 387 432 214 198 301 280
mfw <- model.frame(pages ~ year, books, weights = pages / 100)
names(mfw)
#> [1] "pages" "year" "(weights)"
w <- model.weights(mfw)
head(w, 3)
#> [1] 3.87 4.32 2.14
# The terms object travels with the frame, so model.matrix() can reuse it.
class(attr(mf, "terms"))
#> [1] "terms" "formula"
formula(attr(mf, "terms"))
#> pages ~ year + genre
# get_all_vars() only collects the columns: no transformation, no NA drop.
dim(get_all_vars(log(pages) ~ year, data = books))
#> [1] 16 2
# A name that exists nowhere is an error with a clear message.
missing_var <- tryCatch(
model.frame(pages ~ year + chapters, data = books),
error = \(e) conditionMessage(e)
)
missing_var
#> [1] "object 'chapters' not found"Source: R documentation: model.frame
Extra named arguments passed to model.frame() (the weights and offset of a model, for instance) become further columns with parenthesised names, such as "(weights)" above; model.response(), model.weights() and model.offset() extract them (R documentation: model.extract). get_all_vars() returns only the columns the formula refers to, without calling the functions in it, which is useful when the raw variables are wanted. Only variables of type raw, logical, integer, real, complex or character can be in a model frame: that includes factors, and excludes lists. Note also that every variable named in the formula is in the frame, even if it is preceded by -. R 4.5.0 made model.frame() give more informative errors in some cases when variables in the formula are not found; for a column that does not exist at all, as above, the message is the plain "object 'chapters' not found", so check the spelling of the columns before suspecting anything else.
The model frame is what modelling functions keep with model = TRUE and what expand.model.frame() extends; nrow(model.frame(fit)) is the number of observations that were used in the fit, the number to report.
The model matrix
model.matrix() converts the terms and a model frame into the numeric design matrix, the matrix X of the equation in An Introduction to R, in which each term defines a collection of columns. Character variables are coerced to factors; after that, all variables on the right-hand side must be logical, integer, numeric or factor. A factor of k levels contributes k - 1 columns coded by the contrasts in force (unordered factors by treatment contrasts by default, as the next section explains), unless the intercept is missing, in which case the first factor is coded by dummy variables for all its levels. An interaction contributes the products of the columns of its components. The matrix has two attributes: "assign", an integer vector giving for each column the term that created it (0 is the intercept and positive values count the terms in the order of term.labels), and "contrasts", a named list that says which contrasts were used for each factor.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
three <- books[books$genre %in% c("classic", "fiction", "poetry"), ]
# A character column becomes a factor: one dummy column per level but the first.
mm <- model.matrix(~genre, data = three)
mm
#> (Intercept) genrefiction genrepoetry
#> 2 1 0 0
#> 3 1 1 0
#> 4 1 1 0
#> 7 1 0 1
#> 8 1 0 0
#> 10 1 1 0
#> 11 1 0 1
#> 15 1 1 0
#> attr(,"assign")
#> [1] 0 1 1
#> attr(,"contrasts")
#> attr(,"contrasts")$genre
#> [1] "contr.treatment"
#>
attr(mm, "assign")
#> [1] 0 1 1
attr(mm, "contrasts")
#> $genre
#> [1] "contr.treatment"
#>
# Without the intercept every level gets a column (full dummy coding).
model.matrix(~ 0 + genre, data = three)[1:4, ]
#> genreclassic genrefiction genrepoetry
#> 2 1 0 0
#> 3 0 1 0
#> 4 0 1 0
#> 7 0 0 1
# A numeric variable stays one column; an interaction multiplies columns.
mm2 <- model.matrix(pages ~ genre * year, data = three)
colnames(mm2)
#> [1] "(Intercept)" "genrefiction" "genrepoetry"
#> [4] "year" "genrefiction:year" "genrepoetry:year"
attr(mm2, "assign")
#> [1] 0 1 1 2 3 3
head(mm2, 3)
#> (Intercept) genrefiction genrepoetry year genrefiction:year genrepoetry:year
#> 2 1 0 0 1813 0 0
#> 3 1 1 0 2019 2019 0
#> 4 1 1 0 2021 2021 0
# Missing values are dropped by model.frame() first, so rows can disappear.
dim(model.matrix(pages ~ year, data = books))
#> [1] 13 2
dim(model.matrix(~year, data = books))
#> [1] 15 2
dim(model.matrix(~genre, data = books))
#> [1] 16 6Source: R documentation: model.matrix
The help page states a rule that explains the column order of interactions: the variable whose levels vary fastest is the first one to appear in the formula, not in the term, so in ~ a + b + b:a the interaction has a varying fastest. It also says that, by convention, a response variable that also appears on the right-hand side is dropped from there with a warning, although interactions involving it are retained. The number of rows of the last three calls is the practical thing to remember: pages ~ year loses the three rows with a missing pages or year, the one-sided ~year loses only the row with a missing year, and ~genre loses none because every genre is present. When lm() reports fewer observations than the data frame has rows, this is the reason. For a very large model with mostly zero columns, the Matrix package provides sparse.model.matrix(); the help page points to it.
Contrasts
How is a factor with k levels turned into columns? An Introduction to R explains it (An Introduction to R: Contrasts): an unordered factor generates k - 1 columns, the indicators of its second to last levels, so that the implicit parameterisation contrasts the response at each level with the response at the first; an ordered factor generates k - 1 columns that are the orthogonal polynomials on 1, ..., k without the constant term. The setting is options("contrasts"), a pair of function names for unordered and ordered factors, with the default c("contr.treatment", "contr.poly"). The manual adds a difference from the S language that older textbooks may follow: S used Helmert contrasts for unordered factors, and R’s default is a deliberate choice because treatment contrasts are thought easier for newcomers to interpret.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
genre <- factor(books$genre)
# The defaults: treatment contrasts for unordered factors, polynomial for
# ordered ones.
options("contrasts")
#> $contrasts
#> unordered ordered
#> "contr.treatment" "contr.poly"
#>
# The contrast matrix of a factor: one row per level, one column per
# coefficient.
levels(genre)
#> [1] "classic" "fiction" "mystery" "non-fiction"
#> [5] "poetry" "science fiction"
contrasts(genre)
#> 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
#> non-fiction 0 0 1 0 0
#> poetry 0 0 0 1 0
#> science fiction 0 0 0 0 1
# The built-in coding functions, for three levels.
contr.treatment(3)
#> 2 3
#> 1 0 0
#> 2 1 0
#> 3 0 1
contr.SAS(3)
#> 1 2
#> 1 1 0
#> 2 0 1
#> 3 0 0
contr.sum(3)
#> [,1] [,2]
#> 1 1 0
#> 2 0 1
#> 3 -1 -1
contr.helmert(3)
#> [,1] [,2]
#> 1 -1 -1
#> 2 1 -1
#> 3 0 2
round(contr.poly(3), 3)
#> .L .Q
#> [1,] -0.707 0.408
#> [2,] 0.000 -0.816
#> [3,] 0.707 0.408
contr.treatment(c("low", "mid", "high"), base = 2)
#> low high
#> low 1 0
#> mid 0 0
#> high 0 1Source: R documentation: contrast
The coding functions are all documented in the same help page. Each takes the number of levels (or the vector of level names) and returns a matrix with one row per level and one column per coefficient: contr.treatment() compares each level with a base level (the first by default; the base argument chooses another), contr.SAS() is the same with the last level as base, contr.sum() compares each level with the grand mean and makes the last level the negative sum of the others, contr.helmert() compares each level with the mean of the previous ones, and contr.poly() gives orthogonal polynomials, labelled .L, .Q, .C, ^4 and so on for linear, quadratic, cubic and higher degrees.
The coding changes what the coefficients mean, but not the model: the fitted values, residuals and tests of the whole factor are the same. The script below fits one model with three codings and compares the coefficients against the group means.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
three <- books[books$genre %in% c("classic", "fiction", "poetry"), ]
three$genre <- factor(three$genre)
# Same model, three codings of the factor: different coefficients, same fit.
fit_treatment <- lm(pages ~ genre, data = three)
fit_sum <- lm(pages ~ genre, data = three,
contrasts = list(genre = "contr.sum"))
fit_helmert <- lm(pages ~ genre, data = three,
contrasts = list(genre = "contr.helmert"))
round(rbind(
treatment = coef(fit_treatment),
sum = coef(fit_sum),
helmert = coef(fit_helmert)
), 2)
#> (Intercept) genrefiction genrepoetry
#> treatment 356.00 -133.33 -264.00
#> sum 223.56 132.44 -0.89
#> helmert 223.56 -66.67 -65.78
all.equal(fitted(fit_treatment), fitted(fit_sum))
#> [1] TRUE
all.equal(fitted(fit_treatment), fitted(fit_helmert))
#> [1] TRUE
# Group means, to read the coefficients against.
tapply(three$pages, three$genre, mean, na.rm = TRUE)
#> classic fiction poetry
#> 356.0000 222.6667 92.0000
# Setting the coding on the factor itself.
contrasts(three$genre) <- contr.sum(3)
attr(three$genre, "contrasts")
#> [,1] [,2]
#> classic 1 0
#> fiction 0 1
#> poetry -1 -1
coef(lm(pages ~ genre, data = three))
#> (Intercept) genre1 genre2
#> 223.5555556 132.4444444 -0.8888889
# C() does it in the formula and leaves the data alone.
three$genre <- factor(three$genre)
coef(lm(pages ~ C(genre, sum), data = three))
#> (Intercept) C(genre, sum)1 C(genre, sum)2
#> 223.5555556 132.4444444 -0.8888889
coef(lm(pages ~ C(genre, base = 3), data = three))
#> (Intercept) C(genre, base = 3)classic C(genre, base = 3)fiction
#> 92.0000 264.0000 130.6667
# model.matrix(contrasts.arg =) overrides what is stored on the factor.
model.matrix(~genre, three, contrasts.arg = list(genre = "contr.sum"))[1:3, ]
#> (Intercept) genre1 genre2
#> 2 1 1 0
#> 3 1 0 1
#> 4 1 0 1Source: R documentation: contrasts
Read the coefficients against the group means (classic 356.0, fiction 222.7, poetry 92.0): the treatment coding gives the first mean as the intercept and the differences from it (222.7 - 356.0 = -133.3); the sum coding gives the average of the three means (223.6) as the intercept and each level’s difference from it, with the last level implied; the Helmert coding gives differences from the running mean of the earlier levels. A coefficient in a model with an interaction or in a model without a baseline cannot be interpreted without knowing the coding, which is the point of having "contrasts" attached to the model matrix. The model matrix of the previous section is built with the contrasts of each factor, and there are three ways to change them, from the widest to the narrowest:
-
options(contrasts = c(unordered, ordered))changes the default for every later model fit in the session. -
contrasts(f) <- valuestores a contrast matrix (or the name of a coding function) in the factorfitself; the"contrasts"attribute travels with the factor into the model frame. Ahow.manyargument limits the number of columns. -
C(f, contr, how.many)sets the contrasts offinside the formula, leaving the data untouched (C(genre, sum)orC(genre, base = 3)), andcontrasts.arg = list(f = ...)passes them tomodel.matrix()and to the fitting functions (lm(..., contrasts = list(...))).
The argument contrasts.arg of model.matrix() overrides both the default and any contrasts stored on the factor. R 4.3.0 made C() and contrasts<- preserve the level names when given a function object as well as a function name (before that only a name did), and R 3.6.0 started to warn about an invalid contrasts.arg.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
three <- books[books$genre %in% c("classic", "fiction", "poetry"), ]
# options(contrasts =) changes the default for every later model.
old <- options(contrasts = c("contr.sum", "contr.poly"))
getOption("contrasts")
#> [1] "contr.sum" "contr.poly"
coef(lm(pages ~ genre, data = three))
#> (Intercept) genre1 genre2
#> 223.5555556 132.4444444 -0.8888889
options(old)
getOption("contrasts")
#> unordered ordered
#> "contr.treatment" "contr.poly"
coef(lm(pages ~ genre, data = three))
#> (Intercept) genrefiction genrepoetry
#> 356.0000 -133.3333 -264.0000
# An ordered factor gets polynomial contrasts by default.
size <- factor(
c("short", "medium", "long", "medium"),
levels = c("short", "medium", "long"),
ordered = TRUE
)
colnames(contrasts(size))
#> [1] ".L" ".Q"
round(contrasts(size), 3)
#> .L .Q
#> [1,] -0.707 0.408
#> [2,] 0.000 -0.816
#> [3,] 0.707 0.408
# An invalid contrasts.arg that is not a list is ignored with a warning; a list
# naming a function that does not exist is an error.
tryCatch(
model.matrix(~genre, three, contrasts.arg = "contr.nothing"),
warning = \(w) conditionMessage(w)
)
#> [1] "non-list contrasts argument ignored"
tryCatch(
model.matrix(~genre, three, contrasts.arg = list(genre = "contr.nothing")),
error = \(e) conditionMessage(e)
)
#> [1] "object 'contr.nothing' of mode 'function' was not found"Source: R documentation: C
The two error cases show the difference: an argument that is not a list is ignored with a warning (the model is built with the default contrasts), while a list entry that names a function that does not exist is an error from get(). The ordered factor size uses contr.poly, which is the default for ordered factors; options(contrasts = ...) changes the first element for unordered factors and the second for ordered ones. dummy.coef() (R documentation: dummy.coef) translates the coefficients back to a coefficient per level for a fitted model, whatever the coding; since R 4.4.0 it also works for models with character predictors.
How modelling functions use a formula
The functions of the lm() family do not interpret the formula themselves. They all follow the same recipe, which is worth knowing because a package that offers a new kind of model almost always copies it: capture the call, keep the arguments that model.frame() understands, call model.frame() with them, extract the response with model.response(), build the design matrix with model.matrix(), and hand the numbers to the fitting code. The script below is a complete model function that fits a least-squares line that way; its first lines are the standard idiom (match.call(), match() on the argument names, mf[[1L]] <- quote(stats::model.frame) and eval(mf, parent.frame())), which also explains why subset and na.action are arguments of lm() and why subset can use the names of the columns. The machinery of match.call(), eval() and parent.frame() is the subject of Evaluation and Non-Standard Evaluation.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
# The skeleton of nearly every modelling function: capture the call, keep the
# formula-related arguments, turn the call into model.frame(), evaluate it in
# the caller's frame, then pull out the response and the design matrix.
slope_fit <- function(formula, data, subset,
na.action) { # nolint: object_name_linter.
cl <- match.call()
mf <- match.call(expand.dots = FALSE)
keep <- match(c("formula", "data", "subset", "na.action"), names(mf), 0L)
mf <- mf[c(1L, keep)]
mf[[1L]] <- quote(stats::model.frame)
mf <- eval(mf, parent.frame())
y <- model.response(mf)
x <- model.matrix(attr(mf, "terms"), mf)
beta <- solve(crossprod(x), crossprod(x, y))[, 1L]
list(call = cl, coefficients = beta, n = length(y))
}
res <- slope_fit(pages ~ year, data = books)
res$call
#> slope_fit(formula = pages ~ year, data = books)
res$n
#> [1] 13
res$coefficients
#> (Intercept) year
#> 2067.0926829 -0.8945014
all.equal(res$coefficients, coef(lm(pages ~ year, data = books)))
#> [1] TRUE
# subset reaches model.frame() as an unevaluated expression, so it is looked
# up among the columns of `data`.
slope_fit(pages ~ year, data = books, subset = year > 2000)$n
#> [1] 10Source: R documentation: model.frame
lm() is that recipe plus least squares, and what it stores is what the formula produced. The object keeps terms, model.frame() can rebuild the frame, model.matrix() the design matrix, na.action records the dropped rows, and predict() uses the stored terms to build the design matrix of new data.
source("scenario.R")
fit <- lm(pages ~ year, data = books)
fit
#>
#> Call:
#> lm(formula = pages ~ year, data = books)
#>
#> Coefficients:
#> (Intercept) year
#> 2067.0927 -0.8945
#>
coef(fit)
#> (Intercept) year
#> 2067.0926829 -0.8945014
round(confint(fit), 3)
#> 2.5 % 97.5 %
#> (Intercept) 51.275 4082.910
#> year -1.904 0.115
nobs(fit)
#> [1] 13
# What lm() stored from the formula.
class(fit$terms)
#> [1] "terms" "formula"
formula(fit)
#> pages ~ year
dim(model.frame(fit))
#> [1] 13 2
fit$na.action
#> 6 7 15
#> 6 7 15
#> attr(,"class")
#> [1] "omit"
head(model.matrix(fit), 3)
#> (Intercept) year
#> 1 1 1974
#> 2 1 1813
#> 3 1 2019
round(head(residuals(fit), 3), 2)
#> 1 2 3
#> 85.65 -13.36 -47.09
# predict() reuses the terms object with the response removed.
predict(fit, newdata = data.frame(year = c(1990, 2020)))
#> 1 2
#> 287.0350 260.1999
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.07707An Introduction to R lists the generic functions that extract information from a fitted "lm" object (An Introduction to R: Generic functions for extracting model information): coef(), formula(), residuals(), predict(object, newdata = data.frame), summary(), anova() and the others. predict() is the one that depends most on the formula, because it must reproduce the same columns from new data. For that the fitted object stores the levels of every factor (xlevels) and the terms object, and predict.lm() applies delete.response() to the terms before building the matrix. The new data must hold every variable of the right-hand side with the same names, and a factor level that was not seen in fitting is an error:
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
fit <- lm(pages ~ year + genre, data = books)
# The factor levels seen at fitting time are stored with the model.
fit$xlevels
#> $genre
#> [1] "classic" "fiction" "mystery" "non-fiction"
#> [5] "poetry" "science fiction"
#>
# predict() rebuilds the design matrix from the stored terms and levels.
predict(fit, newdata = data.frame(year = 2020, genre = "poetry"))
#> 1
#> 75.61725
# A level that was not in the training data is an error.
unseen <- tryCatch(
predict(fit, newdata = data.frame(year = 2020, genre = "cookery")),
error = \(e) conditionMessage(e)
)
unseen
#> [1] "factor genre has new level cookery"
# A missing column in newdata is an error too.
no_year <- tryCatch(
predict(fit, newdata = data.frame(genre = "poetry")),
error = \(e) conditionMessage(e)
)
no_year
#> [1] "object 'year' not found"Source: R documentation: predict.lm
offset()
A term wrapped in offset() has its coefficient fixed at one. Modelling functions such as glm() use it for exposure (for instance offset(log(population)) in a Poisson model). In terms() the offset has an entry in the "offset" attribute but not in term.labels; model.frame() keeps its column under the label of the call; model.matrix() has no column for it; and model.offset() extracts the sum of the offsets. The offset argument of the fitting function is an equivalent way, and the two add up.
source("scenario.R")
# offset() is a term whose coefficient is fixed at 1.
tt <- terms(pages ~ year + offset(0.5 * year))
attr(tt, "term.labels")
#> [1] "year"
attr(tt, "offset")
#> [1] 3
mf <- model.frame(pages ~ year + offset(0.5 * year), data = books)
names(mf)
#> [1] "pages" "year" "offset(0.5 * year)"
head(model.offset(mf), 3)
#> [1] 987.0 906.5 1009.5
# The design matrix has no column for the offset; the fit subtracts it.
colnames(model.matrix(pages ~ year + offset(0.5 * year), data = books))
#> [1] "(Intercept)" "year"
fit_plain <- lm(pages ~ year, data = books)
fit_offset <- lm(pages ~ year + offset(0.5 * year), data = books)
coef(fit_plain)
#> (Intercept) year
#> 2067.0926829 -0.8945014
coef(fit_offset)
#> (Intercept) year
#> 2067.092683 -1.394501
# Same fitted values, because the slope moved by exactly the offset's 0.5.
all.equal(fitted(fit_plain), fitted(fit_offset))
#> [1] TRUESource: R documentation: offset
Because the offset is a fixed part of the linear predictor, the model with an offset of 0.5 * year has a year coefficient lower by exactly that amount and identical fitted values, which the last two lines show. R 4.1.3 made model.weights() and model.offset() extract their components more carefully, and R 4.3.0 made predict.lm() look for the offset in the correct environment.
Updating formulas and fitted models
update() is a generic function with two ways to use it that share an idea. For a formula, update.formula(old, new) takes the new formula as a template: a dot on the left of the new formula stands for the left-hand side of the old one, and a dot on the right for its right-hand side. The result is simplified with terms.formula(simplify = TRUE), so terms that are added twice appear once and - removes a term; either argument may also be a length-one character vector, which is coerced with as.formula(), and the environment of the result is the environment of the old formula. For a fitted model, update(model, new.formula, ...) edits the call that created the model and evaluates it again, which is a convenient way to fit a model that differs from a previous one by a few added or removed terms (An Introduction to R: Updating fitted models).
source("scenario.R")
# `.` in the new formula means "what was there before", left and right.
f <- pages ~ year
update(f, . ~ . + genre)
#> pages ~ year + genre
update(f, log(.) ~ .)
#> log(pages) ~ year
update(f, ~ . + I(year^2))
#> pages ~ year + I(year^2)
update(f, . ~ . - 1)
#> pages ~ year - 1
update(pages ~ genre * year, . ~ . - genre:year)
#> pages ~ genre + year
update(f, sqrt(.) ~ . + genre)
#> sqrt(pages) ~ year + genre
# The right-hand side is simplified the way terms(simplify = TRUE) does.
update(pages ~ year + genre, . ~ . - year)
#> pages ~ genre
update(pages ~ year, . ~ . + year)
#> pages ~ year
# The new formula may be a string; it is coerced with as.formula().
update(f, ". ~ . + genre")
#> pages ~ year + genre
# The environment of the result is the environment of the old formula.
make <- function() pages ~ year
old <- make()
identical(environment(update(old, . ~ . + genre)), environment(old))
#> [1] TRUE
identical(environment(update(f, . ~ . + genre)), globalenv())
#> [1] TRUESource: R documentation: update.formula
The rule to remember is the one the manual and the help page share: the dot means whatever was there before, and only in this context. Everything on the right-hand side is passed through terms.formula(simplify = TRUE), which is why . ~ . + year on a formula that already has year does not add it again, and why the result shows genre + year when the interaction is removed from genre * year. Passing a string as the new formula works through as.formula(). For fitted models the code of update.default() takes the call stored in the model (getCall()), changes the arguments you name (and adds new ones, like subset below), and evaluates it in the parent frame; it has the argument evaluate = FALSE to return the call instead of running it. Since R 4.1.0 update.default() calls the generic update() on the formula, which makes it work for models with extended formulas.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
fit <- lm(pages ~ year, data = books)
fit$call
#> lm(formula = pages ~ year, data = books)
formula(fit)
#> pages ~ year
# update() on a fitted model edits the call and evaluates it again.
fit2 <- update(fit, . ~ . + genre)
fit2$call
#> lm(formula = pages ~ year + genre, data = books)
formula(fit2)
#> pages ~ year + genre
# Other arguments of the call can be changed or added in the same call.
fit3 <- update(fit, subset = year > 2000)
fit3$call
#> lm(formula = pages ~ year, data = books, subset = year > 2000)
nobs(fit)
#> [1] 13
nobs(fit3)
#> [1] 10
# evaluate = FALSE returns the new call without running it.
update(fit, . ~ . + genre, evaluate = FALSE)
#> lm(formula = pages ~ year + genre, data = books)
# The anova() comparison is the classic use.
anova(fit, update(fit, . ~ . + I(year^2)))
#> Analysis of Variance Table
#>
#> Model 1: pages ~ year
#> Model 2: pages ~ year + I(year^2)
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 11 88909
#> 2 10 85845 1 3064 0.3569 0.5635Source: R documentation: update
Because update() re-evaluates the call, it needs the variables the call refers to to be reachable again from the place where update() is called: the data frame books, a formula held in a variable and so on. A model that was fitted inside a function from a local data frame cannot be updated from the console, which is a common reason for an object 'x' not found error in an update() call. Refitting with lm(formula(fit), data = ...) and the data frame at hand is always available.
Building and inspecting formulas
A formula can be built from text, from term labels or by editing another formula. reformulate() is the base tool for the second case: it takes a character vector of term labels and, optionally, a response, and returns a formula in the environment of the caller. Its help page describes the arguments: the labels are joined with + when there are several, a zero-length vector gives an intercept-only (or empty) right-hand side, response is a character string, a symbol or a call, intercept = FALSE appends -1, and env sets the environment of the result (added in R 3.6.0). Non-syntactic names must be protected with backticks by the caller. Anything else that terms() can parse, such as "year * genre", can be a label.
source("scenario.R")
predictors <- c("year", "genre")
# reformulate() builds a formula from term labels; no string pasting.
reformulate(predictors, response = "pages")
#> pages ~ year + genre
reformulate(predictors)
#> ~year + genre
reformulate(predictors, response = "pages", intercept = FALSE)
#> pages ~ year + genre - 1
reformulate("year * genre", response = "pages")
#> pages ~ year * genre
# The response can be a symbol or a call.
reformulate("year", response = quote(log(pages)))
#> log(pages) ~ year
reformulate(predictors, response = as.name("pages"))
#> pages ~ year + genre
# Non-syntactic names are protected with backticks.
cols <- c("page count", "year")
tryCatch(reformulate(cols), error = \(e) "Syntax error.")
#> [1] "Syntax error."
reformulate(sprintf("`%s`", cols), response = "pages")
#> pages ~ `page count` + year
# The string-pasting way produces the same formula, with more room for typos.
identical(
as.formula(paste("pages ~", paste(predictors, collapse = " + "))),
reformulate(predictors, response = "pages")
)
#> [1] TRUE
# The environment is the caller's by default and can be set.
environmentName(environment(reformulate(predictors, "pages")))
#> [1] "R_GlobalEnv"
shared <- new.env()
identical(environment(reformulate(predictors, "pages", env = shared)), shared)
#> [1] TRUE
# Zero term labels give an intercept-only model.
reformulate(character(), response = "pages")
#> pages ~ 1Source: R documentation: reformulate
The older way of writing the same thing, as.formula(paste("y ~", paste(x, collapse = " + "))), still works and the script confirms that it produces an identical object. reformulate() avoids the pasting and the parsing of a response string; the help page records that a response that cannot be parsed still works for backward compatibility with a deprecation warning, that R 3.6.0 stopped reformulate("x", response = "sin(y)") producing extra backquotes, that R 4.1.3 gave a helpful error for a response of length greater than one, and that R 4.5.0 improved its handling together with drop.terms().
all.vars() and all.names() return the symbols in an expression (R documentation: all.names). all.names() returns every name, including the functions and operators, in order of appearance and with repeats; all.vars() returns the variables only, without functions, and without repeats unless unique = FALSE. With a formula all.vars(f) is the quick answer to "which columns does this model use?", and setdiff(all.vars(f), names(data)) finds the names that are not columns before model.frame() fails with a less direct message. Note that the dot is reported as a name, since it is a symbol until terms() expands it.
source("scenario.R")
f <- log(pages) ~ year + I(year^2) + genre
all.vars(f)
#> [1] "pages" "year" "genre"
all.names(f)
#> [1] "~" "log" "pages" "+" "+" "year" "I" "^" "year"
#> [10] "genre"
all.names(f, functions = FALSE)
#> [1] "pages" "year" "year" "genre"
all.vars(~ year)
#> [1] "year"
all.vars(pages ~ .)
#> [1] "pages" "."
# Typical use: check a formula against the data before fitting.
f <- pages ~ year + chapters
setdiff(all.vars(f), names(books))
#> [1] "chapters"
# Variables on one side only.
all.vars(f[[2]])
#> [1] "pages"
all.vars(f[[3]])
#> [1] "year" "chapters"
# Repeated variables are reported once unless unique = FALSE.
all.vars(pages ~ year + I(year^2))
#> [1] "pages" "year"
all.vars(pages ~ year + I(year^2), unique = FALSE)
#> [1] "pages" "year" "year"
# Select just the columns a formula uses.
head(books[intersect(names(books), all.vars(pages ~ year + genre))], 3)
#> year genre pages
#> 1 1974 Science fiction 387
#> 2 1813 Classic 432
#> 3 2019 Fiction 214A formula that is held in a variable and passed to lm() leaves a trace in the model’s call, and the next script shows the other edits: building formulas for a family of models, the pitfall of lm(formula = forms[[1L]], ...), and editing a formula as the language object it is. bquote() is the tool for a formula template whose parts are computed; the result is a call, not a formula, until it is evaluated, which is when ~ attaches the class and the environment. The detailed treatment of bquote(), substitute() and the manipulation of calls is in Language Objects and Quoting.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
# A formula built from character vectors, one model per predictor set.
predictor_sets <- c("year", "genre", "year + genre")
forms <- lapply(predictor_sets, \(p) reformulate(p, response = "pages"))
# Each formula remembers the frame of the anonymous function that built it, so
# printing it would add an environment line, so deparse1() shows the text.
vapply(forms, deparse1, character(1L))
#> [1] "pages ~ year" "pages ~ genre" "pages ~ year + genre"
r2 <- vapply(
forms,
\(f) summary(lm(f, data = books))$r.squared,
numeric(1L)
)
round(setNames(r2, predictor_sets), 3)
#> year genre year + genre
#> 0.257 0.770 0.814
# Pitfall: the call of a model fitted from a variable records the variable.
fit <- lm(forms[[1L]], data = books)
fit$call
#> lm(formula = forms[[1L]], data = books)
# ... so update() has to find `forms` and the data again when it re-runs it.
deparse1(formula(fit))
#> [1] "pages ~ year"
deparse1(formula(update(fit, . ~ . + genre)))
#> [1] "pages ~ year + genre"
# bquote() splices a response expression into a formula template. The result
# is a call until it is evaluated; evaluating attaches class and environment.
resp <- quote(log(pages))
cl <- bquote(.(resp) ~ year)
class(cl)
#> [1] "call"
f <- eval(cl)
class(f)
#> [1] "formula"
f
#> log(pages) ~ year
# Programmatic edits of an existing formula go through the language objects.
f2 <- pages ~ year + genre
rhs <- f2[[3]]
rhs
#> year + genre
class(rhs)
#> [1] "call"
rhs[[1]]
#> `+`
as.list(rhs)
#> [[1]]
#> `+`
#>
#> [[2]]
#> year
#>
#> [[3]]
#> genre
#>
f2[[3]] <- rhs[[2]]
f2
#> pages ~ yearFormulas outside modelling
The formula is a convenient piece of syntax for "this variable by these groups", and several functions of base R take it for that purpose, with data and sometimes subset and na.action arguments that work as in model.frame(). All of them build their model frame the same way, so the same rules apply: the variables are looked up in data and then in the environment of the formula, and rows with missing values are dropped by default.
-
aggregate(formula, data, FUN)appliesFUNto the left-hand variable(s) within the groups of the right-hand side.cbind(a, b) ~ ghas two responses; a dot on the left means all other columns; a dot on the right means all other columns as groups. Since R 4.2.0 the first argument of the formula method is calledx, like the generic’s, so portable code does not name it. -
xtabs(formula, data)builds contingency tables: a one-sided formula counts rows, a left-hand side sums that variable per cell. Missing values in the left-hand side make the sumNAunlessna.rm = TRUE(and since R 4.4.0xtabs(Freq ~ .)consistently usesna.pass); R 3.4.0 added theaddNAargument to count missing values. -
boxplot(formula, data)groups the response by the factor(s) on the right and draws one box per group. Withplot = FALSEit returns the statistics it would have drawn (a list withstats,n,conf,out,groupandnames) and opens no graphics device, which is the form used here; the drawing itself is in Base Graphics and Devices.plot(y ~ x)andplot(y ~ g)use the same formula for scatter plots and box plots. -
The formula methods of the classical tests:
t.test(pages ~ genre),wilcox.test(),kruskal.test(),bartlett.test(),var.test()andcor.test(~ a + b)(a one-sided formula with two variables); the paired tests usePair(x1, x2) ~ 1, which came with R 4.0.0, when the formula interface oft.test()andwilcox.test()was extended to one-sample and paired tests; since R 4.4.0 these formula methods reject apairedargument withcannot use 'paired' in formula method. R 4.3.0 fixed the formula method ofcor.test()whenenvironment(formula)was not the calling environment. -
sort_by(x, y)(new in R 4.4.0) sorts a data frame by one or more columns, and itsycan be a one-sided formula such as~year(R documentation: sort_by). -
lm(),glm(),aov(),nls()and the rest of the modelling functions, covered in Statistical Models.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
# aggregate(): response ~ groups, one summary per group.
aggregate(pages ~ genre, data = books, FUN = mean)
#> genre pages
#> 1 classic 356.0000
#> 2 fiction 222.6667
#> 3 mystery 272.5000
#> 4 non-fiction 263.5000
#> 5 poetry 92.0000
#> 6 science fiction 373.0000
# Two responses with cbind(), and two grouping variables with +.
aggregate(cbind(pages, year) ~ genre, data = books, FUN = median)
#> genre pages year
#> 1 classic 356.0 1906
#> 2 fiction 214.0 2021
#> 3 mystery 272.5 2017
#> 4 non-fiction 263.5 2016
#> 5 poetry 88.0 2005
#> 6 science fiction 387.0 2014
books$recent <- books$year >= 2015
aggregate(pages ~ genre + recent, data = books, FUN = length)
#> genre recent pages
#> 1 classic FALSE 2
#> 2 non-fiction FALSE 1
#> 3 poetry FALSE 1
#> 4 science fiction FALSE 2
#> 5 fiction TRUE 3
#> 6 mystery TRUE 2
#> 7 non-fiction TRUE 1
#> 8 science fiction TRUE 1
# `.` on the left means "all the other columns"; rows with an NA in any used
# column are dropped first, because the default na.action is na.omit.
aggregate(. ~ genre, data = books[c("pages", "year", "genre")], FUN = mean)
#> genre pages year
#> 1 classic 356.0000 1906.000
#> 2 fiction 222.6667 2021.000
#> 3 mystery 272.5000 2017.000
#> 4 non-fiction 263.5000 2016.000
#> 5 poetry 88.0000 2005.000
#> 6 science fiction 373.0000 2001.667
# xtabs(): contingency tables; a one-sided formula counts rows.
xtabs(~genre, data = books)
#> genre
#> classic fiction mystery non-fiction poetry
#> 2 4 2 3 2
#> science fiction
#> 3
xtabs(~ genre + recent, data = books)
#> recent
#> genre FALSE TRUE
#> classic 2 0
#> fiction 0 4
#> mystery 0 2
#> non-fiction 2 1
#> poetry 1 0
#> science fiction 2 1
# With a left-hand side it adds that variable up per cell. Two genres have a
# book without a page count, so their sums are NA, which print() shows blank.
xtabs(pages ~ genre, data = books)
#> genre
#> classic fiction mystery non-fiction poetry
#> 712 545 184
#> science fiction
#> 1119
xtabs(pages ~ genre, data = books, na.rm = TRUE)
#> genre
#> classic fiction mystery non-fiction poetry
#> 712 668 545 527 184
#> science fiction
#> 1119
# boxplot(): the formula splits the response by the group. plot = FALSE returns
# the numbers without opening a graphics device.
b <- boxplot(pages ~ genre, data = books, plot = FALSE)
names(b)
#> [1] "stats" "n" "conf" "out" "group" "names"
b$names
#> [1] "classic" "fiction" "mystery" "non-fiction"
#> [5] "poetry" "science fiction"
b$n
#> [1] 2 3 2 2 2 3
b$stats
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 280 198 244.0 175.0 88 322.0
#> [2,] 280 206 244.0 175.0 88 354.5
#> [3,] 356 214 272.5 263.5 92 387.0
#> [4,] 432 235 301.0 352.0 96 398.5
#> [5,] 432 256 301.0 352.0 96 410.0
# One-sided and two-sided formulas in tests and in sorting.
cor.test(~ pages + year, data = books)$estimate
#> cor
#> -0.5069084
sort_by(books[c("title", "year")], ~year)[1:3, ]
#> title year
#> 2 Pride and Prejudice 1813
#> 1 The Dispossessed 1974
#> 8 The Lamp at Noon 1999
two <- books[books$genre %in% c("fiction", "mystery"), ]
unname(t.test(pages ~ genre, data = two)$statistic)
#> [1] -1.494832The two sums that print as blanks come from the NA in pages for one fiction and one non-fiction book; print() of an xtabs table leaves NA empty, which is easy to read as zero, so check for missing values or use na.rm = TRUE on purpose. The first aggregate() calls also show an effect of the default na.action = na.omit: the means are of the non-missing books only, and the counts in the length call add up to 13, not to the 16 rows of books, because pages is missing for two books and recent, computed from year, is NA for the book without a year. A boxplot summary is the five numbers of stats per group: lower whisker, lower hinge, median, upper hinge and upper whisker.
Packages repeat the convention and extend it. lattice (a recommended package that ships with R, 0.22-9) uses y ~ x | g to draw one panel per group, and ggplot2 4.0.3 uses facet_wrap(~ genre) with a one-sided formula (ggplot2 package documentation: facet_wrap); neither is run on this page. purrr 1.2.2 accepts a one-sided formula such as ~ .x + 2 as a short anonymous function; its documentation marks that form as "no longer recommended" (purrr package documentation: as_mapper), and since R 4.1.0 the base lambda \(x) x + 2 does the same without a package. That use of the tilde has nothing to do with models: the formula is only a quoted expression that the package turns into a function, a use of the same "quote, keep the environment" mechanism described at the top of this page. The tidyverse side of formulas is in tidyr, purrr and stringr.
Pitfalls
Most formula problems come from a few places. Two of them have entries in the R FAQ. The model without an intercept changes the definition of R squared, so the two numbers in the output are not comparable, and the sequential sums of squares of anova() depend on the order of the terms.
source("scenario.R")
with_intercept <- lm(pages ~ year, data = books)
through_origin <- lm(pages ~ year - 1, data = books)
identical(coef(lm(pages ~ 0 + year, data = books)), coef(through_origin))
#> [1] TRUE
# summary() uses a different formula for R squared when there is no intercept
# (R FAQ 7.41), so the two numbers are not comparable.
round(summary(with_intercept)$r.squared, 3)
#> [1] 0.257
round(summary(through_origin)$r.squared, 3)
#> [1] 0.887
colnames(model.matrix(through_origin))
#> [1] "year"The FAQ explains that with a zero intercept summary.lm() compares the fit with a model whose mean is zero, not with the mean of the response, because otherwise R squared could be negative; if you do not know a priori that the response is zero when the predictor is zero, you probably should not be fitting a line through the origin. The second entry is about the order of the terms.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
# anova() reports sequential sums of squares, so the order of the terms in
# the formula changes the table when the data are unbalanced (R FAQ 7.18).
anova(lm(pages ~ year + genre, data = books))
#> Analysis of Variance Table
#>
#> Response: pages
#> Df Sum Sq Mean Sq F value Pr(>F)
#> year 1 30746 30746.1 8.310 0.02796 *
#> genre 5 66710 13342.0 3.606 0.07495 .
#> Residuals 6 22199 3699.9
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
anova(lm(pages ~ genre + year, data = books))
#> Analysis of Variance Table
#>
#> Response: pages
#> Df Sum Sq Mean Sq F value Pr(>F)
#> genre 5 84854 16970.7 4.5868 0.04545 *
#> year 1 12602 12602.3 3.4061 0.11449
#> Residuals 6 22199 3699.9
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# terms() keeps the order of the main effects as written.
attr(terms(pages ~ year + genre), "term.labels")
#> [1] "year" "genre"
attr(terms(pages ~ genre + year), "term.labels")
#> [1] "genre" "year"Because the data are unbalanced (and the terms correlated), the year line of the first table tests year alone and the year line of the second tests year after genre; the sums of squares and the p-values differ. The FAQ points out that anova(model1, model2) compares any two models directly and drop1(model1) shows the sums of squares of dropping each term once. The last two lines show that terms() keeps the main effects in the order written, which is why the order in the formula is the order of the sequential table.
The next script shows the dropped rows at work. The two models are fitted on 14 and 13 rows, so anova() refuses to compare them and gives the real message; fitting both on the complete rows fixes it; na.exclude keeps a slot for each dropped row; and a factor level without data stays in the model frame and the model matrix, while lm() drops it.
source("scenario.R")
books$genre <- tolower(trimws(books$genre))
# Two models on different rows cannot be compared with anova().
small <- lm(pages ~ genre, data = books)
large <- lm(pages ~ genre + year, data = books)
c(small = nobs(small), large = nobs(large))
#> small large
#> 14 13
mismatch <- tryCatch(
anova(small, large),
error = \(e) conditionMessage(e)
)
mismatch
#> [1] "models were not all fitted to the same size of dataset"
# Fitting both on the rows that have every variable makes them comparable.
complete <- books[complete.cases(books[c("pages", "genre", "year")]), ]
anova(
lm(pages ~ genre, data = complete),
lm(pages ~ genre + year, data = complete)
)
#> Analysis of Variance Table
#>
#> Model 1: pages ~ genre
#> Model 2: pages ~ genre + year
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 7 34802
#> 2 6 22199 1 12602 3.4061 0.1145
# na.exclude keeps a slot for each dropped row in residuals() and predictions.
fit <- lm(pages ~ year, data = books, na.action = na.exclude)
length(residuals(fit))
#> [1] 16
sum(is.na(residuals(fit)))
#> [1] 3
# An unused factor level is kept by model.frame() and model.matrix() as a
# column of zeros; lm() asks model.frame() to drop unused levels.
two <- books[books$genre %in% c("classic", "poetry"), ]
two$genre <- factor(two$genre, levels = c("classic", "poetry", "mystery"))
colSums(model.matrix(~genre, data = two))
#> (Intercept) genrepoetry genremystery
#> 4 2 0
nlevels(model.frame(pages ~ genre, data = two)$genre)
#> [1] 3
nlevels(model.frame(pages ~ genre, data = two, drop.unused.levels = TRUE)$genre)
#> [1] 2
coef(lm(pages ~ genre, data = two))
#> (Intercept) genrepoetry
#> 356 -264Source: R documentation: model.frame
The other pitfalls, with what R does:
-
Missing values vanish silently. The default
na.omitdrops every row with anNAin any variable in the model, so a model with an extra covariate can be fitted on fewer rows than the one without it, andanova(fit1, fit2)then fails with the messagemodels were not all fitted to the same size of dataset, as in the script above.nobs(fit)andfit$na.actiontell you what happened;na.action = na.excludekeeps a slot, filled withNA, for each dropped row in the residuals and predictions. -
The power operator.
y ~ x^2is not a square; useI(x^2)(R FAQ 3.3.2). In the same wayy ~ x + zis two terms,y ~ I(x + z)is one, andy ~ log(x + 1)is fine because the operator is inside a function call. -
Variables that are not columns. A name that is not a column of
datais looked up in the environment of the formula (the global environment for a console formula), which may be a different variable of the same name than you meant. When a function builds a formula, pass the data explicitly, keep the variables as columns, and avoid formulas that depend on a variable defined in the caller; if the formula is passed on to another function,environment(f)is the thing to check. -
The formula travels with its environment. A large environment is kept alive and serialized with every object that stores the formula, including fitted models;
lm()objects from a function that held a big data frame can become very large for that reason, and one remedy is to build the formula withenv = globalenv()or to assignenvironment(f) <- globalenv(). -
Data-dependent terms.
poly(),scale(),bs()andmean(x)inside a formula are computed on the full data beforesubsetis applied, andpoly()rejectsNA. -
Names with spaces or operators. Use backticks in the formula and in
reformulate(), or rename the column withmake.names()first. -
Factors with unused levels.
model.frame()keeps a level with no observations (drop.unused.levels = FALSEis its default), somodel.matrix()gives it a column of zeros;lm()andglm()askmodel.frame()to drop unused levels, and other code can usedroplevels(). -
The character vector of length greater than one passed to
formula()is deprecated (R 4.0.0): build one string withpaste(x, collapse = " ").
Version notes
The formula machinery is stable and old; the changes since R 3.4 are mostly fixes and a few additions. The entries come from the NEWS files of the corresponding releases (R NEWS).
| Version | Change |
|---|---|
R 3.4.0 |
|
R 3.6.0 |
|
R 4.0.0 |
|
R 4.1.0 |
|
R 4.1.3 |
|
R 4.2.0 |
The formula method of |
R 4.3.0 |
|
R 4.3.3 |
|
R 4.4.0 |
New generic |
R 4.5.0 |
|
R 4.5.3 |
|
R 4.6.0 |
|
The entries for R 3.4.0 and 3.6.0 are in the R 3.x NEWS (R NEWS: Changes in R 3.x). The R-devel NEWS (the upcoming 4.7.0, never run here) says that model.frame.glm() now uses the same logic as model.frame.lm(); it lists no change to the formula language itself at the time of writing.
Where older tutorials differ
| Older tutorial writes | Current practice | Why | Since |
|---|---|---|---|
|
|
|
always |
|
|
no pasting, no parsing of a response string, and an |
R 3.6.0 ( |
|
|
the vector form silently truncated bad input and is deprecated |
R 4.0.0 |
|
|
a formula or a list of columns reads as the intent |
R 4.4.0 |
|
|
purrr marks the formula lambda "no longer recommended" and the base lambda needs no package |
R 4.1.0 (lambda) |
|
|
a global option changes every later model in the session, including those in other packages |
always |
|
treatment contrasts (the default), with |
R and S have different defaults for unordered factors; the manual calls R’s default a deliberate difference |
always |
|
|
|
always |
|
refit with |
|
always |
A formula built with |
|
|
R 3.6.0 ( |
Source: R FAQ 3.3.2: Models
Where to go next
-
Statistical Models for fitting
lm()andglm(), readingsummary(),anova()and diagnostics, and usingpredict()on new data. -
Language Objects and Quoting for formulas as calls:
quote(),bquote(),substitute(),match.call()and editing a call as a list. -
Evaluation and Non-Standard Evaluation for
eval(), data masks (model.frame()is the formula version of the same idea) andwith(). -
Scoping, Environments and Closures for the environment a formula carries and
parent.frame(). -
Data Frames for
aggregate(),xtabs(),sort_by()andsubset()on a data frame. -
Factors for levels,
droplevels(),relevel()and the ordered factors that get polynomial contrasts. -
Descriptive Statistics and Tests for
t.test(),cor.test()and the other formula interfaces of the classical tests. -
Base Graphics and Devices and Grid Graphics and Lattice for
plot(y ~ x),boxplot()and the lattice formula with|. -
Attributes and Names for the
classand.Environmentattributes that make the formula, andtermsattributes.
References
-
An Introduction to R: 11 Statistical models in R; 11.1 Defining statistical models; formulae, 11.1.1 Contrasts, 11.2 Linear models, 11.3 Generic functions for extracting model information, 11.5 Updating fitted models
-
The R Language Definition: Operators; The call stack; Computing on the language
-
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.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 3.4.0 to 4.6.1; R 3.x NEWS; R-devel NEWS
-
contr.treatment, contr.sum, contr.helmert, contr.poly and contr.SAS