QK Box 8.15

We use a simple hypothetical data set to illustrate some issues when covariate ranges do not match well. The data set consists of two treatment groups, a single response variable (Y) and a single covariate (X) Thie data file is here

Preliminaries

First, load the required packages (car)

Import covrange data file

covrange <- read.csv("../data/covrange.csv")
head(covrange,10)

#Make group a factor
covrange$group<-factor(covrange$group)

Fit full model with interaction to evaluate homogeneous slopes

All SS types produce same result for interaction

covrange.aov1 <- aov(y~group+x+group*x, data = covrange)
plot(covrange.aov1)

Anova(covrange.aov1, type = 3)

Anova Table (Type III tests)

Response: y
            Sum Sq Df F value    Pr(>F)    
(Intercept) 35.700  1 40.8864 8.885e-06 ***
group        5.016  1  5.7442   0.02911 *  
x            7.191  1  8.2357   0.01112 *  
group:x      2.445  1  2.8005   0.11367    
Residuals   13.970 16                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Not strong evidence for heterogenous slopes

Run simpler model

covrange.aov2 <- aov(y~group+x, data = covrange)
Anova(covrange.aov2, type = 3)

Anova Table (Type III tests)

Response: y
            Sum Sq Df  F value    Pr(>F)    
(Intercept) 99.269  1 102.8038 1.264e-08 ***
group        4.681  1   4.8472   0.04179 *  
x            5.157  1   5.3405   0.03364 *  
Residuals   16.416 17                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get adjusted means

adjmeans <- effect("group", covrange.aov2, se = TRUE)
summary(adjmeans)

 group effect
group
        1         2 
 8.951958 10.367042 

 Lower 95 Percent Confidence Limits
group
       1        2 
8.130592 9.545677 

 Upper 95 Percent Confidence Limits
group
        1         2 
 9.773323 11.188408 

Use Huitema adjustment

Instead of using raw x values as covariates, Huitema suggests using the x-residuals from each group. This is done already in the data file (xres), but it can be done by fitting a lm with x as the response and group as the predictor, and using the calculated residuals. We’ll skip the comparison of slopes and just run the pooled slopes model and get adjusted means

covrange.aov3 <- aov(y~group+xres, data = covrange)
Anova(covrange.aov3, type = 3)

Anova Table (Type III tests)

Response: y
            Sum Sq Df  F value    Pr(>F)    
(Intercept) 707.28  1 732.4627 2.033e-15 ***
group        31.23  1  32.3367 2.677e-05 ***
xres          5.16  1   5.3405   0.03364 *  
Residuals    16.42 17                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

adjmeans <- effect("group", covrange.aov3, se = TRUE)
summary(adjmeans)

 group effect
group
     1      2 
 8.410 10.909 

 Lower 95 Percent Confidence Limits
group
        1         2 
 7.754387 10.253387 

 Upper 95 Percent Confidence Limits
group
        1         2 
 9.065613 11.564613 

Generate graph

df1<-covrange  # Copy dataframe to df1 so generic code can be used
aggregate(x = df1$x,        # Specify data column
     by = list(df1$group),       # Specify group indicator
     FUN = mean) 

p1<-ggplot(data = df1,aes(x = x,y = y,group = group, shape = group, color = group))+
 geom_point()+
 geom_smooth(method = "lm",se = FALSE)+
 geom_vline(xintercept = 3.872, linetype = "dashed", color = lf)+
 geom_segment(x = 2.803,y = 6,xend = 2.803,yend = 7, arrow = arrow(length = unit(3, "mm")),color = lc)+
 geom_segment(x = 4.94,y = 6,xend = 4.94,yend = 7, arrow = arrow(length = unit(3, "mm")),color = lc)+
 labs(x = "X", y = NULL)+
 scale_color_grey(end = 0.6)+
 theme_qk()+
 theme(
  legend.position = "none"
 )
p2<-ggplot(data = df1,aes(x = xres,y = y,group = group, shape = group, color = group))+
 geom_point()+
 geom_smooth(method = "lm",se = FALSE)+
 geom_vline(xintercept = 0, linetype = "dashed", color = lf)+
 labs(x = "X residual", y = NULL)+
 scale_color_grey(end = 0.6)+
 theme_qk()+
 theme(
  legend.position = "none"
 )
p3<-ggplot(data = df1,aes(x = x,y = y))+
 geom_point(color = sc)+
 geom_smooth(method = "lm",se = FALSE, color = sc,linetype = "dashed")+
 geom_smooth(method = "loess", color = lc, linetype = "solid", se = FALSE, span = 2)+
 labs(x = "X", y = NULL)+
 scale_color_grey(end = 0.6)+
 theme_qk()+
 theme(
  legend.position = "none"
 )
p1+p2+p3