Contrasts

Making pointed comparisons between treatment levels in ANOVA

Kris Sankaran true
09-23-2021

Readings 3.5, Rmarkdown

  1. When we reject the null in ANOVA, we know at least one of the treatments deviates from the global average. But which one(s)? Contrasts address this question. A contrast is a linear combination of the means, \[ \Gamma(c)=\sum_{i=1}^{a} c_{i} \mu_{i} \] For any particular \(c\), we test \[ \begin{aligned} &H_{0}: \Gamma(c)=0 \\ &H_{1}: \Gamma(c) \neq 0 \end{aligned} \]
Is there a difference between group 2 and the average of groups 1 and 3?

Figure 1: Is there a difference between group 2 and the average of groups 1 and 3?

  1. To motivate this, suppose we have 4 different means, \(\mu_1, \mu_2, \mu_3, \mu_4\),
    • c = (1, 1, -1, -1): Are the first two means different from the last two, on average?
    • c = (1, -1, 0, 0): Are the first two means equal to each other?
    • etc.

Testing Contrasts

  1. Recall the hypothesis testing recipe. We need, (a) a test statistic and (b) a reference distribution for that statistic. Our best guess at \(\mu_i\) is \(\bar{y}_i\), so a reasonable statistic is, \[ \hat{\Gamma}(c)=\sum_{i=1}^{a} c_{i} \bar{y}_{i} \]

  2. How will we find its reference distribution? Under the null, this statistic is normally distributed with mean 0 and variance, \[ \begin{aligned} \operatorname{Var}(\hat{\Gamma}(c)) &=\sum_{i=1}^{a} c_{i}^{2} \operatorname{Var}\left(\bar{y}_{i}\right) \\ &=\frac{\sigma^{2}}{n} \sum_{i=1}^{a} c_{i}^{2} \end{aligned} \] Standardizing our original statistic, we obtain,

    \[\begin{aligned} \frac{\hat{\Gamma}(c)}{\sqrt{\operatorname{Var}(\hat{\Gamma}(c))}} &=\frac{\sum_{i=1}^{a} c_{i} \bar{y}_{i}}{\sqrt{\frac{\sigma^{2}}{n} \sum_{i=1}^{a} c_{i}^{2}}} \\ & \approx \frac{\sum_{i=1}^{a} c_{i} \bar{y}_{i}}{\sqrt{\frac{\hat{\sigma}^{2}}{n} \sum_{i=1}^{a} c_{i}^{2}}} \end{aligned}\]

  3. To estimate \(\sigma^2\), we can use \(\hat{\sigma}^2 := MS_E\). This is a good choice, because it remains valid even when the null is untrue.

  4. Since we plugged-in the estimate \(\hat{\sigma}^2\), we have divided our normal distribution by the square root of a chi-square. Therefore, the reference distribution is a t-distribution with \(N - a\) df.

Confidence Intervals for Contrasts

  1. If we make the same computations as above, but without assuming that the null is true, we would find that, \[ \mathbf{P}\left(\frac{\sum_{i=1}^{a} c_{i} \bar{y}_{i}-\sum_{i=1}^{a} c_{i} \mu_{i}}{\sqrt{\frac{\hat{\sigma}^{2}}{n} \sum_{i=1}^{a} c_{i}^{2}}} \in\left[t_{\text {left }}, t_{\text {right }}\right]\right)=0.95 \] where we choose \(t_{\text{left}}\) and \(t_{\text{right}}\) to be the 0.025 and 0.975 quantiles of a \(t\)-distribution with \(N - a\) df.

  2. The resulting confidence interval is, \[ \left[\sum_{i=1}^{a} c_{i} \bar{y}_{i}-t_{\text {right }} \sqrt{\frac{\hat{\sigma}^{2}}{n} \sum_{i=1}^{a} c_{i}^{2}}, \sum_{i=1}^{a} c_{i} \bar{y}_{i}+t_{\text {right }} \sqrt{\left.\frac{\hat{\sigma}^{2}}{n} \sum_{i=1}^{a} c_{i}^{2}\right]}\right. \] This is an explicit formula that you can use in your computations, but don’t let the complexity of the symbols here confuse you. Returning to our original definitions, this is just, \[ \left[\hat{\Gamma}(c)-t_{\mathrm{left}} \sqrt{\widehat{\operatorname{Var}}(\hat{\Gamma}(c))}, \hat{\Gamma}(c)+t_{\mathrm{left}} \sqrt{\widehat{\operatorname{Var}}(\hat{\Gamma}(c))}\right] \] which is our point estimate plus and minus \(t_{\text{left}}\) standard deviations (we’re writing \(\hat{\text{Var}}\) instead of \(\text{Var}\) because we don’t know the true variance and have plugged in the estimate \(\hat{\sigma}^2)\).

Code Example

  1. We will continue the etch rate example. Let’s refit the same model from before.
library(readr)
etch_rate <- read_csv("https://uwmadison.box.com/shared/static/vw3ldbgvgn7rupt4tz3ditl1mpupw44h.csv")
etch_rate$power <- as.factor(etch_rate$power) # consider as discrete levels
fit <- lm(rate ~ power, data = etch_rate)
  1. Let’s say we’re interested in the contrast between power levels of 160 and 180. We have to encode the contrast as a vector, and then pass it to the fit.contrast function. The choice c(1, -1, 0, 0) comes from the fact that we want the difference between the first two power levels (160 and 180) and are ignoring all the rest (the two 0’s at the end).
library(gmodels)
contrast <- c(1, -1, 0, 0)
fit.contrast(fit, "power", contrast)

                     Estimate Std. Error   t value    Pr(>|t|)
power c=( 1 -1 0 0 )    -36.2   11.55335 -3.133289 0.006416224
attr(,"class")
[1] "fit_contrast"
  1. We can get a confidence interval using the conf.int parameter.
fit.contrast(fit, "power", contrast, conf.int = 0.95)

                     Estimate Std. Error   t value    Pr(>|t|)
power c=( 1 -1 0 0 )    -36.2   11.55335 -3.133289 0.006416224
                      lower CI  upper CI
power c=( 1 -1 0 0 ) -60.69202 -11.70798
attr(,"class")
[1] "fit_contrast"