Formulas in R
Peter Ralph
Advanced Biological Statistics
Set-up
Example data
Continuing with the “pumpkin data”:
## fertilizer water plot plant weight mean_weight interaction
## 1 low no water 1 1 1.0966062 1 low,no water
## 2 medium no water 1 1 1.7826589 2 medium,no water
## 3 high no water 1 1 3.4566335 3 high,no water
## 4 low water 1 1 4.8966940 4 low,water
## 5 medium water 1 1 5.4983026 5 medium,water
## 6 high water 1 1 9.5537452 9 high,water
## 7 low no water 2 1 0.8609569 1 low,no water
## 8 medium no water 2 1 2.5096027 2 medium,no water
## 9 high no water 2 1 3.0227186 3 high,no water
## 10 low water 2 1 4.7878898 4 low,water
## 11 medium water 2 1 5.1091442 5 medium,water
## 12 high water 2 1 8.4767323 9 high,water
## 13 low no water 3 1 0.8556557 1 low,no water
## 14 medium no water 3 1 2.2407751 2 medium,no water
## 15 high no water 3 1 2.3918118 3 high,no water
## 16 low water 3 1 4.1540684 4 low,water
## 17 medium water 3 1 4.7399108 5 medium,water
## 18 high water 3 1 8.7788431 9 high,water
## 19 low no water 4 1 0.7003436 1 low,no water
## 20 medium no water 4 1 2.6472889 2 medium,no water
## 21 high no water 4 1 3.4176956 3 high,no water
## 22 low water 4 1 3.7169925 4 low,water
## 23 medium water 4 1 5.3942097 5 medium,water
## 24 high water 4 1 8.4170353 9 high,water
## 25 low no water 1 2 0.7345900 1 low,no water
## 26 medium no water 1 2 1.9994706 2 medium,no water
## 27 high no water 1 2 2.7437188 3 high,no water
## 28 low water 1 2 4.6214338 4 low,water
## 29 medium water 1 2 4.6697086 5 medium,water
## 30 high water 1 2 9.0833121 9 high,water
## 31 low no water 2 2 0.7233974 1 low,no water
## 32 medium no water 2 2 2.0490907 2 medium,no water
## 33 high no water 2 2 2.4266671 3 high,no water
## 34 low water 2 2 3.3750364 4 low,water
## 35 medium water 2 2 4.8989442 5 medium,water
## 36 high water 2 2 9.0887502 9 high,water
## 37 low no water 3 2 1.6538237 1 low,no water
## 38 medium no water 3 2 1.3285323 2 medium,no water
## 39 high no water 3 2 3.4305415 3 high,no water
## 40 low water 3 2 4.1164375 4 low,water
## 41 medium water 3 2 5.5211285 5 medium,water
## 42 high water 3 2 9.6645491 9 high,water
## 43 low no water 4 2 1.4729790 1 low,no water
## 44 medium no water 4 2 1.5071872 2 medium,no water
## 45 high no water 4 2 2.3920666 3 high,no water
## 46 low water 4 2 3.1853970 4 low,water
## 47 medium water 4 2 5.0461581 5 medium,water
## 48 high water 4 2 8.9889413 9 high,water
## 49 low no water 1 3 1.2567390 1 low,no water
## 50 medium no water 1 3 1.8852232 2 medium,no water
## 51 high no water 1 3 3.6307580 3 high,no water
## 52 low water 1 3 3.8807822 4 low,water
## 53 medium water 1 3 4.0010937 5 medium,water
## 54 high water 1 3 9.3505380 9 high,water
## 55 low no water 2 3 1.6104095 1 low,no water
## 56 medium no water 2 3 1.5745380 2 medium,no water
## 57 high no water 2 3 2.7181961 3 high,no water
## 58 low water 2 3 4.1386172 4 low,water
## 59 medium water 2 3 5.4853594 5 medium,water
## 60 high water 2 3 9.4763154 9 high,water
## 61 low no water 3 3 0.8015233 1 low,no water
## 62 medium no water 3 3 1.9989023 2 medium,no water
## 63 high no water 3 3 2.8165377 3 high,no water
## 64 low water 3 3 4.0199476 4 low,water
## 65 medium water 3 3 5.9153819 5 medium,water
## 66 high water 3 3 8.4473475 9 high,water
## 67 low no water 4 3 0.5147604 1 low,no water
## 68 medium no water 4 3 1.3333493 2 medium,no water
## 69 high no water 4 3 3.3735033 3 high,no water
## 70 low water 4 3 3.3933107 4 low,water
## 71 medium water 4 3 5.3630038 5 medium,water
## 72 high water 4 3 8.7850376 9 high,water
## 73 low no water 1 4 0.9337121 1 low,no water
## 74 medium no water 1 4 2.5157479 2 medium,no water
## 75 high no water 1 4 3.0034857 3 high,no water
## 76 low water 1 4 4.2167045 4 low,water
## 77 medium water 1 4 5.6436847 5 medium,water
## 78 high water 1 4 7.7908679 9 high,water
## 79 low no water 2 4 0.3136743 1 low,no water
## 80 medium no water 2 4 1.7603000 2 medium,no water
## 81 high no water 2 4 3.0055068 3 high,no water
## 82 low water 2 4 4.3143720 4 low,water
## 83 medium water 2 4 5.9880317 5 medium,water
## 84 high water 2 4 9.2307791 9 high,water
## 85 low no water 3 4 1.2630081 1 low,no water
## 86 medium no water 3 4 0.9749084 2 medium,no water
## 87 high no water 3 4 2.7245832 3 high,no water
## 88 low water 3 4 3.7103625 4 low,water
## 89 medium water 3 4 5.2640720 5 medium,water
## 90 high water 3 4 8.9301231 9 high,water
## 91 low no water 4 4 1.7729689 1 low,no water
## 92 medium no water 4 4 2.5868684 2 medium,no water
## 93 high no water 4 4 3.3199703 3 high,no water
## 94 low water 4 4 3.9757408 4 low,water
## 95 medium water 4 4 4.7828276 5 medium,water
## 96 high water 4 4 8.8588434 9 high,water
## 97 low no water 1 5 1.6698649 1 low,no water
## 98 medium no water 1 5 2.5233203 2 medium,no water
## 99 high no water 1 5 3.3746075 3 high,no water
## 100 low water 1 5 4.2889990 4 low,water
## 101 medium water 1 5 5.1829166 5 medium,water
## 102 high water 1 5 9.4657109 9 high,water
## 103 low no water 2 5 0.0466142 1 low,no water
## 104 medium no water 2 5 2.0259943 2 medium,no water
## 105 high no water 2 5 2.5757501 3 high,no water
## 106 low water 2 5 3.8159163 4 low,water
## 107 medium water 2 5 5.3123470 5 medium,water
## 108 high water 2 5 9.6101296 9 high,water
## 109 low no water 3 5 -0.1198914 1 low,no water
## 110 medium no water 3 5 1.8156623 2 medium,no water
## 111 high no water 3 5 2.9564572 3 high,no water
## 112 low water 3 5 4.6709858 4 low,water
## 113 medium water 3 5 5.8470400 5 medium,water
## 114 high water 3 5 9.3247510 9 high,water
## 115 low no water 4 5 1.1830913 1 low,no water
## 116 medium no water 4 5 1.8458148 2 medium,no water
## 117 high no water 4 5 3.2572836 3 high,no water
## 118 low water 4 5 3.5494395 4 low,water
## 119 medium water 4 5 5.0978294 5 medium,water
## 120 high water 4 5 8.9213303 9 high,waterFormulas
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
## -1.65344 -0.56985 -0.01622 0.66858 1.59462
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.9442 0.1507 26.18 <2e-16 ***
## fertilizerlow -3.5027 0.1845 -18.98 <2e-16 ***
## fertilizermedium -2.4154 0.1845 -13.09 <2e-16 ***
## waterwater 4.1258 0.1507 27.39 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8252 on 116 degrees of freedom
## Multiple R-squared: 0.9067, Adjusted R-squared: 0.9043
## F-statistic: 375.9 on 3 and 116 DF, p-value: < 2.2e-16