1.4.4 The formula algebra
R’s formulas have an algebra of terms, and so do these. A term is a set of variables, the table’s columns: a column name is a term of one variable, and an interaction of several variables is a term whose column in the design matrix is the product of theirs. The operators build terms from terms, as R’s terms builds them, and the design matrix has a column per term, as R’s model.matrix builds it.
1.4.4.1 Interactions and crossing
(: age dose), R’s age:dose, is the interaction of age and dose, and (* age dose), R’s age*dose, crosses them: both main effects and their interaction. The interaction’s column is named as R names it, with a colon:
> (formula-predictor-names (~ bp (* age dose)) patients) '("age" "dose" "age:dose")
> (define X (formula-design-matrix (~ bp (* age dose)) patients)) > (design-matrix->rows (design-matrix-select-rows X '(0 1)))
'((55.0 4.512962181748854 248.21291999618697)
(73.0 0.03958862699438688 2.889969770590242))
Each row’s third entry is the product of the first two. The fit is the matrix fit of that design matrix. In the simulated data blood pressure does not depend on the product, and the lasso leaves it out:
> (coef (formula-fit (~ bp (* age dose)) patients #:lambda 0.5))
'(("(Intercept)" . 93.51132211478343)
("age" . 0.5139586017419673)
("dose" . -1.7735029027253124)
("age:dose" . 0.0))
The operators take any terms, so sums cross and interact term by term: (: (+ age dose) marker) gives "age:marker" and "dose:marker", and (* (+ age dose) marker) adds the three main effects. The prefix forms take more than two operands, folded from the left, and (* age dose marker) crosses all three:
> (formula-predictor-names (~ bp (: (+ age dose) marker)) patients) '("age:marker" "dose:marker")
> (formula-predictor-names (~ bp (* age dose marker)) patients)
'("age"
"dose"
"marker"
"age:dose"
"age:marker"
"dose:marker"
"age:dose:marker")
The terms come out in R’s order: main effects first, then the interactions of two variables, then of three, each in the order they first appear. A term written twice, as (: age dose) and (: dose age) are, counts once.
1.4.4.2 Powers, and R’s x^2
(^ term n), R’s (a + b + c)^n, crosses a sum with itself: every interaction of at most n of its terms. With 2, that is the main effects and every two-way interaction:
> (formula-predictor-names (~ bp (^ (+ age dose marker) 2)) patients) '("age" "dose" "marker" "age:dose" "age:marker" "dose:marker")
A variable crossed with itself is itself, so in a formula x^2 is not a square. R’s y ~ x + x^2 is the same model as y ~ x, and so is (~ bp age (^ age 2)):
> (formula-predictor-names (~ bp age (^ age 2)) patients) '("age")
This surprises people who come to R’s formulas from algebra. R writes a square I(x^2), a transform, and so does a formula here: (I (expt age 2)), or (sqr age) (see Squares, and R’s I()).
1.4.4.3 Removing terms
(- term removed ...) removes the terms that are equal to a term of removed, so it can take an interaction out of a crossing or a column out of all. A term that is not there is ignored, as in R:
> (formula-predictor-names (~ bp (- (* age dose marker) (: age dose marker))) patients) '("age" "dose" "marker" "age:dose" "age:marker" "dose:marker")
> (formula-predictor-names (~ bp (- (* age dose) dose)) patients) '("age" "age:dose")
> (formula-predictor-names (~ bp (- (* age dose) marker)) patients) '("age" "dose" "age:dose")
1.4.4.4 The intercept
A model has an intercept unless its formula says otherwise. 0, or removing 1 with - 1, fits without one, as R’s y ~ 0 + x and y ~ x - 1 do; 1 keeps it. The coefficients still list the intercept, at zero, as R’s coef does:
> (coef (formula-fit (~ bp 0 age dose) patients #:lambda 0.5))
'(("(Intercept)" . 0.0)
("age" . 1.9427491619268988)
("dose" . 0.6144740506407587))
The fit procedures’ #:intercept? then defaults to the formula’s intercept. Given explicitly, it must agree with a formula that writes 1, 0 or - 1, and a contradiction is an error that names both; with a formula that writes none, it decides. The Cox family has no intercept, and accepts 1, 0 and - 1 with any #:intercept?. As in R, 0 and - 1 still change how a factor is coded (see Treatment contrasts).
A formula whose terms leave no predictors, such as (~ bp 1) or (~ bp), is R’s intercept-only model. glmnet cannot fit it, since it needs at least one predictor, and the formula procedures say so.
1.4.4.5 Infix formulas
The operators can also go between terms, as R writes them. The Racket reader’s infix dots make (bp . ~ . age * dose + marker) the same as (~ bp age * dose + marker), so an R formula reads almost as written:
> (define f (bp . ~ . 1 + age * dose + marker)) > f (~ bp 1 + age * dose + marker)
> (formula-predictor-names f patients) '("age" "dose" "marker" "age:dose")
The formula prints as the ~ form that the reader makes of it. The operators have R’s precedence: ^ binds tightest, then a leading - or +, then :, then *, then + and -, and each groups from the left, except that a power cannot be raised again, as in R. A parenthesized group can be infix too, and infix and prefix forms mix:
> (formula-predictor-names (bp . ~ . (age + dose + marker) ^ 2 - age : dose) patients) '("age" "dose" "marker" "age:marker" "dose:marker")
> (formula-predictor-names (bp . ~ . - 1 + (* age dose) + marker) patients) '("age" "dose" "marker" "age:dose")
An operator needs spaces around it. The reader reads age:dose and -age as one name each, -1 as a number, and -0 as 0, dropping the sign that R reads as an operator, and a formula that uses one is a syntax error that says so:
> (~ bp age:dose) eval:56:0: ~: age:dose reads as one name; put spaces around
an operator, as in wt : hp, or write a column with this name
as a string, "age:dose"
at: age:dose
in: (~ bp age:dose)
parsing context:
while parsing a column name
term: age:dose
location: eval:56:0
while parsing a term: a column, all, 1, 0, a transform
such as (log x), (factor x), or a group such as (* x z) or
(x + z)
term: age:dose
location: eval:56:0
while parsing the right-hand side of a formula
term: (age:dose)
location: eval:56:0
R’s -0, which in y ~ x + -0 keeps the intercept, is therefore written (- 0), as in (y . ~ . x + (- 0)).
1.4.4.6 The response on the right-hand side
A response column is not a predictor. R’s model.matrix drops the response when the right-hand side has it as a term of its own, with a warning, and keeps it in an interaction; so do the formula procedures, which log the warning on the 'glmnet topic. The columns of a Cox or multi-response response are each dropped in the same way:
> (formula-predictor-names (bp . ~ . age + bp + age : bp) patients) '("age" "bp:age")
1.4.4.7 Predicting from new data
A formula model keeps the terms it was fitted with. predict builds their design matrix from a new table’s columns, which can come in any order, and does not expand the formula again, so all stands for the columns it stood for in the fit. The table needs only the columns that the terms read, here the age and the dose:
> (define crossed (formula-fit (~ bp (* age dose)) patients #:lambda 0.1)) > (predict crossed (list (cons "dose" '(2.0 8.0)) (cons "age" '(40 70)))) '(110.48450167523144 115.28077863947216)
A table that lacks one of them is an error that names it, as A Gaussian fit shows.
1.4.4.8 From R’s syntax
R |
| Prefix |
| Infix |
y ~ a + b |
| (~ y a b) |
| |
y ~ a:b |
| (~ y (: a b)) |
| (~ y a : b) |
y ~ a*b |
|
| ||
y ~ (a + b + c)^2 |
|
| ||
y ~ a*b - a |
|
| ||
y ~ 0 + a |
| (~ y 0 a) |
| |
y ~ a - 1 |
|
| ||
y ~ . |
| (~ y all) |
| |
y ~ . - a |
|
| ||
y ~ .^2 |
| (~ y (^ all 2)) |
| (~ y all ^ 2) |
Surv(time, status) ~ . |
|
| ||
cbind(y1, y2) ~ a |
| (~ (y1 y2) a) |
|
Each infix form can also be written with the reader’s infix dots, as (y . ~ . a * b). Transforms has the transforms, such as R’s log(a), From R’s factors the factors, such as R’s factor(a), and Matching R lists what the formula language does not have.