Example

• The data set chickwts is available in R, and on the course webpage.
• Newly hatched chicks were randomly allocated into six groups, and each group was given a different feed supplement.
• Their weights in grams after six weeks are given along with feed types, i.e. we have a sample with corresponding measurements of 2 variables:
  • weight: a numeric variable giving the chick weight.
  • feed: a factor giving the feed type.
• Always start with some graphics:

library(mosaic)
gf_boxplot(weight ~ feed, data = chickwts)

The ANOVA Model

• We measure the response \(y\) which in this case is weight.
• We want to study the effect of the factor \(x\) on \(y\). In this case \(x=\)feed and divides the sample in \(g=6\) groups.
• The mean responses within the groups are denoted \(\mu_1,\mu_2,\ldots,\mu_g\).
• We will assume that
  • \(y=\mu_x+\epsilon\), when \(y\) is a response in group \(x\)
  • \(\epsilon\) are a sample from a population with mean zero and standard deviation \(\sigma\).
  • The standard deviation for the population in each group is the same and equals \(\sigma\)
  • The response variable, \(y\), is normal distributed within each group.
• The ANOVA test is a test of independence between the quantitative response variable and the qualitative explanatory variables.

Estimates

• Least squares estimates for population means \(\widehat\mu_x\) is given by the average of the response measurements in group \(x\).
• For a given measured response \(y\) we let \(\widehat y\) denote the model’s prediction of \(y\), i.e. \[\widehat y = \widehat\mu_x\] if \(y\) is a response for an observation in group \(x\).
• We use mean to find the mean, for each group:

mean(weight ~ feed, data = chickwts)

##    casein horsebean   linseed  meatmeal   soybean sunflower 
##  323.5833  160.2000  218.7500  276.9091  246.4286  328.9167

• We can e.g. see that \(\widehat y=323.6\), when feed=casein but \(\widehat y=160.2\), when feed=horsebean.
• Is it a significant difference?

Contrast coding

• In many cases there is a group corresponding to “no treatment” and we are interested in the effect of different treatments.
• In this example we only have different feeds, which are sorted in lexicographical order by R, so casein is the reference.
• We can specify the model via:
  • Intercept corresponding to the mean response for the reference (casein).
  • For each of the other groups we have a contrast, which measures the difference between the mean value for that group and the reference group.
• For a given contrast we can calculate standard error, t-score and p-value, and thereby investigate whether there is a difference between this group and the reference group.
• In Agresti this is referred to as using dummy variables.

Example

model <- lm(weight ~ feed, data = chickwts)
summary(model)

## 
## Call:
## lm(formula = weight ~ feed, data = chickwts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -123.909  -34.413    1.571   38.170  103.091 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    323.583     15.834  20.436  < 2e-16 ***
## feedhorsebean -163.383     23.485  -6.957 2.07e-09 ***
## feedlinseed   -104.833     22.393  -4.682 1.49e-05 ***
## feedmeatmeal   -46.674     22.896  -2.039 0.045567 *  
## feedsoybean    -77.155     21.578  -3.576 0.000665 ***
## feedsunflower    5.333     22.393   0.238 0.812495    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 54.85 on 65 degrees of freedom
## Multiple R-squared:  0.5417, Adjusted R-squared:  0.5064 
## F-statistic: 15.36 on 5 and 65 DF,  p-value: 5.936e-10

• We get information about contrasts and their significance:
• Intercept corresponding to casein has weight different from zero (\(p < 2\times 10^{-16}\)).
• Weight difference between casein and horsebean is extremely significant (p=\(2\times 10^{-9}\)).
• There is no significant weight difference between casein and sunflower (p=\(81\)%).

Graphical representation of models

• We have two alternative explanations of the data.
• Simple model with one parameter (mean): “The feed type doesn’t matter. The weight is just random around a common mean value”.
• Complex model with six parameters (means): “The feed type is important. For each feed type we get a different mean value and the weights are random around these values.”

Hypotheses and test statistic

• Is the complex model significantly better (i.e. is there any effect of the explanatory grouping variable)? We can write the corresponding hypotheses in two different ways \[H_0: \mu_1 = \mu_2 = \ldots=\mu_g \quad \mbox{against} \quad H_a: \mbox{ At least 2 of the population means are different}\]
• Alternatively \[H_0: \mbox{ All contrasts are equal to zero. } \quad H_a: \mbox{ At least one contrast is non-zero}.\]
• We will (indirectly) use \(R^2\) to do the test. If it is large, the complex model has good predictive power compared to the simple model. To judge significance we use \[F_{obs} = \frac{(n-g)R^2}{(g-1)(1-R^2)} = \frac{(TSS-SSE)/(g-1)}{SSE/(n-g)}.\]
• Large values of \(R^2\) implies large values of \(F_{obs}\), which points to the alternative hypothesis.
• I.e. when we have calculated the observed value \(F_{obs}\), then we have to find the probability that a new experiment would result in a larger value.

Interpretation of \(F\) statistic - Variance between/within groups

• It can be shown that the numerator of \(F_{obs}\) is a measure of the variance between the groups, i.e. how much “boxes” vary around the total average (the red line).

• Likewise it can be shown the denominator of \(F_{obs}\) is a measure for the variance within groups, i.e. how “tall” the boxes in the boxplot are.

• If the boxes’ deviations from the red line are to be explained by randomness, then the two types of variances should be of same magnitude. This is measured by the F-test statistic, which can be stated as \[F_{obs} = \frac{\mbox{variance between groups}}{\mbox{variance within groups}}\]

Example

model <- lm(weight ~ feed, data = chickwts)
summary(model)

## 
## Call:
## lm(formula = weight ~ feed, data = chickwts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -123.909  -34.413    1.571   38.170  103.091 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    323.583     15.834  20.436  < 2e-16 ***
## feedhorsebean -163.383     23.485  -6.957 2.07e-09 ***
## feedlinseed   -104.833     22.393  -4.682 1.49e-05 ***
## feedmeatmeal   -46.674     22.896  -2.039 0.045567 *  
## feedsoybean    -77.155     21.578  -3.576 0.000665 ***
## feedsunflower    5.333     22.393   0.238 0.812495    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 54.85 on 65 degrees of freedom
## Multiple R-squared:  0.5417, Adjusted R-squared:  0.5064 
## F-statistic: 15.36 on 5 and 65 DF,  p-value: 5.936e-10

• The last line gives us the value of \(F_{obs} = 15.36\) and the corresponding \(p\)-value (\(5.9 \times 10^{-10}\)). Clearly there is a significant difference between the types of feed.

Additive effects

• The data set ToothGrowth is available in R and on the webpage. For more info about this data, use ?ToothGrowth.

• The data describes the tooth length in guinea pigs where some recvieve vitamin C treatment and others are given orange juice in different dosage.

• A total of \(60\) observations on 3 variables.
  • [,1] len The tooth length
  • [,2] supp The type of the supplement (OJ or VC)
  • [,3] dose The dosage (LO, ME, HI)
• We will study the response len with the predictors supp and dose.
• At first we look at the model with additive effects
  • len=\(\mu\) + "effect of supp"+ "effect of dose" + error
• This is also called the main effects model since it does not contain interaction terms.
• The parameter \(\mu\) corresponds to the Intercept and is the mean tooth length in the reference group (supp OJ, dose LO).
• The effect of supp is the difference in mean when changing from OJ to VC.

• The effect of dose is the difference in mean when changing from LO to eitherME or HI.

Dummy coding

• Let us introduce dummy variables:
  • \(s_C=1\) if supp VC and zero otherwise.
  • \(d_M=1\) if dose is ME and zero otherwise.
  • \(d_H=1\) if dose is HI and zero otherwise.
• Then we state the model \[\mbox{length}=\mu+\beta_1 s_C+\beta_2 d_M+\beta_3 d_H + \mbox{error} .\]
• Interpretation:
  • \(\mu\) is the expected tooth length when supp is OJ and dose is LO (\(s_C=d_M=d_H=0)\)).
  • \(\beta_1\) is the effect of supplement OJ to VC (\(s_C=1\)).
  • \(\beta_2\) is the effect of increasing dosage from LO to ME (\(d_M=1\)).
  • \(\beta_3\) is the effect of increasing dosage from LO to HI (\(d_H=1\)).

Main effect model in R

• The main effects model is fitted by

MainEff <- lm(len ~ supp + dose, data = ToothGrowth)
summary(MainEff)

## 
## Call:
## lm(formula = len ~ supp + dose, data = ToothGrowth)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.085 -2.751 -0.800  2.446  9.650 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  12.4550     0.9883  12.603  < 2e-16 ***
## suppVC       -3.7000     0.9883  -3.744 0.000429 ***
## doseME        9.1300     1.2104   7.543 4.38e-10 ***
## doseHI       15.4950     1.2104  12.802  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.828 on 56 degrees of freedom
## Multiple R-squared:  0.7623, Adjusted R-squared:  0.7496 
## F-statistic: 59.88 on 3 and 56 DF,  p-value: < 2.2e-16

• The model has 4 parameters.
• The \(F\) test at the end compares with the (null) model with only one overall mean parameter. Does it seem like supp and dose has an additive effect?

Testing effect of supp

• Alternative model without effect of supp:

doseEff <- lm(len ~ dose, data = ToothGrowth)

• We can compare \(R^2\) to see if doseEff (Model 1) is sufficent to explain the data compared to MainEff (Model 2). This is done by converting to \(F\)-statistic: \[ F_{obs} = \frac{(R_2^2 - R_1^2)/(df_1 - df_2)}{(1 - R_2^2)/df_2} = \frac{(SSE_1 - SSE_2)/(df_1 - df_2)}{(SSE_2)/df_2}. \]
• In R the calculations are done using anova:

anova(doseEff, MainEff)

## Analysis of Variance Table
## 
## Model 1: len ~ dose
## Model 2: len ~ supp + dose
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1     57 1025.78                                  
## 2     56  820.43  1    205.35 14.017 0.0004293 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Testing effect of dose

• Alternative model without effect of dose:

suppEff <- lm(len ~ supp, data = ToothGrowth)
anova(suppEff, MainEff)

## Analysis of Variance Table
## 
## Model 1: len ~ supp
## Model 2: len ~ supp + dose
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     58 3246.9                                  
## 2     56  820.4  2    2426.4 82.811 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Example

• We will extend the model by introducing an interaction between supp and dose.
• A graphical check for no interaction in the main effects model:

• Interaction plot:

with(ToothGrowth, interaction.plot(dose, supp, len, col = 2:3))

• For each of the supplement types we plot the average number of tooth length as a function of dosage.
• If the main effects model is correct then the difference between supplements is the same for all levels of dosage, i.e. the curves should be parallel - except for noise.
• This does not seem to be the case.

Dummy coding

• The extended model can be formulated as \[ \mathtt{length} = \mu+\beta_1 s_C+\beta_2 d_M+\beta_3 d_H+ \beta_4 s_C d_M+\beta_5 s_C d_H+\mathtt{error} \]
• Interpretation:
  • \(\mu\) is the expected tooth length for supp OJ and dose LO (\(s_C=d_M=d_H=0\)).
  • \(\beta_1\) is the effect of changing from supp OJ to VC, dose is LO (\(s_C=1,d_M=d_H=0\)).
  • \(\beta_2\) is the effect of increasing dose from LO to ME, when supp is OJ (\(s_C=0,d_M=1\)).
  • \(\beta_3\) is the effect of increasing dose from LO to HI, when supp is OJ (\(s_C=0,d_H=1\)).
  • \(\beta_4\) is an additional effect of both changing from supp OJ to VC and increasing dose from LO to ME (\(s_C=1,d_M=1\))
  • \(\beta_5\) is an additional effect of both changing from supp OJ to VC and increasing dose from LO to HI (\(s_C=1,d_H=1\))

Example

• We fit the interaction model by changing plus to multiply in the model expression from before:

Interaction <- lm(len ~ supp*dose, data = ToothGrowth)

• Now we can think of an experiment with 6 groups corresponding to each combination of the predictors.
• Looking at the group averages it looks like, the supplement types behave quite differently depending on dose:

mean(len ~ supp + dose, data = ToothGrowth)

## OJ.LO VC.LO OJ.ME VC.ME OJ.HI VC.HI 
## 13.23  7.98 22.70 16.77 26.06 26.14

• But is that significant?

anova(MainEff, Interaction)

## Analysis of Variance Table
## 
## Model 1: len ~ supp + dose
## Model 2: len ~ supp * dose
##   Res.Df    RSS Df Sum of Sq     F  Pr(>F)  
## 1     56 820.43                             
## 2     54 712.11  2    108.32 4.107 0.02186 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• With a p-value of 2.1860269% there is a significant interaction supp:dose, i.e. the lack of parallel curves in the interaction plot is significant.

summary(Interaction)

## 
## Call:
## lm(formula = len ~ supp * dose, data = ToothGrowth)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -8.20  -2.72  -0.27   2.65   8.27 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     13.230      1.148  11.521 3.60e-16 ***
## suppVC          -5.250      1.624  -3.233  0.00209 ** 
## doseME           9.470      1.624   5.831 3.18e-07 ***
## doseHI          12.830      1.624   7.900 1.43e-10 ***
## suppVC:doseME   -0.680      2.297  -0.296  0.76831    
## suppVC:doseHI    5.330      2.297   2.321  0.02411 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.631 on 54 degrees of freedom
## Multiple R-squared:  0.7937, Adjusted R-squared:  0.7746 
## F-statistic: 41.56 on 5 and 54 DF,  p-value: < 2.2e-16

• The additional effect of both changing from supp OJ to VC and increasing dose from LO to HI (\(\beta_5\)=suppVC:doseHI) is highly significant.