Workshop 11 - The power of Contrasts

23 April 2011

Topics

Contrasts in regression
Contrasts in LMM
Contrasts in GLMM
- Contrasts
- manual calculations
Contrasts in ANOVA
Contrasts in factorial ANOVA

Context

Recent statistical advice has progressively advocated the use of effect sizes and confidence intervals over the traditional p-values of hyothesis tests.
It is argued that the size of p-values are regularly incorrectly inferred to represent the strength of an effect.
Proponents argue that effect sizes (the magnitude of the differences between treatment levels or trajectories) and associated measures of confidence or precision of these estimates provide far more informative summaries.
Furthermore, the resolve of these recommendations is greatly enhanced as model complexity increases (particularly with the introduction of nesting or repeated measures) and appropriate degrees of freedom on which to base p-value calculations become increasingly more tenuous.
Consequently, there have been numerous pieces of work that oppose/criticize/impugn the use of p-values and advocate the use of effect sizes.
Unfortunately, few offer much guidance in how such effects sizes can be calculated, presented and interpreted.
Therefore, the reader is left questioning the way that they perform and present their statistical findings and is confused about how they should proceed in the future knowing that if they continue to present their results in the traditional framework, they are likely to be challenged during the review process.
This paper aims to address these issues by illustrating the following:
- The use of contrast coefficients to calculate effect sizes, marginal means and associated confidence intervals for a range of common statistical designs.
- The techniques apply equally to models of any complexity and model fitting procedure (including Bayesian models).
- Graphical presentation of effect sizes
- All examples are backed up with full R code (in appendix)

All of the following data sets are intentionally artificial. By creating data sets with specific characteristics (samples drawn from exactly known populations), we can better assess the accuracy of the techniques.

Contrasts in prediction in regression analyses

Imagine a situation in which we had just performed simple linear regression, and now we wished to summarize the trend graphically. For this demonstration, we will return to a dataset we used in a previous Workshop (the fertilizer data set): An agriculturalist was interested in the effects of fertilizer load on the yield of grass. Grass seed was sown uniformly over an area and different quantities of commercial fertilizer were applied to each of ten 1 m² randomly located plots. Two months later the grass from each plot was harvested, dried and weighed. The data are in the file fertilizer.csv.

Download Fertilizer data set

Format of fertilizer.csv data files

FERTILIZER	YIELD
25	84
50	80
75	90
100	154
125	148
...	...

FERTILIZER	Mass of fertilizer (g.m^-2) - Predictor variable
YIELD	Yield of grass (g.m^-2) - Response variable

Open the fertilizer dataset

Show code

> fert <- read.table("../downloads/data/fertilizer.csv", header = T, 
+     sep = ",", strip.white = T)
> fert

   FERTILIZER YIELD
1          25    84
2          50    80
3          75    90
4         100   154
5         125   148
6         150   169
7         175   206
8         200   244
9         225   212
10        250   248

Linear model parameterization for regression looks like:

y_i = β₀ + β_ix₁ + ε_i
ε_i ∼ Normal(0,σ²)	#Residual 'random' effects

Lets fit the linear model relating grass yield to fertilizer concentration:

Show code

> fert.lm <- lm(YIELD ~ FERTILIZER, data = fert)
> fert.lm

Call:
lm(formula = YIELD ~ FERTILIZER, data = fert)

Coefficients:
(Intercept)   FERTILIZER  
     51.933        0.811

When we fit a linear model, we estimate the values of the parameters (β₀ and β₁) such that the deviations in the predicted values of y_i are minimized (OLS). Once we have the estimates of the intercept and slope parameters, we can solve the equation for new values of x so as to predict new values of y.

For example, if the estimated intercept (β₀) and slope (β₀) where 2 and 3 respectively, we could predict a new value of y when x was equal to 5.

y_i = β₀ + β₁x_i

y_i = 2 + 3× x_i

y = 2 + 3× 5

y = 17

In other words, what we are saying, is that to estimate a new value of y, we want to add one unit of 2 to five units of 3. So the contrasts associated with predicting a new value of y from a new x value of 5 would be 1 and 5

y_i =

β₀

β₁

y_i =

So what we are after is what is called the outer product of the two matrices. The outer product multiplies the first entry of the left matrix by the first value (first row) of the second matrix and the second values together and so on, and then sums them all up. In R, the outer product of two matrices is performed by the %*%.

> # Create a sequence of x-values that span the range of the Fertilizer
> #   variable
> coefs <- c(2, 3)
> cmat <- c(1, 5)
> coefs %*% cmat

     [,1]
[1,]   17

In this way, we can predict multiple new y values from multiple new x values. When the second matrix contains multiple columns, each column represents a separate contrast. Therefore the outer products are performed separately for each column of the second (contrast) matrix.

y_i =

1	1	1
5	5.5	8

> # Create a sequence of x-values that span the range of the Fertilizer
> #   variable
> coefs <- c(2, 3)
> cmat <- cbind(c(1, 5), c(1, 5.5), c(1, 8))
> cmat

     [,1] [,2] [,3]
[1,]    1  1.0    1
[2,]    5  5.5    8

> coefs %*% cmat

     [,1] [,2] [,3]
[1,]   17 18.5   26

Use the predict() function to extract the line of best fit and standard errors

For simple examples such as this, there is a function called predict() that can be used to derive estimates (and standard errors) from simple linear models. This function is essentially a clever wrapper for the series of actions that we will perform manually. As it is well tested under simple conditions, it provides a good comparison and reference point for our manual calculations.

> # Create a sequence of x-values that span the range of the Fertilizer
> #   variable
> newX <- seq(min(fert$FERTILIZER), max(fert$FERTILIZER), l = 100)
> pred <- predict(fert.lm, newdata = data.frame(FERTILIZER = newX), 
+     se = T)

Show the list of predicted fits and standard errors

> pred

$fit
     1      2      3      4      5      6      7      8      9     10     11     12     13     14 
 72.22  74.06  75.91  77.75  79.59  81.44  83.28  85.13  86.97  88.81  90.66  92.50  94.35  96.19 
    15     16     17     18     19     20     21     22     23     24     25     26     27     28 
 98.04  99.88 101.72 103.57 105.41 107.26 109.10 110.94 112.79 114.63 116.48 118.32 120.16 122.01 
    29     30     31     32     33     34     35     36     37     38     39     40     41     42 
123.85 125.70 127.54 129.38 131.23 133.07 134.92 136.76 138.60 140.45 142.29 144.14 145.98 147.83 
    43     44     45     46     47     48     49     50     51     52     53     54     55     56 
149.67 151.51 153.36 155.20 157.05 158.89 160.73 162.58 164.42 166.27 168.11 169.95 171.80 173.64 
    57     58     59     60     61     62     63     64     65     66     67     68     69     70 
175.49 177.33 179.17 181.02 182.86 184.71 186.55 188.40 190.24 192.08 193.93 195.77 197.62 199.46 
    71     72     73     74     75     76     77     78     79     80     81     82     83     84 
201.30 203.15 204.99 206.84 208.68 210.52 212.37 214.21 216.06 217.90 219.74 221.59 223.43 225.28 
    85     86     87     88     89     90     91     92     93     94     95     96     97     98 
227.12 228.96 230.81 232.65 234.50 236.34 238.19 240.03 241.87 243.72 245.56 247.41 249.25 251.09 
    99    100 
252.94 254.78 

$se.fit
     1      2      3      4      5      6      7      8      9     10     11     12     13     14 
11.167 11.007 10.848 10.690 10.534 10.378 10.224 10.070  9.918  9.768  9.619  9.471  9.325  9.180 
    15     16     17     18     19     20     21     22     23     24     25     26     27     28 
 9.037  8.896  8.757  8.619  8.484  8.351  8.220  8.091  7.965  7.842  7.721  7.603  7.488  7.376 
    29     30     31     32     33     34     35     36     37     38     39     40     41     42 
 7.267  7.162  7.060  6.962  6.868  6.778  6.692  6.611  6.534  6.461  6.394  6.331  6.274  6.222 
    43     44     45     46     47     48     49     50     51     52     53     54     55     56 
 6.175  6.134  6.098  6.069  6.045  6.027  6.015  6.009  6.009  6.015  6.027  6.045  6.069  6.098 
    57     58     59     60     61     62     63     64     65     66     67     68     69     70 
 6.134  6.175  6.222  6.274  6.331  6.394  6.461  6.534  6.611  6.692  6.778  6.868  6.962  7.060 
    71     72     73     74     75     76     77     78     79     80     81     82     83     84 
 7.162  7.267  7.376  7.488  7.603  7.721  7.842  7.965  8.091  8.220  8.351  8.484  8.619  8.757 
    85     86     87     88     89     90     91     92     93     94     95     96     97     98 
 8.896  9.037  9.180  9.325  9.471  9.619  9.768  9.918 10.070 10.224 10.378 10.534 10.690 10.848 
    99    100 
11.007 11.167 

$df
[1] 8

$residual.scale
[1] 19

Use these predicted values to draw a line of best fit and standard error (and confidence) bands

> plot(YIELD ~ FERTILIZER, fert, pch = 16)
> points(pred$fit ~ newX, type = "l")
> points(pred$fit - pred$se.fit ~ newX, type = "l", col = "grey")
> points(pred$fit + pred$se.fit ~ newX, type = "l", col = "grey")

> #Now for 95% confidence intervals
> plot(YIELD ~ FERTILIZER, fert, pch = 16)
> points(pred$fit ~ newX, type = "l")
> lwrCI <- pred$fit + qnorm(0.025) * pred$se.fit
> uprCI <- pred$fit + qnorm(0.975) * pred$se.fit
> points(lwrCI ~ newX, type = "l", col = "blue")
> points(uprCI ~ newX, type = "l", col = "blue")

Use contrasts and matrix algebra to derive predictions and errors

First retrieve the fitted model parameters (intercept and slope)

> coefs <- fert.lm$coef
> coefs

(Intercept)  FERTILIZER 
    51.9333      0.8114

Generate a model matrix appropriate for the predictor variable range. The more numbers there are over the range the smoother (less jagged) the fitted curve.
> cc <- as.matrix(expand.grid(1, newX)) > head(cc)
Var1 Var2 [1,] 1 25.00 [2,] 1 27.27 [3,] 1 29.55 [4,] 1 31.82 [5,] 1 34.09 [6,] 1 36.36
Note that the contrast matrix that we just generated has the individual contrasts in rows whereas (recall) they need to be in columns. We can either transpose this matrix (with the t() function) and restore it, or else just use the t() function inline when calculating the outer product of the coefficients and contrast matrices.
Essentially, what we want to do is:

y_i =

51.933 0.8114

×

1.000 1.000 1.000 ...

25.000 27.273 29.545 ...

Calculate (derive) the predicted yield values for each of the new Fertilizer concentrations.

> newY <- coefs %*% t(cc)
> newY

      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] 72.22 74.06 75.91 77.75 79.59 81.44 83.28 85.13 86.97 88.81 90.66  92.5 94.35 96.19 98.04
     [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30]
[1,] 99.88 101.7 103.6 105.4 107.3 109.1 110.9 112.8 114.6 116.5 118.3 120.2   122 123.9 125.7
     [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45]
[1,] 127.5 129.4 131.2 133.1 134.9 136.8 138.6 140.4 142.3 144.1   146 147.8 149.7 151.5 153.4
     [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60]
[1,] 155.2   157 158.9 160.7 162.6 164.4 166.3 168.1   170 171.8 173.6 175.5 177.3 179.2   181
     [,61] [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75]
[1,] 182.9 184.7 186.6 188.4 190.2 192.1 193.9 195.8 197.6 199.5 201.3 203.1   205 206.8 208.7
     [,76] [,77] [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90]
[1,] 210.5 212.4 214.2 216.1 217.9 219.7 221.6 223.4 225.3 227.1   229 230.8 232.7 234.5 236.3
     [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100]
[1,] 238.2   240 241.9 243.7 245.6 247.4 249.2 251.1 252.9  254.8

Calculate (derive) the standard errors of the predicted yield values for each of the new Fertilizer concentrations.

> vc <- vcov(fert.lm)
> se <- sqrt(diag(cc %*% vc %*% t(cc)))
> head(se)

[1] 11.17 11.01 10.85 10.69 10.53 10.38

Reconstruct the summary plot with the derived fit, standard errors and confidence interval

> plot(YIELD ~ FERTILIZER, fert, pch = 16)
> points(newY[1, ] ~ newX, type = "l")
> points(newY[1, ] - se ~ newX, type = "l", col = "grey")
> points(newY[1, ] + se ~ newX, type = "l", col = "grey")
> lwrCI <- newY[1, ] + qnorm(0.025) * se
> uprCI <- newY[1, ] + qnorm(0.975) * se
> points(lwrCI ~ newX, type = "l", col = "blue")
> points(uprCI ~ newX, type = "l", col = "blue")

Contrasts in linear mixed effects models

For the next demonstration, I will use completely fabricated data. Perhaps one day I will seek out a good example with real data, but for now, the artificial data will do.

Create the artificial data

This data set consists of fish counts over time on seven reefs (you can tell it is artificial, because they have increased ;^)

> set.seed(1)
> x <- NULL
> y <- NULL
> for (i in 1:7) {
+     xx <- 2006:2010
+     #yy <- round(3+2*x + rnorm(5,0,1) + rnorm(1,0,3),0)+1
+     
+     yy <- round(-8300 + 4.14 * x + rnorm(5, 0, 1) + rnorm(1, 0, 3), 0) + 1
+     yy[yy < 0] <- 0
+     x <- c(x, xx)
+     y <- c(y, yy)
+ }
> data <- data.frame(x, y, reef = gl(7, 5, 35, LETTERS[1:7]))
> #use a scatterplot to explore the trend for each reef
> library(car)
> scatterplot(y ~ x | reef, data)

Lets fit the generalized linear mixed effects model relating fish counts (y) to year (x). To do so, we will employ the glmmPQL() function from within the MASS package. The glmmPQL fits generalized linear mixed effects models with penalized quasi-likelihood estimation. We will initially use a Gaussian (Normal) error distribution.

> library(MASS)
> data.lme <- glmmPQL(y ~ x, random = ~1 | reef, data, family = "gaussian")
> summary(data.lme)

Linear mixed-effects model fit by maximum likelihood
 Data: data 
  AIC BIC logLik
   NA  NA     NA

Random effects:
 Formula: ~1 | reef
        (Intercept) Residual
StdDev:       1.362    1.744

Variance function:
 Structure: fixed weights
 Formula: ~invwt 
Fixed effects: y ~ x 
            Value Std.Error DF t-value p-value
(Intercept) -7875    244.02 97  -32.27       0
x               4      0.12 97   32.33       0
 Correlation: 
  (Intr)
x -1    

Standardized Within-Group Residuals:
    Min      Q1     Med      Q3     Max 
-1.7777 -0.5180  0.1723  0.6385  1.8105 

Number of Observations: 105
Number of Groups: 7

Use the predict() function to extract the line of best fit and standard errors

> # Create a sequence of x-values that span the range of the x variable
> newX <- seq(min(data$x), max(data$x), l = 100)
> pred <- predict(data.lme, newdata = data.frame(x = newX), level = 0, 
+     se = T)
> pred

  [1]  5.543  5.702  5.860  6.019  6.178  6.337  6.495  6.654  6.813  6.971  7.130  7.289  7.448
 [14]  7.606  7.765  7.924  8.083  8.241  8.400  8.559  8.717  8.876  9.035  9.194  9.352  9.511
 [27]  9.670  9.829  9.987 10.146 10.305 10.463 10.622 10.781 10.940 11.098 11.257 11.416 11.575
 [40] 11.733 11.892 12.051 12.210 12.368 12.527 12.686 12.844 13.003 13.162 13.321 13.479 13.638
 [53] 13.797 13.956 14.114 14.273 14.432 14.590 14.749 14.908 15.067 15.225 15.384 15.543 15.702
 [66] 15.860 16.019 16.178 16.337 16.495 16.654 16.813 16.971 17.130 17.289 17.448 17.606 17.765
 [79] 17.924 18.083 18.241 18.400 18.559 18.717 18.876 19.035 19.194 19.352 19.511 19.670 19.829
 [92] 19.987 20.146 20.305 20.463 20.622 20.781 20.940 21.098 21.257
attr(,"label")
[1] "Predicted values"

You may well have noticed that although we asked for standard errors of the predictions, we did not receive them. The overloaded predict() function only provides standard errors for simple linear and generalized linear models (not mixed effects models). This is one of the reasons we are exploring the use of contrasts - so that we can derive these things manually.

Use contrasts and matrix algebra to derive predictions and errors

First retrieve the fitted model parameters (intercept and slope) for the fixed effects.

> coefs <- data.lme$coef$fixed
> coefs

(Intercept)           x 
  -7875.171       3.929

Generate a model matrix appropriate for the predictor variable range. The more numbers there are over the range the smoother (less jagged) the fitted curve.

> newX <- seq(2006, 2010, l = 100)
> cc <- as.matrix(expand.grid(1, newX))
> head(cc)

     Var1 Var2
[1,]    1 2006
[2,]    1 2006
[3,]    1 2006
[4,]    1 2006
[5,]    1 2006
[6,]    1 2006

Calculate (derive) the predicted fish counts over time across all reefs.

> newY <- coefs %*% t(cc)
> newY

      [,1]  [,2] [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
[1,] 5.543 5.702 5.86 6.019 6.178 6.337 6.495 6.654 6.813 6.971  7.13 7.289 7.448 7.606 7.765 7.924
     [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31]
[1,] 8.083 8.241   8.4 8.559 8.717 8.876 9.035 9.194 9.352 9.511  9.67 9.829 9.987 10.15  10.3
     [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46]
[1,] 10.46 10.62 10.78 10.94  11.1 11.26 11.42 11.57 11.73 11.89 12.05 12.21 12.37 12.53 12.69
     [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61]
[1,] 12.84    13 13.16 13.32 13.48 13.64  13.8 13.96 14.11 14.27 14.43 14.59 14.75 14.91 15.07
     [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76]
[1,] 15.23 15.38 15.54  15.7 15.86 16.02 16.18 16.34  16.5 16.65 16.81 16.97 17.13 17.29 17.45
     [,77] [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90] [,91]
[1,] 17.61 17.77 17.92 18.08 18.24  18.4 18.56 18.72 18.88 19.03 19.19 19.35 19.51 19.67 19.83
     [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100]
[1,] 19.99 20.15  20.3 20.46 20.62 20.78 20.94  21.1  21.26

Calculate (derive) the standard errors of the predicted fish counts.

> vc <- vcov(data.lme)
> se <- sqrt(diag(cc %*% vc %*% t(cc)))
> head(se)

[1] 0.5931 0.5912 0.5893 0.5874 0.5855 0.5837

Make a data frame to hold the predicted values, standard errors and confidence intervals

> mat <- data.frame(x = newX, fit = newY[1, ], se = se, lwr = newY[1, 
+     ] + qnorm(0.025) * se, upr = newY[1, ] + qnorm(0.975) * se)
> head(mat)

     x   fit     se   lwr   upr
1 2006 5.543 0.5931 4.380 6.705
2 2006 5.702 0.5912 4.543 6.860
3 2006 5.860 0.5893 4.705 7.015
4 2006 6.019 0.5874 4.868 7.170
5 2006 6.178 0.5855 5.030 7.325
6 2006 6.337 0.5837 5.192 7.481

Finally, lets plot the data with the predicted trend and confidence intervals overlaid

> plot(y ~ x, data, pch = as.numeric(data$reef), col = "grey")
> lines(fit ~ x, mat, col = "black")
> lines(lwr ~ x, mat, col = "grey")
> lines(upr ~ x, mat, col = "grey")

Contrasts in generalized linear mixed effects models

The next demonstration, will use very similar data to the previous example, except that a Poisson error distribution would be more appropriate (due to the apparent relationship between mean and variance).

Create the artificial data

This data set consists of fish counts over time on seven reefs (you can tell it is artificial, because they have increased ;^)

> set.seed(2)
> x <- NULL
> y <- NULL
> for (i in 1:7) {
+     xx <- 2006:2010
+     yy <- round(-8300 + 4.14 * x + rpois(5, 2) * 2 + rnorm(1, 0, 3), 0) + 1
+     yy[yy < 0] <- 0
+     x <- c(x, xx)
+     y <- c(y, yy)
+ }
> data <- data.frame(x, y, reef = gl(7, 5, 35, LETTERS[1:7]))
> scatterplot(y ~ x | reef, data)

> hist(yy)

> #use a scatterplot to explore the trend for each reef
> library(car)
> scatterplot(y ~ x | reef, data)

Lets fit the generalized linear mixed effects model relating fish counts (y) to year (x). To do so, we will employ the glmmPQL() function from within the MASS package. The glmmPQL fits generalized linear mixed effects models with penalized quasi-likelihood estimation. This time we will use a Poisson distribution distribution.

> library(MASS)
> data.lme <- glmmPQL(y ~ x, random = ~1 | reef, data, family = "poisson")
> summary(data.lme)

Linear mixed-effects model fit by maximum likelihood
 Data: data 
  AIC BIC logLik
   NA  NA     NA

Random effects:
 Formula: ~1 | reef
        (Intercept) Residual
StdDev:     0.02314    1.142

Variance function:
 Structure: fixed weights
 Formula: ~invwt 
Fixed effects: y ~ x 
             Value Std.Error DF t-value p-value
(Intercept) -479.6     37.74 97  -12.71       0
x              0.2      0.02 97   12.79       0
 Correlation: 
  (Intr)
x -1    

Standardized Within-Group Residuals:
    Min      Q1     Med      Q3     Max 
-1.8834 -0.7651  0.1044  0.5940  2.0398 

Number of Observations: 105
Number of Groups: 7

Use contrasts and matrix algebra to derive predictions and errors

First retrieve the fitted model parameters (intercept and slope) for the fixed effects.

> coefs <- data.lme$coef$fixed
> coefs <- coefs
> coefs

(Intercept)           x 
  -479.5766      0.2403

Generate a model matrix appropriate for the predictor variable range. The more numbers there are over the range the smoother (less jagged) the fitted curve.

> newX <- seq(2006, 2010, l = 100)
> cc <- as.matrix(expand.grid(1, newX))
> head(cc)

     Var1 Var2
[1,]    1 2006
[2,]    1 2006
[3,]    1 2006
[4,]    1 2006
[5,]    1 2006
[6,]    1 2006

Calculate (derive) the predicted fish counts over time across all reefs. As is, these will be on the log scale, since the canonical link function for a Poisson distribution is a log function. Exponentiation will put this back in scale of the response (fish counts)

> newY <- exp(coefs %*% t(cc))
> newY

      [,1]  [,2] [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
[1,] 11.28 11.39 11.5 11.62 11.73 11.84 11.96 12.08 12.19 12.31 12.43 12.55 12.68  12.8 12.92 13.05
     [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31]
[1,] 13.18 13.31 13.44 13.57  13.7 13.83 13.97  14.1 14.24 14.38 14.52 14.66 14.81 14.95  15.1
     [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46]
[1,] 15.24 15.39 15.54 15.69 15.85    16 16.16 16.32 16.48 16.64  16.8 16.96 17.13 17.29 17.46
     [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61]
[1,] 17.63 17.81 17.98 18.15 18.33 18.51 18.69 18.87 19.06 19.24 19.43 19.62 19.81 20.01  20.2
     [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76]
[1,]  20.4  20.6  20.8    21 21.21 21.41 21.62 21.83 22.05 22.26 22.48  22.7 22.92 23.14 23.37
     [,77] [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90] [,91]
[1,]  23.6 23.83 24.06 24.29 24.53 24.77 25.01 25.25  25.5 25.75    26 26.25 26.51 26.77 27.03
     [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100]
[1,] 27.29 27.56 27.83  28.1 28.38 28.65 28.93 29.21   29.5

Calculate (derive) the standard errors of the predicted fish counts. Again, these need to be corrected if they are to be on the response scale. This is done by multiplying the standard error estimates by the exponential fitted estimates.
> vc <- vcov(data.lme) > se <- sqrt(diag(cc %*% vc %*% t(cc))) > se <- se * newY[1, ] > head(se)
[1] 0.6003 0.5988 0.5972 0.5956 0.5939 0.5922

Make a data frame to hold the predicted values, standard errors and confidence intervals

> mat <- data.frame(x = newX, fit = newY[1, ], se = se, lwr = newY[1, 
+     ] + qnorm(0.025) * se, upr = newY[1, ] + qnorm(0.975) * se)
> head(mat)

     x   fit     se   lwr   upr
1 2006 11.28 0.6003 10.11 12.46
2 2006 11.39 0.5988 10.22 12.57
3 2006 11.50 0.5972 10.33 12.67
4 2006 11.62 0.5956 10.45 12.78
5 2006 11.73 0.5939 10.56 12.89
6 2006 11.84 0.5922 10.68 13.00

Finally, lets plot the data with the predicted trend and confidence intervals overlaid

> plot(y ~ x, data, pch = as.numeric(data$reef), col = "grey")
> lines(fit ~ x, mat, col = "black")
> lines(lwr ~ x, mat, col = "grey")
> lines(upr ~ x, mat, col = "grey")

Contrasts in Single Factor ANOVA

Linear model effects parameterization for single factor ANOVA looks like:

y_i = μ + β_ix_i + ε_i
ε_i ∼ Normal(0,σ²)	#Residual 'random' effects

Motivating example

Lets say we had measured a response (e.g. weight) from four individuals (replicates) treated in one of three ways (three treatment levels).

Show code

> set.seed(1)
> baseline <- 40
> all.effects <- c(baseline, 10, 15)
> sigma <- 3  #standard deviation
> n <- 3
> reps <- 4
> FactA <- gl(n = n, k = 4, len = n * reps, lab = paste("a", 1:n, sep = ""))
> X <- as.matrix(model.matrix(~FactA))
> eps <- round(rnorm(n * reps, 0, sigma), 2)
> Response <- as.numeric(as.matrix(X) %*% as.matrix(all.effects) + 
+     eps)
> data <- data.frame(Response, FactA)
> write.table(data, file = "fish.csv", quote = F, row.names = F, sep = ",")

Download fish data set

Q4-1 Open and summarize the data

Show code

> fish <- read.table("../downloads/data/fish.csv", header = T, sep = ",", 
+     strip.white = T)
> fish

   Response FactA
1     38.12    a1
2     40.55    a1
3     37.49    a1
4     44.79    a1
5     50.99    a2
6     47.54    a2
7     51.46    a2
8     52.21    a2
9     56.73    a3
10    54.08    a3
11    59.54    a3
12    56.17    a3

Raw data

Treatment means




     FactA: a1    38.12, 40.55, 37.49, 44.79  
     FactA: a2    50.99, 47.54, 51.46, 52.21  
     FactA: a3    56.73, 54.08, 59.54, 56.17




     FactA: a1   FactA: a2   FactA: a3   
     Mean    40.24   50.55   56.63

Q4-2 Linear models (ANOVA)

y_i = μ + β_ix_i + ε_i
y = μ + β₁x₁ + β₂x₂ + β₃x₃ + ε_i	# over-parameterized model
y = β₁x₁ + β₂x₂ + β₃x₃ + ε_i	# means parameterization model
y = Intercept + β₂x₂ + β₃x₃ + ε_i	# effects parameterization model

Two alternative model parameterizations (based on treatment contrasts) would be:

Show code

> # Means parameterization<br>
> data.frame(Response, FactA, model.matrix(~-1 + FactA, data = fish))

   Response FactA FactAa1 FactAa2 FactAa3
1     38.12    a1       1       0       0
2     40.55    a1       1       0       0
3     37.49    a1       1       0       0
4     44.79    a1       1       0       0
5     50.99    a2       0       1       0
6     47.54    a2       0       1       0
7     51.46    a2       0       1       0
8     52.21    a2       0       1       0
9     56.73    a3       0       0       1
10    54.08    a3       0       0       1
11    59.54    a3       0       0       1
12    56.17    a3       0       0       1

> # Effects parameterization<br>
> data.frame(Response, FactA, model.matrix(~FactA, data = fish))

   Response FactA X.Intercept. FactAa2 FactAa3
1     38.12    a1            1       0       0
2     40.55    a1            1       0       0
3     37.49    a1            1       0       0
4     44.79    a1            1       0       0
5     50.99    a2            1       1       0
6     47.54    a2            1       1       0
7     51.46    a2            1       1       0
8     52.21    a2            1       1       0
9     56.73    a3            1       0       1
10    54.08    a3            1       0       1
11    59.54    a3            1       0       1
12    56.17    a3            1       0       1

Means parameterization

Effects parameterization

Data		Model parameters
Response	FactA	FactAa1	FactAa2	FactAa3
38.12	a1	1.00	0.00	0.00
40.55	a1	1.00	0.00	0.00
37.49	a1	1.00	0.00	0.00
44.79	a1	1.00	0.00	0.00
50.99	a2	0.00	1.00	0.00
47.54	a2	0.00	1.00	0.00
51.46	a2	0.00	1.00	0.00
52.21	a2	0.00	1.00	0.00
56.73	a3	0.00	0.00	1.00
54.08	a3	0.00	0.00	1.00
59.54	a3	0.00	0.00	1.00
56.17	a3	0.00	0.00	1.00

Data		Model parameters
Response	FactA	Intercept	FactAa2	FactAa3
38.12	a1	1.00	0.00	0.00
40.55	a1	1.00	0.00	0.00
37.49	a1	1.00	0.00	0.00
44.79	a1	1.00	0.00	0.00
50.99	a2	1.00	1.00	0.00
47.54	a2	1.00	1.00	0.00
51.46	a2	1.00	1.00	0.00
52.21	a2	1.00	1.00	0.00
56.73	a3	1.00	0.00	1.00
54.08	a3	1.00	0.00	1.00
59.54	a3	1.00	0.00	1.00
56.17	a3	1.00	0.00	1.00

Whilst the above two are equivalent for simple designs (such as that illustrated), the means parameterization (that does not include an intercept, or overall mean) is limited to simple single factor designs. Designs that involve additional predictors can only be accommodated via effects parameterization.

Lets now fit the above means and effects parameterized models and contrasts the estimated parameters:

Means parameterization	Effects parameterization
> summary(lm(Response ~ -1 + FactA, fish))$coef Estimate Std. Error t value Pr(>\|t\|) FactAa1 40.24 1.3 30.94 1.886e-10 FactAa2 50.55 1.3 38.87 2.453e-11 FactAa3 56.63 1.3 43.55 8.871e-12	> summary(lm(Response ~ FactA, fish))$coef Estimate Std. Error t value Pr(>\|t\|) (Intercept) 40.24 1.300 30.940 1.886e-10 FactAa2 10.31 1.839 5.607 3.312e-04 FactAa3 16.39 1.839 8.913 9.243e-06

Means parameterization

Effects parameterization

> summary(lm(Response ~ -1 + FactA, fish))$coef

        Estimate Std. Error t value  Pr(>|t|)
FactAa1    40.24        1.3   30.94 1.886e-10
FactAa2    50.55        1.3   38.87 2.453e-11
FactAa3    56.63        1.3   43.55 8.871e-12

> summary(lm(Response ~ FactA, fish))$coef

            Estimate Std. Error t value  Pr(>|t|)
(Intercept)    40.24      1.300  30.940 1.886e-10
FactAa2        10.31      1.839   5.607 3.312e-04
FactAa3        16.39      1.839   8.913 9.243e-06

The model on the left omits the intercept and thus estimates the treatment means. The parameters estimated by the model on the right are the mean of the first treatment level and the differences between each of the other treatment means and the first treatment mean. If there were only two treatment levels or first treatment was a control (two which you wish to compare the other treatment means), then the effects model is providing what the statisticians are advocating - estimates of the effect sizes as well as estimates of precision (standard error from which confidence intervals can easily be derived).

Derived parameters - simple contrasts

These parameters represent model coefficients that can be substituted back into a linear model. They represent the set of linear constants in the linear model. We can subsequently multiply these parameters to produce alternative linear combinations in order to estimate other values from the linear model.

y_i = μ + β_ix_i

y_i = μ + β₂x₂ + β₂x₂

y_i = 40.2375 × a1 + 10.3125 × a2 + 16.3925 × a3

For example, if we wanted to estimate the effect size for the difference between treatments a2 and a3, we could make a1=0, a2=1 and a3=-1. That is, if we multiplied our model parameters by a contrast matrix of 0,1,-1;

y_i =

β₂

β₃

The matrix;

is a matrix of contrast coefficients. Each row corresponds to a column in the parameter matrix (vector). The contrast matrix can have multiple columns, each corresponding to a different derived comparison.

y_i = 40.2375 × a1 + 10.3125 × a2 + 16.3925 × a3

y_i = 40.2375 × 0 + 10.3125 × 1 + 16.3925 × -1

y_i = 10.3125 - 16.3925

y_i = -6.080

Planned comparisons

Some researchers might be familiar with the use of contrasts in performing so called 'planned comparisons' or 'planned contrasts'. This involves redefining the actual parameters used in the linear model to reflect more meaningful comparisons amongst groups. For example, rather than compare each group mean against the mean of one (the first) treatment group, we many wish to include a comparison of the first treatment group (perhaps a control group) against the average two other treatment groups.

> contrasts(fish$FactA) <- cbind(`a1 vs (a2 & a3)` = c(1, -0.5, -0.5))
> summary(aov(Response ~ FactA, fish), split = list(FactA = list(1, 
+     2)))

            Df Sum Sq Mean Sq F value  Pr(>F)    
FactA        2    549     275    40.6 3.1e-05 ***
  FactA: C1  1    475     475    70.3 1.5e-05 ***
  FactA: C2  1     74      74    10.9  0.0091 ** 
Residuals    9     61       7                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> summary(aov(Response ~ FactA, fish), split = list(FactA = list(`a1 vs (a2 & a3)` = 1)))

                         Df Sum Sq Mean Sq F value  Pr(>F)    
FactA                     2    549     275    40.6 3.1e-05 ***
  FactA: a1 vs (a2 & a3)  1    475     475    70.3 1.5e-05 ***
Residuals                 9     61       7                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> contrasts(fish$FactA) <- contr.treatment

> contrasts(fish$FactA) <- cbind(`a2 vs a3` = c(0, 1, -1))
> summary(aov(Response ~ FactA, fish), split = list(FactA = list(1, 
+     2)))

            Df Sum Sq Mean Sq F value  Pr(>F)    
FactA        2    549     275    40.6 3.1e-05 ***
  FactA: C1  1     74      74    10.9  0.0091 ** 
  FactA: C2  1    475     475    70.3 1.5e-05 ***
Residuals    9     61       7                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> summary(aov(Response ~ FactA, fish), split = list(FactA = list(`a2 vs a3` = 1)))

                  Df Sum Sq Mean Sq F value  Pr(>F)    
FactA              2    549   274.7    40.6 3.1e-05 ***
  FactA: a2 vs a3  1     74    73.9    10.9  0.0091 ** 
Residuals          9     61     6.8                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> contrasts(fish$FactA) <- contr.treatment

In so doing, this technique provides a relatively easy way to test other non-standard hypotheses (linear effects).

However, a greater understanding of contrasts has even greater benefits when it comes to calculations of effect sizes, marginal means and their associated confidence intervals - features that are considered of greater interpretive value than hypothesis tests. Contrasts can also be used to derive a range of other derived parameters from the fitted model.

Calculating treatment means (and confidence intervals) from the effects model

A simple way to illustrate the use of contrasts to derive new parameters from a linear model is in the estimation of population treatment means from the effects model.

y_i =

β₂

β₃

1	1	1
0	1	0
0	0	1

> contrasts(fish$FactA) <- contr.treatment
> fish.lm <- lm(Response ~ FactA, fish)
> cmat <- cbind(c(1, 0, 0), c(1, 1, 0), c(1, 0, 1))
> colnames(cmat) <- c("a1", "a2", "a3")
> fish.mu <- fish.lm$coef %*% cmat
> fish.mu

        a1    a2    a3
[1,] 40.24 50.55 56.63

> fish.se <- sqrt(diag(t(cmat) %*% vcov(fish.lm) %*% cmat))
> fish.se

 a1  a2  a3 
1.3 1.3 1.3

> fish.ci <- as.numeric(fish.mu) + outer(fish.se, qt(c(0.025, 0.975), 
+     df = fish.lm$df.resid))
> fish.ci

    [,1]  [,2]
a1 37.30 43.18
a2 47.61 53.49
a3 53.69 59.57

Calculating all the pairwise effects (and confidence intervals) from the effects model

y_i =

β₂

β₃

0	0	0
1	0	-1
0	1	1

> contrasts(data$FactA) <- contr.treatment
> data.lm <- lm(Response ~ FactA, data)
> cmat <- cbind(c(0, 1, 0), c(0, 0, 1), c(0, -1, 1))
> colnames(cmat) <- c("a1 vs a2", "a1 vs a3", "a2 vs a3")
> data.mu <- data.lm$coef %*% cmat
> data.mu

     a1 vs a2 a1 vs a3 a2 vs a3
[1,]    10.31    16.39     6.08

> data.se <- sqrt(diag(t(cmat) %*% vcov(data.lm) %*% cmat))
> data.se

a1 vs a2 a1 vs a3 a2 vs a3 
   1.839    1.839    1.839

> data.ci <- as.numeric(data.mu) + outer(data.se, qt(c(0.025, 0.975), 
+     df = data.lm$df.resid))
> data.ci

           [,1]  [,2]
a1 vs a2  6.152 14.47
a1 vs a3 12.232 20.55
a2 vs a3  1.920 10.24

Contrasts in factorial ANOVA

Motivating example

We now extend this example such that each individual fish was subjected to one of five levels of one treatment and one of three levels of another treatment. There were four replicate fish in each treatment combination.

Note, this is an artificial data set. Most of the data sets used in this Workshop series have been real data sets. Real datasets are useful as they provide 'real world' examples of how to apply statistical principles and tools. However, for understanding some techniques, fabricated data allows us to explore techniques under controlled conditions. Furthermore, they allow us to compare statistical outcomes and patterns derived from statistical techniques to the 'truth'. Recall that the purpose of statistics is to 'estimate' population parameters. Fabricating data allows us to know the true values of the parameters that we are trying to estimate and thus allow us to evaluate the accuracy of our outcomes and understanding. I apologise if the use of artificial data causes you any pain or personal hardship.

In the creation of these data, I will be defining the 'true' (theoretical) effect sizes and level of 'noise'.

Show code

> set.seed(1)
> n.FactA <- 5
> n.FactB <- 3
> nsample <- 12
> n <- n.FactA * nsample
> #create factors
> FactA <- gl(n = n.FactA, k = nsample, length = n, lab = paste("a", 
+     1:n.FactA, sep = ""))
> FactB <- gl(n = n.FactB, k = nsample/n.FactB, length = n, lab = paste("b", 
+     1:n.FactB, sep = ""))
> #Effects
> baseline <- 40
> FactA.effects <- c(-10, -5, 5, 10)
> FactB.effects <- c(5, 10)
> interaction.effects <- c(-2, 3, 0, 4, 4, 0, 3, -2)
> all.effects <- c(baseline, FactA.effects, FactB.effects, interaction.effects)
> X <- as.matrix(model.matrix(~FactA * FactB))

Theoretical means

We will confirm the theoretical means of this design matrix (based on the defined effects).

> RESPONSE1 <- as.numeric(as.matrix(X) %*% as.matrix(all.effects))
> theoryMeans <- tapply(RESPONSE1, list(FactA, FactB), mean)
> theoryMeans

   b1 b2 b3
a1 40 45 50
a2 30 33 44
a3 35 43 45
a4 45 50 58
a5 50 59 58

> library(xtable)
> print(xtable(theoryMeans), type = "html", sanitize.rownames.function = function(x) paste("<b>", 
+     x, "</b>"), html.table.attributes = list("border='3' cellpadding='5' cellspacing='0'"))




     b1   b2   b3   
     a1    40.00   45.00   50.00  
     a2    30.00   33.00   44.00  
     a3    35.00   43.00   45.00  
     a4    45.00   50.00   58.00  
     a5    50.00   59.00   58.00

	b1	b2	b3
a1	40.00	45.00	50.00
a2	30.00	33.00	44.00
a3	35.00	43.00	45.00
a4	45.00	50.00	58.00
a5	50.00	59.00	58.00

And now we can add noise to the response to create data with variability. We will add noise by adding value to each observation that is randomly drawn from a normal distribution with a mean of zero and a standard deviation of 3.

> sigma <- 3
> eps <- rnorm(n, 0, sigma)
> X <- as.matrix(model.matrix(~FactA * FactB))
> RESPONSE <- as.numeric(as.matrix(X) %*% as.matrix(all.effects) + 
+     eps)
> data <- data.frame(RESPONSE, FactB, FactA)
> write.table(data, file = "fish1.csv", quote = F, row.names = F, sep = ",")

Lets just examine the first six rows of these data to confirm that a little noise has been added

> head(data)

  RESPONSE FactB FactA
1    38.12    b1    a1
2    40.55    b1    a1
3    37.49    b1    a1
4    44.79    b1    a1
5    45.99    b2    a1
6    42.54    b2    a1

Download fish1 data set

Q5-1 Open and summarize the data

Show code

> fish1 <- read.table("fish1.csv", header = T, sep = ",", strip.white = T)
> head(fish1)

  RESPONSE FactB FactA
1    38.12    b1    a1
2    40.55    b1    a1
3    37.49    b1    a1
4    44.79    b1    a1
5    45.99    b2    a1
6    42.54    b2    a1

Simple sample means

We will now calculate the sample means of each FactA by FactB combination. These are called 'cell means' as they are the means of each cell in a table in which the levels of FactA and FactB are the rows and columns (as was depicted in the table above for the theoretical means). If the samples are collected in an unbiased manner (whereby the collected samples reflect the populations), then the cell means should roughly approximate the true means of the underlying populations from whcih the observations were drawn.

Raw sample group means:

> groupMeans <- with(fish1, tapply(RESPONSE, list(FactA, FactB), mean))
> groupMeans

      b1    b2    b3
a1 40.24 45.55 51.63
a2 28.68 34.76 43.84
a3 34.20 43.90 43.91
a4 46.06 50.63 57.80
a5 50.42 60.97 57.27

And the differences between the theoretical means (which DID NOT contain any observational variance) and these raw sample group means:

> theory_groupMeans <- theoryMeans - groupMeans
> theory_groupMeans

        b1      b2      b3
a1 -0.2376 -0.5511 -1.6290
a2  1.3170 -1.7571  0.1603
a3  0.7971 -0.8968  1.0936
a4 -1.0572 -0.6268  0.1974
a5 -0.4161 -1.9689  0.7325

Here is a side-by-side comparison of the theoretical means, the raw sample means and the cell differences:

Theoretical cell means

Sample cell means

Differences




     b1   b2   b3   
     a1    40.00   45.00   50.00  
     a2    30.00   33.00   44.00  
     a3    35.00   43.00   45.00  
     a4    45.00   50.00   58.00  
     a5    50.00   59.00   58.00




     b1   b2   b3   
     a1    40.24   45.55   51.63  
     a2    28.68   34.76   43.84  
     a3    34.20   43.90   43.91  
     a4    46.06   50.63   57.80  
     a5    50.42   60.97   57.27




     b1   b2   b3   
     a1    -0.24   -0.55   -1.63  
     a2    1.32   -1.76   0.16  
     a3    0.80   -0.90   1.09  
     a4    -1.06   -0.63   0.20  
     a5    -0.42   -1.97   0.73

The sum of the deviations (differences) of the sample means from the theoretical means is:

> sum(theory_groupMeans)

[1] -4.843

Modeled population means

Alternatively, we can estimate the cell means based on a full multiplicative linear model.

y_ij = µ + α_i + β_j + αβ_ij

Note that this model is over-parameterized - that is, there are more parameters to estimate than there are groups of data from which to independently estimate them.

The effects parameterization of this model via treatment contrasts incorporates the mean of the first cell (µ) followed by a series of differences. In all, there are 15 groups (five levels of FactA multiplied by three levels of FactB) and correspondingly, there are 15 parameters being estimated:

y_ij = µ + β1 + β2 + β3 + β4 + β5 + β6 + β7 + β8 + β9 + β10 + β11 + β12 + β13 + β14

These parameters correspond to the following:




  Model term   Summarized model parameter
(name given in R)   
Interpretation   
    µ   (Intercept)   (a1b1)  
    β1   FactAa2   (a2b1 - a1b1)  
    β2   FactAa3   (a3b1 - a1b1)  
    β3   FactAa4   (a4b1 - a1b1)  
    β4   FactAa5   (a5b1 - a1b1)  
    β5   FactBb2   (a1b2 - a1b1)  
    β6   FactBb3   (a1b3 - a1b1)  
    β7   FactAa2:FactBb2   (a2b1+a1b2) - (a1b1+a2b2)  
    β8   FactAa3:FactBb2   (a3b1+a1b2) - (a1b1+a3b2)  
    β9   FactAa4:FactBb2   (a4b1+a1b2) - (a1b1+a4b2)  
    β10   FactAa5:FactBb2   (a5b1+a1b2) - (a1b1+a5b2)  
    β11   FactAa2:FactBb3   (a2b1+a1b3) - (a1b1+a2b3)  
    β12   FactAa3:FactBb3   (a3b1+a1b3) - (a1b1+a3b3)  
    β13   FactAa4:FactBb3   (a4b1+a1b3) - (a1b1+a4b3)  
    β14   FactAa5:FactBb3   (a5b1+a1b3) - (a1b1+a5b3)

In R, the linear model is fitted via the lm() function.

> fish1.lm <- lm(RESPONSE ~ FactA * FactB, fish1)

Whilst these parameters satisfy the model parameterization, interpretating and utilizing the estimated parameters can take some getting used to. Lets begin by summarizing the estimates of the 15 model parameters.

> fish1.lm$coef

    (Intercept)         FactAa2         FactAa3         FactAa4         FactAa5         FactBb2 
        40.2376        -11.5546         -6.0348          5.8196         10.1785          5.3135 
        FactBb3 FactAa2:FactBb2 FactAa3:FactBb2 FactAa4:FactBb2 FactAa5:FactBb2 FactAa2:FactBb3 
        11.3914          0.7606          4.3804         -0.7439          5.2393          3.7653 
FactAa3:FactBb3 FactAa4:FactBb3 FactAa5:FactBb3 
        -1.6879          0.3541         -4.5400

But this is where the good stuff begins. As the parameters are defined to encapsulate all of the patterns amongst the groups, the linear combination of parameters can now be used to calculate a whole range of related derived estimates (and their standard deviations etc) including:

population cell means
marginal means (row and column means)
effect sizes (differences between marginal means)
effect sizes (deviations from overall mean)

In order to achieve this wizardry, we must indicate which combination of the parameter estimates are relevant for each of the derived estimates required. And this is where contrast coefficients come in. The contrast coefficients are a set of numbers (equal in length to the number of model parameters) that determine the combinations and weights of the model parameters that are appropriate for each of the derived estimates. Hence when the set of model parameters are multiplied by a single set of contrast coefficients and then summed up, the result is a new derived estimate.

By way of example, the following contrast would specify that the derived estimate should include only (and all of) the first model parameter (and thus according to the table above, estimate the mean of cell FactA level a1 and FactB level b1). In R, the %*% is a special operator that performs matrix multiplication (each element of one matrix is multiplied by the corresponding elements of the other matrix) and returns their sum.

> fish1.lm$coef

    (Intercept)         FactAa2         FactAa3         FactAa4         FactAa5         FactBb2 
        40.2376        -11.5546         -6.0348          5.8196         10.1785          5.3135 
        FactBb3 FactAa2:FactBb2 FactAa3:FactBb2 FactAa4:FactBb2 FactAa5:FactBb2 FactAa2:FactBb3 
        11.3914          0.7606          4.3804         -0.7439          5.2393          3.7653 
FactAa3:FactBb3 FactAa4:FactBb3 FactAa5:FactBb3 
        -1.6879          0.3541         -4.5400

> c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

 [1] 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0

> fish1.lm$coef %*% c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

      [,1]
[1,] 40.24

In the above, we indicated that only the first model parameter should be included in deriving the specified estimate, as all the others parameters are multiplied by zeros.

We can calculate multiple derived estimates at once, by multiplying the model parameters by multiple sets of contrast coefficients (represented as a matrix). That is, matrix multiplication. Hence to also calculate the cell mean of FactA level a2 and FactB b1 (a2b1), we would specify the first (a1b1 cell mean) and second (difference between a2b1 and a1b1 cell means) contrasts as ones (1) and the others zeros.

The cbind() function binds together vectors into the columns of a matrix.

> mat <- cbind(a1b1 = c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 
+     a2b2 = c(1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
> mat

      a1b1 a2b2
 [1,]    1    1
 [2,]    0    1
 [3,]    0    0
 [4,]    0    0
 [5,]    0    0
 [6,]    0    0
 [7,]    0    0
 [8,]    0    0
 [9,]    0    0
[10,]    0    0
[11,]    0    0
[12,]    0    0
[13,]    0    0
[14,]    0    0
[15,]    0    0

> fish1.lm$coef %*% mat

      a1b1  a2b2
[1,] 40.24 28.68

At this point we could go on and attempt to work out which combinations of model parameters would be required to define each of the cell means. Whilst it is reasonably straight forward to do this for cells a1b1, a2b1, a3b1, a4b1, a5b1, a1b2, a1b3 (those along the top and down the left column), it requires more mental gymnastics to arrive at the entire set. However, there is a function in R called model.matrix() that helps us achieve this in a relatively pain-free manner. The model.matrix() function gives the design matrix that is used in linear modelling to indicate group membership of categorical variables and in general provides the link between the observed data (response) and the linear combination of predictor variables.

As just eluded the means

to, model.matrix relates each individual observation. However, we are not interested in estimating ('predicting') each observation, we want to predict of groups of observations. Hence, we need to reduce the matrix down to the level of the groups. id="Q5-2">Show code for entire design matrix

	mu	se	lwr	upr
a1b1	40.24	1.35	37.52	42.96
a1b2	45.55	1.35	42.83	48.27
a1b3	51.63	1.35	48.91	54.35
a2b1	28.68	1.35	25.96	31.40
a2b2	34.76	1.35	32.04	37.48
a2b3	43.84	1.35	41.12	46.56
a3b1	34.20	1.35	31.48	36.92
a3b2	43.90	1.35	41.18	46.62
a3b3	43.91	1.35	41.19	46.63
a4b1	46.06	1.35	43.34	48.78
a4b2	50.63	1.35	47.91	53.35
a4b3	57.80	1.35	55.08	60.52
a5b1	50.42	1.35	47.70	53.14
a5b2	60.97	1.35	58.25	63.69
a5b3	57.27	1.35	54.55	59.99

Workshop 11 - The power of Contrasts

Context

Contrasts in prediction in regression analyses

Download Fertilizer data set

Use the predict() function to extract the line of best fit and standard errors

Use contrasts and matrix algebra to derive predictions and errors

Contrasts in linear mixed effects models

Create the artificial data

Use the predict() function to extract the line of best fit and standard errors

Use contrasts and matrix algebra to derive predictions and errors

Contrasts in generalized linear mixed effects models

Create the artificial data

Use contrasts and matrix algebra to derive predictions and errors

Contrasts in Single Factor ANOVA

Motivating example

Download fish data set

Derived parameters - simple contrasts

Planned comparisons

Calculating treatment means (and confidence intervals) from the effects model

Calculating all the pairwise effects (and confidence intervals) from the effects model

Contrasts in factorial ANOVA

Motivating example

Theoretical means

Download fish1 data set

Simple sample means

Modeled population means

Conditions

Standard deviations, standard errors and confidence intervals

Interaction plot

Marginal means

Effect sizes - differences between means

Effect sizes - for simple main effects

Effect sizes - for marginal effects (no interactions)

Simple main effects - effects associated with one factor at specific levels of another factor

Model term	Summarized model parameter (name given in R)	Interpretation
µ	(Intercept)	(a1b1)
β1	FactAa2	(a2b1 - a1b1)
β2	FactAa3	(a3b1 - a1b1)
β3	FactAa4	(a4b1 - a1b1)
β4	FactAa5	(a5b1 - a1b1)
β5	FactBb2	(a1b2 - a1b1)
β6	FactBb3	(a1b3 - a1b1)
β7	FactAa2:FactBb2	(a2b1+a1b2) - (a1b1+a2b2)
β8	FactAa3:FactBb2	(a3b1+a1b2) - (a1b1+a3b2)
β9	FactAa4:FactBb2	(a4b1+a1b2) - (a1b1+a4b2)
β10	FactAa5:FactBb2	(a5b1+a1b2) - (a1b1+a5b2)
β11	FactAa2:FactBb3	(a2b1+a1b3) - (a1b1+a2b3)
β12	FactAa3:FactBb3	(a3b1+a1b3) - (a1b1+a3b3)
β13	FactAa4:FactBb3	(a4b1+a1b3) - (a1b1+a4b3)
β14	FactAa5:FactBb3	(a5b1+a1b3) - (a1b1+a5b3)