Formulas in R
Peter Ralph
Advanced Biological Statistics
Set-up
Example data
Continuing with the “pumpkin data”:
## fertilizer water plot plant weight mean_weight
## 1 low no water 1 1 5.289819 5.0
## 2 medium no water 1 1 4.847977 5.5
## 3 high no water 1 1 7.369901 6.0
## 4 low water 1 1 8.690082 6.0
## 5 medium water 1 1 7.994908 6.5
## 6 high water 1 1 7.661236 6.0
## 7 low no water 2 1 4.582871 5.0
## 8 medium no water 2 1 7.028808 5.5
## 9 high no water 2 1 6.068156 6.0
## 10 low water 2 1 8.363669 6.0
## 11 medium water 2 1 6.827433 6.5
## 12 high water 2 1 4.430197 6.0
## 13 low no water 3 1 4.566967 5.0
## 14 medium no water 3 1 6.222325 5.5
## 15 high no water 3 1 4.175435 6.0
## 16 low water 3 1 6.462205 6.0
## 17 medium water 3 1 5.719733 6.5
## 18 high water 3 1 5.336529 6.0
## 19 low no water 4 1 4.101031 5.0
## 20 medium no water 4 1 7.441867 5.5
## 21 high no water 4 1 7.253087 6.0
## 22 low water 4 1 5.150977 6.0
## 23 medium water 4 1 7.682629 6.5
## 24 high water 4 1 4.251106 6.0
## 25 low no water 1 2 4.203770 5.0
## 26 medium no water 1 2 5.498412 5.5
## 27 high no water 1 2 5.231156 6.0
## 28 low water 1 2 7.864301 6.0
## 29 medium water 1 2 5.509126 6.5
## 30 high water 1 2 6.249936 6.0
## 31 low no water 2 2 4.170192 5.0
## 32 medium no water 2 2 5.647272 5.5
## 33 high no water 2 2 4.280001 6.0
## 34 low water 2 2 4.125109 6.0
## 35 medium water 2 2 6.196832 6.5
## 36 high water 2 2 6.266250 6.0
## 37 low no water 3 2 6.961471 5.0
## 38 medium no water 3 2 3.485597 5.5
## 39 high no water 3 2 7.291625 6.0
## 40 low water 3 2 6.349313 6.0
## 41 medium water 3 2 8.063386 6.5
## 42 high water 3 2 7.993647 6.0
## 43 low no water 4 2 6.418937 5.0
## 44 medium no water 4 2 4.021562 5.5
## 45 high no water 4 2 4.176200 6.0
## 46 low water 4 2 3.556191 6.0
## 47 medium water 4 2 6.638474 6.5
## 48 high water 4 2 5.966824 6.0
## 49 low no water 1 3 5.770217 5.0
## 50 medium no water 1 3 5.155670 5.5
## 51 high no water 1 3 7.892274 6.0
## 52 low water 1 3 5.642347 6.0
## 53 medium water 1 3 3.503281 6.5
## 54 high water 1 3 7.051614 6.0
## 55 low no water 2 3 6.831229 5.0
## 56 medium no water 2 3 4.223614 5.5
## 57 high no water 2 3 5.154588 6.0
## 58 low water 2 3 6.415852 6.0
## 59 medium water 2 3 7.956078 6.5
## 60 high water 2 3 7.428946 6.0
## 61 low no water 3 3 4.404570 5.0
## 62 medium no water 3 3 5.496707 5.5
## 63 high no water 3 3 5.449613 6.0
## 64 low water 3 3 6.059843 6.0
## 65 medium water 3 3 9.246146 6.5
## 66 high water 3 3 4.342042 6.0
## 67 low no water 4 3 3.544281 5.0
## 68 medium no water 4 3 3.500048 5.5
## 69 high no water 4 3 7.120510 6.0
## 70 low water 4 3 4.179932 6.0
## 71 medium water 4 3 7.589011 6.5
## 72 high water 4 3 5.355113 6.0
## 73 low no water 1 4 4.801136 5.0
## 74 medium no water 1 4 7.047244 5.5
## 75 high no water 1 4 6.010457 6.0
## 76 low water 1 4 6.650114 6.0
## 77 medium water 1 4 8.431054 6.5
## 78 high water 1 4 2.372604 6.0
## 79 low no water 2 4 2.941023 5.0
## 80 medium no water 2 4 4.780900 5.5
## 81 high no water 2 4 6.016520 6.0
## 82 low water 2 4 6.943116 6.0
## 83 medium water 2 4 9.464095 6.5
## 84 high water 2 4 6.692337 6.0
## 85 low no water 3 4 5.789024 5.0
## 86 medium no water 3 4 2.424725 5.5
## 87 high no water 3 4 5.173750 6.0
## 88 low water 3 4 5.131088 6.0
## 89 medium water 3 4 7.292216 6.5
## 90 high water 3 4 5.790369 6.0
## 91 low no water 4 4 7.318907 5.0
## 92 medium no water 4 4 7.260605 5.5
## 93 high no water 4 4 6.959911 6.0
## 94 low water 4 4 5.927222 6.0
## 95 medium water 4 4 5.848483 6.5
## 96 high water 4 4 5.576530 6.0
## 97 low no water 1 5 7.009595 5.0
## 98 medium no water 1 5 7.069961 5.5
## 99 high no water 1 5 7.123822 6.0
## 100 low water 1 5 6.866997 6.0
## 101 medium water 1 5 7.048750 6.5
## 102 high water 1 5 7.397133 6.0
## 103 low no water 2 5 2.139843 5.0
## 104 medium no water 2 5 5.577983 5.5
## 105 high no water 2 5 4.727250 6.0
## 106 low water 2 5 5.447749 6.0
## 107 medium water 2 5 7.437041 6.5
## 108 high water 2 5 7.830389 6.0
## 109 low no water 3 5 1.640326 5.0
## 110 medium no water 3 5 4.946987 5.5
## 111 high no water 3 5 5.869372 6.0
## 112 low water 3 5 8.012957 6.0
## 113 medium water 3 5 9.041120 6.5
## 114 high water 3 5 6.974253 6.0
## 115 low no water 4 5 5.549274 5.0
## 116 medium no water 4 5 5.037444 5.5
## 117 high no water 4 5 6.771851 6.0
## 118 low water 4 5 4.648319 6.0
## 119 medium water 4 5 6.793488 6.5
## 120 high water 4 5 5.763991 6.0Formulas
Anatomy of a formula
weight ~ fertilizer + wateris read something like
mean weight is determined by additive effects of fertlizer and of water
\[\begin{equation}\begin{split} y &\sim x \qquad \text{means} \\ y &= a + b x + \text{(mean-zero noise)} . \end{split}\end{equation}\]
Intercepts
The intercept is included implicitly, so these are equivalent:
weight ~ fertilizer + water
weight ~ 1 + fertilizer + water… so if you don’t want an intercept, do
weight ~ 0 + fertilizer + waterInteractions
To assign an effect to each element of a crossed design, use
:, e.g.
weight ~ fertilizer + water + fertilizer:waterwhich is the same as
weight ~ fertilizer * watersince lower-order effects are included implicitly.
A translation table
~: depends on+: and also, independently:: in combination with*: and alsoI(x+y): actuallyxplusyI(x^2): actuallyxsquared
Trickier things:
1: an intercept0: but not an intercept-: but not.: all columns not otherwise in the formulax/y:x, andynested withinx(same asx + x:y)
The secret to formulas
If you want to know what a formula is really doing, look at
its model.matrix( ), whose columns correspond to the
coefficients of the resulting model.
## (Intercept) fertilizerlow fertilizermedium
## 1 1 1 0
## 2 1 0 1
## 3 1 0 0
## 4 1 1 0
## 5 1 0 1
## 6 1 0 0
## 7 1 1 0
## 8 1 0 1
## 9 1 0 0
## 10 1 1 0
## 11 1 0 1
## 12 1 0 0
## 13 1 1 0
## 14 1 0 1
## 15 1 0 0
## 16 1 1 0
## 17 1 0 1
## 18 1 0 0
## 19 1 1 0
## 20 1 0 1
## 21 1 0 0
## 22 1 1 0
## 23 1 0 1
## 24 1 0 0
## 25 1 1 0
## 26 1 0 1
## 27 1 0 0
## 28 1 1 0
## 29 1 0 1
## 30 1 0 0
## 31 1 1 0
## 32 1 0 1
## 33 1 0 0
## 34 1 1 0
## 35 1 0 1
## 36 1 0 0
## 37 1 1 0
## 38 1 0 1
## 39 1 0 0
## 40 1 1 0
## 41 1 0 1
## 42 1 0 0
## 43 1 1 0
## 44 1 0 1
## 45 1 0 0
## 46 1 1 0
## 47 1 0 1
## 48 1 0 0
## 49 1 1 0
## 50 1 0 1
## 51 1 0 0
## 52 1 1 0
## 53 1 0 1
## 54 1 0 0
## 55 1 1 0
## 56 1 0 1
## 57 1 0 0
## 58 1 1 0
## 59 1 0 1
## 60 1 0 0
## 61 1 1 0
## 62 1 0 1
## 63 1 0 0
## 64 1 1 0
## 65 1 0 1
## 66 1 0 0
## 67 1 1 0
## 68 1 0 1
## 69 1 0 0
## 70 1 1 0
## 71 1 0 1
## 72 1 0 0
## 73 1 1 0
## 74 1 0 1
## 75 1 0 0
## 76 1 1 0
## 77 1 0 1
## 78 1 0 0
## 79 1 1 0
## 80 1 0 1
## 81 1 0 0
## 82 1 1 0
## 83 1 0 1
## 84 1 0 0
## 85 1 1 0
## 86 1 0 1
## 87 1 0 0
## 88 1 1 0
## 89 1 0 1
## 90 1 0 0
## 91 1 1 0
## 92 1 0 1
## 93 1 0 0
## 94 1 1 0
## 95 1 0 1
## 96 1 0 0
## 97 1 1 0
## 98 1 0 1
## 99 1 0 0
## 100 1 1 0
## 101 1 0 1
## 102 1 0 0
## 103 1 1 0
## 104 1 0 1
## 105 1 0 0
## 106 1 1 0
## 107 1 0 1
## 108 1 0 0
## 109 1 1 0
## 110 1 0 1
## 111 1 0 0
## 112 1 1 0
## 113 1 0 1
## 114 1 0 0
## 115 1 1 0
## 116 1 0 1
## 117 1 0 0
## 118 1 1 0
## 119 1 0 1
## 120 1 0 0
## attr(,"assign")
## [1] 0 1 1
## attr(,"contrasts")
## attr(,"contrasts")$fertilizer
## [1] "contr.treatment"
## (Intercept) fertilizerlow fertilizermedium waterwater
## 1 1 1 0 0
## 2 1 0 1 0
## 3 1 0 0 0
## 4 1 1 0 1
## 5 1 0 1 1
## 6 1 0 0 1
## 7 1 1 0 0
## 8 1 0 1 0
## 9 1 0 0 0
## 10 1 1 0 1
## 11 1 0 1 1
## 12 1 0 0 1
## 13 1 1 0 0
## 14 1 0 1 0
## 15 1 0 0 0
## 16 1 1 0 1
## 17 1 0 1 1
## 18 1 0 0 1
## 19 1 1 0 0
## 20 1 0 1 0
## 21 1 0 0 0
## 22 1 1 0 1
## 23 1 0 1 1
## 24 1 0 0 1
## 25 1 1 0 0
## 26 1 0 1 0
## 27 1 0 0 0
## 28 1 1 0 1
## 29 1 0 1 1
## 30 1 0 0 1
## 31 1 1 0 0
## 32 1 0 1 0
## 33 1 0 0 0
## 34 1 1 0 1
## 35 1 0 1 1
## 36 1 0 0 1
## 37 1 1 0 0
## 38 1 0 1 0
## 39 1 0 0 0
## 40 1 1 0 1
## 41 1 0 1 1
## 42 1 0 0 1
## 43 1 1 0 0
## 44 1 0 1 0
## 45 1 0 0 0
## 46 1 1 0 1
## 47 1 0 1 1
## 48 1 0 0 1
## 49 1 1 0 0
## 50 1 0 1 0
## 51 1 0 0 0
## 52 1 1 0 1
## 53 1 0 1 1
## 54 1 0 0 1
## 55 1 1 0 0
## 56 1 0 1 0
## 57 1 0 0 0
## 58 1 1 0 1
## 59 1 0 1 1
## 60 1 0 0 1
## 61 1 1 0 0
## 62 1 0 1 0
## 63 1 0 0 0
## 64 1 1 0 1
## 65 1 0 1 1
## 66 1 0 0 1
## 67 1 1 0 0
## 68 1 0 1 0
## 69 1 0 0 0
## 70 1 1 0 1
## 71 1 0 1 1
## 72 1 0 0 1
## 73 1 1 0 0
## 74 1 0 1 0
## 75 1 0 0 0
## 76 1 1 0 1
## 77 1 0 1 1
## 78 1 0 0 1
## 79 1 1 0 0
## 80 1 0 1 0
## 81 1 0 0 0
## 82 1 1 0 1
## 83 1 0 1 1
## 84 1 0 0 1
## 85 1 1 0 0
## 86 1 0 1 0
## 87 1 0 0 0
## 88 1 1 0 1
## 89 1 0 1 1
## 90 1 0 0 1
## 91 1 1 0 0
## 92 1 0 1 0
## 93 1 0 0 0
## 94 1 1 0 1
## 95 1 0 1 1
## 96 1 0 0 1
## 97 1 1 0 0
## 98 1 0 1 0
## 99 1 0 0 0
## 100 1 1 0 1
## 101 1 0 1 1
## 102 1 0 0 1
## 103 1 1 0 0
## 104 1 0 1 0
## 105 1 0 0 0
## 106 1 1 0 1
## 107 1 0 1 1
## 108 1 0 0 1
## 109 1 1 0 0
## 110 1 0 1 0
## 111 1 0 0 0
## 112 1 1 0 1
## 113 1 0 1 1
## 114 1 0 0 1
## 115 1 1 0 0
## 116 1 0 1 0
## 117 1 0 0 0
## 118 1 1 0 1
## 119 1 0 1 1
## 120 1 0 0 1
## attr(,"assign")
## [1] 0 1 1 2
## attr(,"contrasts")
## attr(,"contrasts")$fertilizer
## [1] "contr.treatment"
##
## attr(,"contrasts")$water
## [1] "contr.treatment"
## fertilizerhigh fertilizerlow fertilizermedium waterwater fertilizerlow:waterwater fertilizermedium:waterwater
## 1 0 1 0 0 0 0
## 2 0 0 1 0 0 0
## 3 1 0 0 0 0 0
## 4 0 1 0 1 1 0
## 5 0 0 1 1 0 1
## 6 1 0 0 1 0 0
## 7 0 1 0 0 0 0
## 8 0 0 1 0 0 0
## 9 1 0 0 0 0 0
## 10 0 1 0 1 1 0
## 11 0 0 1 1 0 1
## 12 1 0 0 1 0 0
## 13 0 1 0 0 0 0
## 14 0 0 1 0 0 0
## 15 1 0 0 0 0 0
## 16 0 1 0 1 1 0
## 17 0 0 1 1 0 1
## 18 1 0 0 1 0 0
## 19 0 1 0 0 0 0
## 20 0 0 1 0 0 0
## 21 1 0 0 0 0 0
## 22 0 1 0 1 1 0
## 23 0 0 1 1 0 1
## 24 1 0 0 1 0 0
## 25 0 1 0 0 0 0
## 26 0 0 1 0 0 0
## 27 1 0 0 0 0 0
## 28 0 1 0 1 1 0
## 29 0 0 1 1 0 1
## 30 1 0 0 1 0 0
## 31 0 1 0 0 0 0
## 32 0 0 1 0 0 0
## 33 1 0 0 0 0 0
## 34 0 1 0 1 1 0
## 35 0 0 1 1 0 1
## 36 1 0 0 1 0 0
## 37 0 1 0 0 0 0
## 38 0 0 1 0 0 0
## 39 1 0 0 0 0 0
## 40 0 1 0 1 1 0
## 41 0 0 1 1 0 1
## 42 1 0 0 1 0 0
## 43 0 1 0 0 0 0
## 44 0 0 1 0 0 0
## 45 1 0 0 0 0 0
## 46 0 1 0 1 1 0
## 47 0 0 1 1 0 1
## 48 1 0 0 1 0 0
## 49 0 1 0 0 0 0
## 50 0 0 1 0 0 0
## 51 1 0 0 0 0 0
## 52 0 1 0 1 1 0
## 53 0 0 1 1 0 1
## 54 1 0 0 1 0 0
## 55 0 1 0 0 0 0
## 56 0 0 1 0 0 0
## 57 1 0 0 0 0 0
## 58 0 1 0 1 1 0
## 59 0 0 1 1 0 1
## 60 1 0 0 1 0 0
## 61 0 1 0 0 0 0
## 62 0 0 1 0 0 0
## 63 1 0 0 0 0 0
## 64 0 1 0 1 1 0
## 65 0 0 1 1 0 1
## 66 1 0 0 1 0 0
## 67 0 1 0 0 0 0
## 68 0 0 1 0 0 0
## 69 1 0 0 0 0 0
## 70 0 1 0 1 1 0
## 71 0 0 1 1 0 1
## 72 1 0 0 1 0 0
## 73 0 1 0 0 0 0
## 74 0 0 1 0 0 0
## 75 1 0 0 0 0 0
## 76 0 1 0 1 1 0
## 77 0 0 1 1 0 1
## 78 1 0 0 1 0 0
## 79 0 1 0 0 0 0
## 80 0 0 1 0 0 0
## 81 1 0 0 0 0 0
## 82 0 1 0 1 1 0
## 83 0 0 1 1 0 1
## 84 1 0 0 1 0 0
## 85 0 1 0 0 0 0
## 86 0 0 1 0 0 0
## 87 1 0 0 0 0 0
## 88 0 1 0 1 1 0
## 89 0 0 1 1 0 1
## 90 1 0 0 1 0 0
## 91 0 1 0 0 0 0
## 92 0 0 1 0 0 0
## 93 1 0 0 0 0 0
## 94 0 1 0 1 1 0
## 95 0 0 1 1 0 1
## 96 1 0 0 1 0 0
## 97 0 1 0 0 0 0
## 98 0 0 1 0 0 0
## 99 1 0 0 0 0 0
## 100 0 1 0 1 1 0
## 101 0 0 1 1 0 1
## 102 1 0 0 1 0 0
## 103 0 1 0 0 0 0
## 104 0 0 1 0 0 0
## 105 1 0 0 0 0 0
## 106 0 1 0 1 1 0
## 107 0 0 1 1 0 1
## 108 1 0 0 1 0 0
## 109 0 1 0 0 0 0
## 110 0 0 1 0 0 0
## 111 1 0 0 0 0 0
## 112 0 1 0 1 1 0
## 113 0 0 1 1 0 1
## 114 1 0 0 1 0 0
## 115 0 1 0 0 0 0
## 116 0 0 1 0 0 0
## 117 1 0 0 0 0 0
## 118 0 1 0 1 1 0
## 119 0 0 1 1 0 1
## 120 1 0 0 1 0 0
## attr(,"assign")
## [1] 1 1 1 2 3 3
## attr(,"contrasts")
## attr(,"contrasts")$fertilizer
## [1] "contr.treatment"
##
## attr(,"contrasts")$water
## [1] "contr.treatment"
Note:
For fine control of factors in linear models, either
relevel()them, or- manually set their
contrasts().
Exercise:
Make formulas that give you estimates of
A global mean (\(\mu\)), two fertilizer effects (\(\alpha_\text{medium}\) and \(\alpha_\text{high}\)), and one water effect (\(\beta_\text{water}\)).
Three fertilizer effects (\(\alpha_\text{low}\), \(\alpha_\text{medium}\) and \(\alpha_\text{high}\)), and one water effect (\(\beta_\text{water}\)).
Two fertilizer effects (\(\alpha_\text{medium}\) and \(\alpha_\text{high}\)), and two water effect (\(\beta_\text{no water}\) and \(\beta_\text{water}\)).
A single mean for each of the six conditions (\(\gamma_\text{high, water}\), \(\gamma_\text{medium, water}\), \(\gamma_\text{low, water}\), \(\gamma_\text{high, no water}\), \(\gamma_\text{medium, no water}\), \(\gamma_\text{low, no water}\)).
Data: pumkpins.tsv
Example:
- A global mean (\(\mu\)), two fertilizer effects (\(\alpha_\text{medium}\) and \(\alpha_\text{high}\)), and one water effect (\(\gamma_\text{water}\)).
Example:
##
## Call:
## lm(formula = weight ~ fertilizer + water, data = pumpkins)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.1705 -0.8935 -0.0265 1.1688 2.6672
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.4992 0.2663 20.649 < 2e-16 ***
## fertilizerlow -0.5081 0.3262 -1.558 0.12199
## fertilizermedium 0.2538 0.3262 0.778 0.43805
## waterwater 1.0439 0.2663 3.920 0.00015 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.459 on 116 degrees of freedom
## Multiple R-squared: 0.1534, Adjusted R-squared: 0.1315
## F-statistic: 7.009 on 3 and 116 DF, p-value: 0.0002254