Adjusted Means

Having run an ANCOVA, and having estimated the relationship between IQ and memory score, we are now in a position to answer the interesting question, What would memory scores have been for placebo and drug groups IF their scores on IQ had been similar?. Of course their scores on IQ were not similar, but ANCOVA allows us to answer the question, IF they had been similar, then taking account what we know about how memory score varies with IQ (independent of any experimental manipulation, like placebo vs drug), what would their scores have been?

We are going to ask R to generate predicted values of memory score given the full model that includes both \(\alpha_{j}\) and \(\beta\) parameters. We will then ask R to generate predictions of what memory score would have been, for both placebo and drug groups, if their scores on IQ had been the same. We will have to choose some arbitrary value of IQ, the usual is to simply choose the grand mean of IQ across both groups. We will use the predict() command to generate predictions from our full model.

group.predict <- factor(c("placebo","drug"), levels=c("placebo","drug"))
iq.predict <- rep(mean(iq), 2)
data.predict <- data.frame(group=group.predict, iq=iq.predict)
print(data.predict)

##     group    iq
## 1 placebo 115.2
## 2    drug 115.2

adjmeans <- predict(mod.full, data.predict)
data.predict = data.frame(data.predict, memoryscore=adjmeans)
print(data.predict)

##     group    iq memoryscore
## 1 placebo 115.2    9.861678
## 2    drug 115.2   12.298322

So we see that our full model predicts that if the placebo and drug groups had been equal on IQ (for both groups, IQ = 115.2), then the mean memory score for the placebo group would have been 9.86 and for the drug group would have been 12.3. Compare that to the original unadjusted means (placebo = 8.36 and drug = 13.8) and you see that in this case, arithmetically adjusting for the group differences in IQ, resulted in less of a difference in memory score.

Graphically we can see how the prediction works, by drawing a vertical line at IQ = 115.2 and determining the value of memoryscore that intercepts that line for each group.

plot(iq,memoryscore, pch=c(rep(1,5),rep(24,5)), bg=c(rep(1,5),rep("red",5)))
legend("topleft", c("placebo","drug"), pch=c(1,24), pt.bg=c(1,"red"))
for (j in 1:2) {
    if (j==1) {
        w <- c(1,2,3,4,5)
        cl <- 1
    }
    if (j==2) {
        w <- c(6,7,8,9,10)
        cl <- "red"
    }
    lines(c(min(iq[w]),max(iq[w])),
        c(min(placebo.pred[w]), max(placebo.pred[w])), col=cl)
    abline(v=mean(iq), lty=3, col="blue")
    abline(h=adjmeans[1], lty=2, col="black")
    abline(h=adjmeans[2], lty=2, col="red")
}

By including a covariate, you can change your conclusions in many ways (or fail to change them at all). An apparent effect can become weaker, or reversed, or be nullified altogether. An effect can become much stronger. Think about the conditions that must exist in your dataset for these different kinds of scenarios to happen.

ANCOVA Assumptions

Scores on the dependent variable must be normally distributed, and they must be so at all values of the covariate.

The homogeneity of regression assumption is that the slope of the relationship between the dependent variable and the covariate is the same for all groups. In the equation for our full model this is apparent as \(\beta\) does not have a subscript (indicating different values for each group, i.e. each level of the factor).

Having said that, there is nothing stopping you from fitting a full model that has different slopes for each group. It would look like this:

\[ Y_{ijk} = \mu + \beta_{j} X_{ij} + \alpha_{j} + \epsilon_{ijk} \]

The cost is, you lose a degree of freedom for every extra slope you have to estimate. In addition however, it raises the troublesome question, why do different groups have different relationships between the dependent variable and the covariate? This may not be a problem. Then again it may be. It depends on the situation and the meaning of the dependent variable, the covariate, and the experimental design.

You can perform a test to see if the slopes are different very easily, by comparing the following two models:

\[ \begin{eqnarray} full: &Y_{ijk} &= \mu + \beta_{j} X_{ij} + \alpha_{j} + \epsilon_{ijk}\\ restricted: &Y_{ijk} &= \mu + \beta X_{ij} + \alpha_{j} + \epsilon_{ijk} \end{eqnarray} \]

In R you can simply compare a model that includes an interaction term with one that does not. In other words, in the full model, the effect of IQ on memory score (the slope) varies depending on the level of group (this is the definition of an interaction). In the restricted model, it does not.

mod.F <- lm(memoryscore ~ iq + group + iq:group)
mod.R <- lm(memoryscore ~ iq + group)
anova(mod.F, mod.R)

## Analysis of Variance Table
## 
## Model 1: memoryscore ~ iq + group + iq:group
## Model 2: memoryscore ~ iq + group
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1      6 19.519                           
## 2      7 20.913 -1   -1.3943 0.4286  0.537

Or we could simply list the anova table for the full model:

anova(mod.F)

## Analysis of Variance Table
## 
## Response: memoryscore
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## iq         1 115.743 115.743 35.5788 0.0009948 ***
## group      1  10.360  10.360  3.1848 0.1245828    
## iq:group   1   1.394   1.394  0.4286 0.5369540    
## Residuals  6  19.519   3.253                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

What is the Meaning of Adjusted Means?

Although ANCOVA allows you to arithmetically adjust scores to predict what they would be if groups were equal on a covariate, if the fact of the matter is, in your experiment, they are not equal on the covariate, then it may be difficult to interpret the adjusted scores, in particular when your experiment involves a treatment of some kind. For example, in the fake data we have presented here on memory, drug vs placebo, and IQ, what if the drug in question actually has different effects on memory in low IQ versus high IQ individuals? The only way to truly deal with this problem is to equate the subjects on the covariate in the first place, prior to the treatment. You also may want to consider turning your covariate into a real factor (e.g. crossing your design with levels of IQ, e.g. low, mid, high).

Followup tests on individual means

After a significant omnibus test on the ANCOVA, you may wish to perform follow-up tests on individual means. Tests are done on adjusted means. Unfortunately there is no button nor is there a simple R command to perform followup tests after an ANCOVA.

One way of doing followup tests, and in many ways the most flexible, is to compute an F-ratio as described in the Maxwell & Delaney text. In my hardcover 2nd edition — for pairwise comparisons, pg. 435, equations 37 and 38; and for any arbitrary contrasts: pg. 436, equations 39, 41 and 42. The denominator of the F-ratio is no longer just read off of the ANOVA output table as the mean-squared error, it’s mean-squared error multiplied by a rather ugly looking term that you need to compute. It wouldn’t be a big deal for you to write an R function to compute the denominator. Note that these are tests uncorrected for Type-I error. You can do Tukey corrections by converting the F-value you get to a q using \(q = \sqrt{2F}\), and then compute a probability based on the q distribution using \(df_{num}=1\) and \(df_{denom}=N-a-1\) where \(a\) is the number of groups.

Below is an example of how you would do this, for some made-up data. First, here are the fake data:

mydv <- c(2,3,5,4,5,6,6,8,9,7,13,15,16,17,18)
mygroup <- c(1,1,1,1,1,2,2,2,2,2,3,3,3,3,3)
mygroupf <- factor(mygroup)
mycovariate <- c(10,11,13,12,10,13,14,14,13,12,18,18,21,23,22)
rawmeans <- tapply(mydv, list(mygroupf), mean)
print(rawmeans)

##    1    2    3 
##  3.8  7.2 15.8

plot(mycovariate, mydv, pch=mygroup)
legend("topleft", legend=levels(mygroupf), pch=unique(mygroup))

Now let’s run the ANCOVA model.

mymod <- lm(mydv ~ mycovariate + mygroupf)
mymod.anova <- anova(mymod)
print(mymod.anova)

## Analysis of Variance Table
## 
## Response: mydv
##             Df Sum Sq Mean Sq  F value    Pr(>F)    
## mycovariate  1 379.14  379.14 239.5787 8.186e-09 ***
## mygroupf     2  14.39    7.19   4.5452   0.03641 *  
## Residuals   11  17.41    1.58                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Now let’s calculate the adjusted means.

# generate adjusted means
mygroup.predict <- factor(c(1,2,3))
mycovariate.predict <- rep(mean(mycovariate), 3)
mydata.predict <- data.frame(mygroupf=mygroup.predict,
        mycovariate=mycovariate.predict)
print(mydata.predict)

##   mygroupf mycovariate
## 1        1    14.93333
## 2        2    14.93333
## 3        3    14.93333

adjmeans <- predict(mymod, mydata.predict)
print(adjmeans)

##         1         2         3 
##  6.030303  8.235498 12.534199

Now let’s do a contrast comparing adjusted means 1 and 2.

# compare adjusted mean 1 and adjusted mean 2
Y1 <- adjmeans[1]   # adjusted mean 1
Y2 <- adjmeans[2]   # adjusted mean 2
a <- length(levels(mygroupf))   # total number of groups
N <- length(mydv)   # total number of subjects in whole experiment
n <- N/a # number of subjects in each group
mse <- mymod.anova$`Mean Sq`[3] # error term from the ANOVA output table
n1 <- n
n2 <- n
# the extra junk in equation 37, Maxwell & Delaney
i1 <- which(mygroup==1)
i2 <- which(mygroup==2)
i3 <- which(mygroup==3)
sse1 <- sum((mydv[i1]-mean(mydv[i1]))^2)
sse2 <- sum((mydv[i2]-mean(mydv[i2]))^2)
sse3 <- sum((mydv[i3]-mean(mydv[i3]))^2)
denom_mult <- (1/n1) + (1/n2) + (((Y1-Y2)^2)/(sse1+sse2+sse3))
s2y1y2 <- mse * denom_mult
dfnum <- 1
dfdenom <- N-a-1
Fobs <- ((Y1-Y2)^2)/s2y1y2
p.value.uncorrected = pf(Fobs, dfnum, dfdenom, lower.tail=FALSE)
print(p.value.uncorrected)

##          1 
## 0.04062419

qobs = sqrt(2*Fobs)
p.value.tukey = ptukey(q=qobs, nmeans=a, df=dfdenom, lower.tail=FALSE)
print(p.value.tukey)

##          1 
## 0.09477965

Another possibility is to use the multcomp package in R (which you will have to download and install). It can perform pairwise Tukey tests on adjusted means, in the following way. Note how the results are presented as a t-test, and the probability is not exactly the same as when using the Maxwell & Delaney equations. I don’t yet know why this is, I am investigating it. My guess is that the internals of glht() are doing something slightly different than what Maxwell & Delaney suggest.

library(multcomp)

## Loading required package: mvtnorm

## Loading required package: survival

## Loading required package: TH.data

## Loading required package: MASS

## 
## Attaching package: 'TH.data'

## The following object is masked from 'package:MASS':
## 
##     geyser

summary(glht(aov(mydv ~ mycovariate + mygroupf), linfct=mcp(mygroupf="Tukey")))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = mydv ~ mycovariate + mygroupf)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)  
## 2 - 1 == 0   2.2052     0.9157   2.408   0.0711 .
## 3 - 1 == 0   6.5039     2.2320   2.914   0.0296 *
## 3 - 2 == 0   4.2987     1.8157   2.368   0.0764 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)