Formulas for linear models

Basic form

  • In R linear models are specified using a model formula, which is an expression that contains a tilde (the ~ character).
  • The response is on the left-hand side of ~, typically as the name of a variable, e.g. optden, but it can also be a function call, e.g. log(BrainWt).
  • The right-hand side of the formula is composed of model terms separated by +.

Fitting linear models is not trivial

  • Although it may seem straightforward, fitting a linear model can be quite involved.
  • Complications arrive from
    • missing data
    • categorical covariates
    • terms that are function calls
    • several possible auxiliary arguments such as subset, na.action, constrasts, …
    • provision for analysis of variance and other testing procedures

Model frames and model matrices

  • Because the same issues are encountered in any statistical model that is based on a model matrix, there is a standard approach.
    • Create a model.frame by examining the formula/data specification, applying a subset specification, handling NAs, evaluating function calls, re-ordering terms if necessary.
    • Create a model.matrix and model_response from the model frame. Associate columns in the model matrix with terms in the model frame.
  • Note the use of the term “model matrix”. Sometimes $\bf X$ in $\bf X\beta$ is called a “design matrix” but that is a misnomer.

A simple example - Formaldehyde

  • The Formaldehyde data in the datasets package is from a calibration experiment in which the optical density in an assay is measure for various carbohydrate concentrations.
str(Formaldehyde)

## 'data.frame':	6 obs. of  2 variables:
##  $ carb  : num  0.1 0.3 0.5 0.6 0.7 0.9
##  $ optden: num  0.086 0.269 0.446 0.538 0.626 0.782

Formaldehyde data plot

plot of chunk formplot

Fitting a simple linear regression model 

summary(m1 <- lm(optden ~ 1 + carb, Formaldehyde))

## 
## Call:
## lm(formula = optden ~ 1 + carb, data = Formaldehyde)
## 
## Residuals:
##         1         2         3         4         5         6 
## -0.006714  0.001029  0.002771  0.007143  0.007514 -0.011743 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.005086   0.007834   0.649    0.552
## carb        0.876286   0.013535  64.744 3.41e-07
## 
## Residual standard error: 0.008649 on 4 degrees of freedom
## Multiple R-squared:  0.999,	Adjusted R-squared:  0.9988 
## F-statistic:  4192 on 1 and 4 DF,  p-value: 3.409e-07

Extractor methods 

class(m1)

## [1] "lm"
  • There are many extractor methods defined for this class. Use methods(class="lm") to see them.

A model frame contains information on terms 

str(model.frame(m1))

## 'data.frame':	6 obs. of  2 variables:
##  $ optden: num  0.086 0.269 0.446 0.538 0.626 0.782
##  $ carb  : num  0.1 0.3 0.5 0.6 0.7 0.9
##  - attr(*, "terms")=Classes 'terms', 'formula'  language optden ~ 1 + carb
##   .. ..- attr(*, "variables")= language list(optden, carb)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "optden" "carb"
##   .. .. .. ..$ : chr "carb"
##   .. ..- attr(*, "term.labels")= chr "carb"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: 0xb819d9e00> 
##   .. ..- attr(*, "predvars")= language list(optden, carb)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. ..- attr(*, "names")= chr [1:2] "optden" "carb"

Model matrix relates columns to terms 

model.matrix(m1)

##   (Intercept) carb
## 1           1  0.1
## 2           1  0.3
## 3           1  0.5
## 4           1  0.6
## 5           1  0.7
## 6           1  0.9
## attr(,"assign")
## [1] 0 1

lm objects have many components 

str(m1)

## List of 12
##  $ coefficients : Named num [1:2] 0.00509 0.87629
##   ..- attr(*, "names")= chr [1:2] "(Intercept)" "carb"
##  $ residuals    : Named num [1:6] -0.00671 0.00103 0.00277 0.00714 0.00751 ...
##   ..- attr(*, "names")= chr [1:6] "1" "2" "3" "4" ...
##  $ effects      : Named num [1:6] -1.12146 0.55996 0.00514 0.00992 0.01069 ...
##   ..- attr(*, "names")= chr [1:6] "(Intercept)" "carb" "" "" ...
##  $ rank         : int 2
##  $ fitted.values: Named num [1:6] 0.0927 0.268 0.4432 0.5309 0.6185 ...
##   ..- attr(*, "names")= chr [1:6] "1" "2" "3" "4" ...
##  $ assign       : int [1:2] 0 1
##  $ qr           :List of 5
##   ..$ qr   : num [1:6, 1:2] -2.449 0.408 0.408 0.408 0.408 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:6] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "(Intercept)" "carb"
##   .. ..- attr(*, "assign")= int [1:2] 0 1
##   ..$ qraux: num [1:2] 1.41 1.15
##   ..$ pivot: int [1:2] 1 2
##   ..$ tol  : num 1e-07
##   ..$ rank : int 2
##   ..- attr(*, "class")= chr "qr"
##  $ df.residual  : int 4
##  $ xlevels      : Named list()
##  $ call         : language lm(formula = optden ~ 1 + carb, data = Formaldehyde)
##  $ terms        :Classes 'terms', 'formula'  language optden ~ 1 + carb
##   .. ..- attr(*, "variables")= language list(optden, carb)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "optden" "carb"
##   .. .. .. ..$ : chr "carb"
##   .. ..- attr(*, "term.labels")= chr "carb"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: 0xb819d9e00> 
##   .. ..- attr(*, "predvars")= language list(optden, carb)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. ..- attr(*, "names")= chr [1:2] "optden" "carb"
##  $ model        :'data.frame':	6 obs. of  2 variables:
##   ..$ optden: num [1:6] 0.086 0.269 0.446 0.538 0.626 0.782
##   ..$ carb  : num [1:6] 0.1 0.3 0.5 0.6 0.7 0.9
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language optden ~ 1 + carb
##   .. .. ..- attr(*, "variables")= language list(optden, carb)
##   .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:2] "optden" "carb"
##   .. .. .. .. ..$ : chr "carb"
##   .. .. ..- attr(*, "term.labels")= chr "carb"
##   .. .. ..- attr(*, "order")= int 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: 0xb819d9e00> 
##   .. .. ..- attr(*, "predvars")= language list(optden, carb)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. .. ..- attr(*, "names")= chr [1:2] "optden" "carb"
##  - attr(*, "class")= chr "lm"

Other extractors

  • model.response goes through model.frame
  • fitted and residuals just extract components
  • the qr component provides a QR decompositon of $\bf X$
  • the effects component, $\bf Q’\bf y$, is used for analysis of variance, anova
  • many extractors such as hatvalues, cooks.distance, dfbeta, dfbetas, rstandard, rstudent and kappa are related to regression diagnostics.

Anova and effects 

effects(m1)^2

##  (Intercept)         carb                                                     
## 1.257668e+00 3.135496e-01 2.641354e-05 9.832048e-05 1.143173e-04 6.014870e-05 
## attr(,"assign")
## [1] 0 1
## attr(,"class")
## [1] "coef"

anova(m1)

## Analysis of Variance Table
## 
## Response: optden
##           Df   Sum Sq  Mean Sq F value    Pr(>F)
## carb       1 0.313550 0.313550  4191.8 3.409e-07
## Residuals  4 0.000299 0.000075

Formula examples

Covariate names in the examples

  • In the formulas below we write the response as y, continuous covariates as x,z,u,... and categorical covariates as f and g.
  • Note that the categorical covariates are assumed to be stored as factors

Simple linear regression 

y ~ x

denotes the simple linear regression model \(y_i = \beta_0+\beta_1 x_i + \epsilon_i, \quad i = 1,\dots,n\)

  • In this formula, the intercept term is implicit. To make it explicit use
y ~ 1 + x
  • For clarity I prefer to use the explicit form.
  • We are still debating the formula terms for the Julia language. Right now the prevailing opinion is not to assume an implicit intercept term.

Regression through the origin

  • To suppress an intercept term use
y ~ 0 + x
  • An alternative formula is
y ~ x - 1
  • Both of these generate the model \(y_i = \beta_1 x_i + \epsilon_i, \quad i = 1,\dots, n\)

Zero intercept for Formaldehyde? 

m2 <- lm(optden ~ 0 + carb, Formaldehyde)
anova(m2,m1)

## Analysis of Variance Table
## 
## Model 1: optden ~ 0 + carb
## Model 2: optden ~ 1 + carb
##   Res.Df        RSS Df  Sum of Sq      F Pr(>F)
## 1      5 0.00033073                            
## 2      4 0.00029920  1 3.1526e-05 0.4215 0.5516

coef(summary(m1))

##                Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept) 0.005085714 0.007833679  0.6492115 5.515953e-01
## carb        0.876285714 0.013534536 64.7444207 3.409192e-07

Multiple linear regression

  • Multiple covariate terms can be given, as in
y ~ 1 + x + z + sqrt(u)

corresponding to \(y_i = \beta_0 + \beta_i x_i+\beta_2 z_i+\beta_3 \sqrt{u_i}+\epsilon_i, \quad i=1,\dots,n\)

Polynomial regression

  • To include polynomial terms you must protect the ^ operator by wrapping it in I().
  • That is, the model \(y_i=\beta_0+\beta_1 x_i+\beta_2 x_i^2 + \beta_3 x_i^3+\epsilon_i,\quad i=1,\dots,n\) is written
y ~ 1 + x + I(x^2) + I(x^3)

Orthogonal polynomial terms

  • Another specification for a polynomial regression model uses the poly() function which generates orthogonal polynomial terms
y ~ poly(x, 3)
  • The fitted responses will be the same from the model using I(x^2), etc. but the coefficients will be different.
  • Orthogonal polynomials allow for backward elimination of higher order polynomial terms without refitting the model. (Not a consideration these days)
  • Fitting high order polynomials is discouraged. Smoothing approaches are preferred.

Quadratic for Formaldehyde? 

m3 <- lm(optden ~ 1 + carb + I(carb^2), Formaldehyde)
coef(summary(m3))

##                Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept) -0.01168268 0.006001455 -1.946642 1.467548e-01
## carb         0.97112338 0.026782523 36.259593 4.613325e-05
## I(carb^2)   -0.09621212 0.026309447 -3.656942 3.532047e-02

m4 <- lm(optden ~ poly(carb,2), Formaldehyde)
coef(summary(m4))

##                   Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept)     0.45783333 0.001745171 262.342918 1.221350e-07
## poly(carb, 2)1  0.55995503 0.004274779 130.990397 9.809825e-07
## poly(carb, 2)2 -0.01563262 0.004274779  -3.656942 3.532047e-02

A single categorical covariate

  • A one-way analysis of variance model for the levels of a factor, f, is written
y ~ 1 + f

  • Often we use the function aov() to fit such a model instead of lm(). The result is the same except for the class which changes the way that the fitted model is summarized.

  • The model being fit is sometimes written \(y_{ij}=\mu+\alpha_i+\epsilon_{ij},\quad i=1,\dots,I\;j=1,\dots,n_i\) although it is not fit in that form.

InsectSprays 

str(InsectSprays)

## 'data.frame':	72 obs. of  2 variables:
##  $ count: num  10 7 20 14 14 12 10 23 17 20 ...
##  $ spray: Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...

plot of chunk insectspraysplot1

More InsectSprays plots

plot of chunk insectspraysplot2 plot of chunk insectspraysplot3 plot of chunk insectspraysplot4

InsectSprays fit 

summary(m5 <- aov(sqrt(count) ~ 1 + spray, InsectSprays))

##             Df Sum Sq Mean Sq F value Pr(>F)
## spray        5  88.44  17.688    44.8 <2e-16
## Residuals   66  26.06   0.395

effects(m5)[1:6]^2

## (Intercept)      sprayB      sprayC      sprayD      sprayE      sprayF 
## 569.5041501  16.3082599  27.5296750   9.5833274  34.6174157   0.3991945

c(sum(effects(m5)[2:6]^2), sum(effects(m5)[-(1:6)]^2))

## [1] 88.43787 26.05798

Contrasts and model matrices 

unique(model.matrix(m5))

##    (Intercept) sprayB sprayC sprayD sprayE sprayF
## 1            1      0      0      0      0      0
## 13           1      1      0      0      0      0
## 25           1      0      1      0      0      0
## 37           1      0      0      1      0      0
## 49           1      0      0      0      1      0
## 61           1      0      0      0      0      1

contrasts(InsectSprays$spray)

##   B C D E F
## A 0 0 0 0 0
## B 1 0 0 0 0
## C 0 1 0 0 0
## D 0 0 1 0 0
## E 0 0 0 1 0
## F 0 0 0 0 1

Terms in m5 

terms(m5)

## sqrt(count) ~ 1 + spray
## attr(,"variables")
## list(sqrt(count), spray)
## attr(,"factors")
##             spray
## sqrt(count)     0
## spray           1
## attr(,"term.labels")
## [1] "spray"
## attr(,"order")
## [1] 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: 0xb819d9e00>
## attr(,"predvars")
## list(sqrt(count), spray)
## attr(,"dataClasses")
## sqrt(count)       spray 
##   "numeric"    "factor"

Two categorical factors

  • The additive model is specified as
y ~ 1 + f + g
  • The analysis of variance for a model with multiple categorical covariates is the sequential anova.
  • For unbalanced data the F statistics for f and g change if you specify the model as y ~ 1 + g + f
  • The two-factor model with interactions is
y ~ 1 + f + g + f:g

or, equivalently,

y ~ 1 + f * g

Examples of additive models 

summary(m6 <- aov(price ~ 1 + carat + cut + color, diamonds))

##                Df    Sum Sq   Mean Sq  F value Pr(>F)
## carat           1 7.291e+11 7.291e+11 355472.6 <2e-16
## cut             4 6.133e+09 1.533e+09    747.5 <2e-16
## color           6 1.260e+10 2.100e+09   1023.7 <2e-16
## Residuals   53928 1.106e+11 2.051e+06

summary(m7 <- aov(price ~ 1 + carat + color + cut, diamonds))

##                Df    Sum Sq   Mean Sq  F value Pr(>F)
## carat           1 7.291e+11 7.291e+11 355472.6 <2e-16
## color           6 1.256e+10 2.093e+09   1020.6 <2e-16
## cut             4 6.171e+09 1.543e+09    752.1 <2e-16
## Residuals   53928 1.106e+11 2.051e+06
This site uses Just the Docs, a documentation theme for Jekyll.