Formulas in R
Peter Ralph
20 October – 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 2.0966062 2
## 2 medium no water 1 1 6.7826589 7
## 3 high no water 1 1 12.4566335 12
## 4 low water 1 1 8.8966940 8
## 5 medium water 1 1 10.4983026 10
## 6 high water 1 1 12.5537452 12
## 7 low no water 2 1 1.8609569 2
## 8 medium no water 2 1 7.5096027 7
## 9 high no water 2 1 12.0227186 12
## 10 low water 2 1 8.7878898 8
## 11 medium water 2 1 10.1091442 10
## 12 high water 2 1 11.4767323 12
## 13 low no water 3 1 1.8556557 2
## 14 medium no water 3 1 7.2407751 7
## 15 high no water 3 1 11.3918118 12
## 16 low water 3 1 8.1540684 8
## 17 medium water 3 1 9.7399108 10
## 18 high water 3 1 11.7788431 12
## 19 low no water 4 1 1.7003436 2
## 20 medium no water 4 1 7.6472889 7
## 21 high no water 4 1 12.4176956 12
## 22 low water 4 1 7.7169925 8
## 23 medium water 4 1 10.3942097 10
## 24 high water 4 1 11.4170353 12
## 25 low no water 1 2 1.7345900 2
## 26 medium no water 1 2 6.9994706 7
## 27 high no water 1 2 11.7437188 12
## 28 low water 1 2 8.6214338 8
## 29 medium water 1 2 9.6697086 10
## 30 high water 1 2 12.0833121 12
## 31 low no water 2 2 1.7233974 2
## 32 medium no water 2 2 7.0490907 7
## 33 high no water 2 2 11.4266671 12
## 34 low water 2 2 7.3750364 8
## 35 medium water 2 2 9.8989442 10
## 36 high water 2 2 12.0887502 12
## 37 low no water 3 2 2.6538237 2
## 38 medium no water 3 2 6.3285323 7
## 39 high no water 3 2 12.4305415 12
## 40 low water 3 2 8.1164375 8
## 41 medium water 3 2 10.5211285 10
## 42 high water 3 2 12.6645491 12
## 43 low no water 4 2 2.4729790 2
## 44 medium no water 4 2 6.5071872 7
## 45 high no water 4 2 11.3920666 12
## 46 low water 4 2 7.1853970 8
## 47 medium water 4 2 10.0461581 10
## 48 high water 4 2 11.9889413 12
## 49 low no water 1 3 2.2567390 2
## 50 medium no water 1 3 6.8852232 7
## 51 high no water 1 3 12.6307580 12
## 52 low water 1 3 7.8807822 8
## 53 medium water 1 3 9.0010937 10
## 54 high water 1 3 12.3505380 12
## 55 low no water 2 3 2.6104095 2
## 56 medium no water 2 3 6.5745380 7
## 57 high no water 2 3 11.7181961 12
## 58 low water 2 3 8.1386172 8
## 59 medium water 2 3 10.4853594 10
## 60 high water 2 3 12.4763154 12
## 61 low no water 3 3 1.8015233 2
## 62 medium no water 3 3 6.9989023 7
## 63 high no water 3 3 11.8165377 12
## 64 low water 3 3 8.0199476 8
## 65 medium water 3 3 10.9153819 10
## 66 high water 3 3 11.4473475 12
## 67 low no water 4 3 1.5147604 2
## 68 medium no water 4 3 6.3333493 7
## 69 high no water 4 3 12.3735033 12
## 70 low water 4 3 7.3933107 8
## 71 medium water 4 3 10.3630038 10
## 72 high water 4 3 11.7850376 12
## 73 low no water 1 4 1.9337121 2
## 74 medium no water 1 4 7.5157479 7
## 75 high no water 1 4 12.0034857 12
## 76 low water 1 4 8.2167045 8
## 77 medium water 1 4 10.6436847 10
## 78 high water 1 4 10.7908679 12
## 79 low no water 2 4 1.3136743 2
## 80 medium no water 2 4 6.7603000 7
## 81 high no water 2 4 12.0055068 12
## 82 low water 2 4 8.3143720 8
## 83 medium water 2 4 10.9880317 10
## 84 high water 2 4 12.2307791 12
## 85 low no water 3 4 2.2630081 2
## 86 medium no water 3 4 5.9749084 7
## 87 high no water 3 4 11.7245832 12
## 88 low water 3 4 7.7103625 8
## 89 medium water 3 4 10.2640720 10
## 90 high water 3 4 11.9301231 12
## 91 low no water 4 4 2.7729689 2
## 92 medium no water 4 4 7.5868684 7
## 93 high no water 4 4 12.3199703 12
## 94 low water 4 4 7.9757408 8
## 95 medium water 4 4 9.7828276 10
## 96 high water 4 4 11.8588434 12
## 97 low no water 1 5 2.6698649 2
## 98 medium no water 1 5 7.5233203 7
## 99 high no water 1 5 12.3746075 12
## 100 low water 1 5 8.2889990 8
## 101 medium water 1 5 10.1829166 10
## 102 high water 1 5 12.4657109 12
## 103 low no water 2 5 1.0466142 2
## 104 medium no water 2 5 7.0259943 7
## 105 high no water 2 5 11.5757501 12
## 106 low water 2 5 7.8159163 8
## 107 medium water 2 5 10.3123470 10
## 108 high water 2 5 12.6101296 12
## 109 low no water 3 5 0.8801086 2
## 110 medium no water 3 5 6.8156623 7
## 111 high no water 3 5 11.9564572 12
## 112 low water 3 5 8.6709858 8
## 113 medium water 3 5 10.8470400 10
## 114 high water 3 5 12.3247510 12
## 115 low no water 4 5 2.1830913 2
## 116 medium no water 4 5 6.8458148 7
## 117 high no water 4 5 12.2572836 12
## 118 low water 4 5 7.5494395 8
## 119 medium water 4 5 10.0978294 10
## 120 high water 4 5 11.9213303 12Formulas
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}\)).
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
## -2.77906 -1.15970 -0.01605 1.17450 2.32947
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.4442 0.2455 42.54 <2e-16 ***
## fertilizerlow -7.0027 0.3007 -23.29 <2e-16 ***
## fertilizermedium -3.4154 0.3007 -11.36 <2e-16 ***
## waterwater 3.1258 0.2455 12.73 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.345 on 116 degrees of freedom
## Multiple R-squared: 0.8586, Adjusted R-squared: 0.855
## F-statistic: 234.9 on 3 and 116 DF, p-value: < 2.2e-16