Analysis of covariance
Linear models with at least one covariate (quantitative regressor) and one factor (qualitative regressor) are historically called Analysis of Covariance models. The ChickWeight dataset contains the factor Diet and covariate Time.
Different models are possible:
Model 1: common slope, common intercept: \(weight_{ij} = \beta_0 + \beta_1 Time_{ij} + \epsilon_{ij}\)
Model 2: common slope, different intercepts: \(weight_{ij} = \beta_{0i} + \beta_1 Time_{ij} + \epsilon_{ij}\)
Model 3: different slopes, different intercepts: \(weight_{ij} = \beta_{0i} + \beta_{1i} Time_{ii} + \epsilon_{ij}\)
Model 4: different slope, common intercept: \(weight_{ij} = \beta_{0} + \beta_{1i} Time_{ii} + \epsilon_{ij}\)
In each model we assume independent errors \(\epsilon_{ij} \sim \mathbf{N}(0,\sigma^2)\).
Graphical display of the models
To display the models, we first store the intercepts and slopes and vectors a and b. We put these vectors in a new dataframe called df. For the first model, we have a single line, for the second model we have four different lines. Hence, the data.frame will have 1 + 3*4 = 13 rows.
df <- data.frame(r = 1:13,a = NA,b = NA,Diet = as.character(c(0,rep(1:4,3))),
model = c('ancW1',rep(c('ancW2','ancW3','ancW4'),each = 4)))
df$a[1] <- coef(ancW1)[1] # common intercept
df$b[1] <- coef(ancW1)[2] # common slope
df$a[2:5] <- coef(ancW2)[1:4] # 4 intercepts in 4 groups
df$b[2:5] <- rep(coef(ancW2)[5],4) # common slope in 4 groups
df$a[6:9] <- coef(ancW3)[1:4] + c(0,rep(coef(ancW3)[1],3)) # 4 intercepts in 4 groups
df$b[6:9] <- coef(ancW3)[5:8] + c(0,rep(coef(ancW3)[5],3)) # 4 slopes in 4 groups
df$a[10:13] <- rep(coef(ancW4)[1],4) # common intercept in 4 groups
df$b[10:13] <- coef(ancW4)[2:5] # 4 slopes in 4 groups
CWd <- data.frame(ChickWeight,ChickWeight,ChickWeight,ChickWeight,
model = rep(c('ancW1','ancW2','ancW3','ancW4'), each = nrow(ChickWeight)))
## Warning in data.frame(ChickWeight, ChickWeight, ChickWeight, ChickWeight, :
## row names were found from a short variable and have been discarded
CWd$Diet <- as.character(CWd$Diet)
ggplot(data = CWd,aes(x = Time,y = weight,colour = Diet,shape = Diet)) +
geom_point() + facet_wrap(~model) +
geom_abline(data = df,size = 1, aes(intercept = a,
slope = b,colour = Diet,shape = Diet))
Testing models
anova(ancW1,ancW2,ancW3,ancW4)
## Analysis of Variance Table
##
## Model 1: weight ~ Time
## Model 2: weight ~ Diet + Time
## Model 3: weight ~ Diet + Time + Diet:Time
## Model 4: weight ~ Diet:Time
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 576 872212
## 2 573 742336 3 129876 37.3020 < 2.2e-16 ***
## 3 570 661532 3 80804 23.2079 3.474e-14 ***
## 4 573 665525 -3 -3993 1.1469 0.3295
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Model ancW2 fits significantly better than model ancW1. Model ancW3 fits significantly better than model ancW2. Model ancW3 does not significantly fit better than model ancW4, though.
Finding the best model
We cannot formally test whether models ancW2 and ancW4 differ significantly, because they are not nested: neither is model ancW2 a special case of ancW4, nor vice versa. A less formal comparison of the models can be made with the anova function: we see that the residual sum of squares of model ancW4 is 76811 lower than that of model ancW2. So, model ancW4 fits better.
## Analysis of Variance Table
##
## Model 1: weight ~ Diet + Time
## Model 2: weight ~ Diet:Time
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 573 742336
## 2 573 665525 0 76811
Reproduce the F-statistic for comparing larger model 3 with smaller model 4
Test statistic F is a corrected difference in RSS:
\(F = \frac{\frac{RSS_4 - RSS_3}{p_3 - p_4}}{\frac{RSS_3}{n - p_3}}\)
\(F \sim F_{p_4-p_3,n -p_3}\)
RSS4 <- deviance(ancW4)
RSS3 <- deviance(ancW3)
p4 <- ancW4$rank
p3 <- ancW3$rank
df3 <- ancW3$df.residual
(F <- ((RSS4 - RSS3)/(p3 - p4))/ (RSS3/df3))
## [1] 1.146906
(p <- 1-pf(F,df1 = p3 - p4, df2 = df3))
## [1] 0.3295058
## [1] FALSE