2

Does anyone know how to do post hoc tests in an ANCOVA model with a factorial design?

I have two vectors consisting of 23 baseline values (covariate) and 23 values after treatment (dependent variable) and I have two factors with both two levels. I created an ANCOVA model and calculated the adjusted means, standard errors and confidence intervals. Example:

library(effects)

baseline = c(0.7672,1.846,0.6487,0.4517,0.5599,0.2255,0.5946,1.435,0.5374,0.4901,1.258,0.5445,1.078,1.142,0.5,1.044,0.7824,1.059,0.6802,0.8003,0.5547,1.003,0.9213)
after_treatment = c(0.4222,1.442,0.8436,0.5544,0.8818,0.08789,0.6291,1.23,0.4093,0.7828,-0.04061,0.8686,0.8525,0.8036,0.3758,0.8531,0.2897,0.8127,1.213,0.05276,0.7364,1.001,0.8974)

age = factor(c(rep(c("Young","Old"),11),"Young")) 
treatment = factor(c(rep("Drug",12),rep("Placebo",11)))

ANC = aov(after_treatment ~ baseline + treatment*age)

effect_treatage = effect("treatment*age",ANC)
data.frame(effect_treatage)

  treatment   age       fit        se     lower     upper
1      Drug   Old 0.8232137 0.1455190 0.5174897 1.1289377
2   Placebo   Old 0.6168641 0.1643178 0.2716452 0.9620831
3      Drug Young 0.5689036 0.1469175 0.2602413 0.8775659
4   Placebo Young 0.7603360 0.1462715 0.4530309 1.0676410

Now I want test if there is a difference between the adjusted means of

Young-Placebo:Young-Drug
Old-Placebo:Old-Drug
Young-Placebo:Old-Drug
Old-Placebo:Young-Drug

So I tried:

library(multcomp)
pH = glht(ANC, linfct = mcp(treatment*age="Tukey"))
# Error: unexpected '=' in "ph = glht(ANC_nback, linfct = mcp(treat*age="

And:

pH = TukeyHSD(ANC)
# Error in rep.int(n, length(means)) : unimplemented type 'NULL' in 'rep3'
# In addition: Warning message:
# In replications(paste("~", xx), data = mf) : non-factors ignored: baseline

Does anyone know how to resolve this?

Many thanks!

PS for more info see

R: How to graphically plot adjusted means, SE, CI ANCOVA

Community
  • 1
  • 1
Alexm
  • 57
  • 2
  • 6

3 Answers3

3

If you wish to use multcomp, then you can take advantage of the wonderful and seamless interface between lsmeans and multcomp packages (see ?lsm), whereas lsmeans provides support for glht(). 

baseline = c(0.7672,1.846,0.6487,0.4517,0.5599,0.2255,0.5946,1.435,0.5374,0.4901,1.258,0.5445,1.078,1.142,0.5,1.044,0.7824,1.059,0.6802,0.8003,0.5547,1.003,0.9213)
after_treatment = c(0.4222,1.442,0.8436,0.5544,0.8818,0.08789,0.6291,1.23,0.4093,0.7828,-0.04061,0.8686,0.8525,0.8036,0.3758,0.8531,0.2897,0.8127,1.213,0.05276,0.7364,1.001,0.8974)
age = factor(c(rep(c("Young","Old"),11),"Young")) 
treatment = factor(c(rep("Drug",12),rep("Placebo",11)))
Treat <- data.frame(baseline, after_treatment, age, treatment)

ANC  <-  aov(after_treatment ~ baseline + treatment*age, data=Treat)

library(multcomp)
library(lsmeans)

summary(glht(ANC, linfct = lsm(pairwise ~ treatment * age)))
## Note: df set to 18
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: aov(formula = after_treatment ~ baseline + treatment * age, data = Treat)
## 
## Linear Hypotheses:
##                                  Estimate Std. Error t value Pr(>|t|)
## Drug,Old - Placebo,Old == 0       0.20635    0.21913   0.942    0.783
## Drug,Old - Drug,Young == 0        0.25431    0.20698   1.229    0.617
## Drug,Old - Placebo,Young == 0     0.06288    0.20647   0.305    0.990
## Placebo,Old - Drug,Young == 0     0.04796    0.22407   0.214    0.996
## Placebo,Old - Placebo,Young == 0 -0.14347    0.22269  -0.644    0.916
## Drug,Young - Placebo,Young == 0  -0.19143    0.20585  -0.930    0.789
## (Adjusted p values reported -- single-step method)

This eliminates the need for reparametrization. You can achieve the same results by using lsmeans alone: 

lsmeans(ANC, list(pairwise ~ treatment * age))
## $`lsmeans of treatment, age`
##  treatment age      lsmean        SE df  lower.CL  upper.CL
##  Drug      Old   0.8232137 0.1455190 18 0.5174897 1.1289377
##  Placebo   Old   0.6168641 0.1643178 18 0.2716452 0.9620831
##  Drug      Young 0.5689036 0.1469175 18 0.2602413 0.8775659
##  Placebo   Young 0.7603360 0.1462715 18 0.4530309 1.0676410
## 
## Confidence level used: 0.95 
## 
## $`pairwise differences of contrast`
##  contrast                       estimate        SE df t.ratio p.value
##  Drug,Old - Placebo,Old       0.20634956 0.2191261 18   0.942  0.7831
##  Drug,Old - Drug,Young        0.25431011 0.2069829 18   1.229  0.6175
##  Drug,Old - Placebo,Young     0.06287773 0.2064728 18   0.305  0.9899
##  Placebo,Old - Drug,Young     0.04796056 0.2240713 18   0.214  0.9964
##  Placebo,Old - Placebo,Young -0.14347183 0.2226876 18  -0.644  0.9162
##  Drug,Young - Placebo,Young  -0.19143238 0.2058455 18  -0.930  0.7893
## 
## P value adjustment: tukey method for comparing a family of 4 estimates
landroni
  • 2,902
  • 1
  • 32
  • 39
  • Is is possible to extract connecting letters (e.g., Tukey style) using this approach? – DirtStats Nov 01 '16 at 13:31
  • Ah, figured it out. For Tukey connecting letters you need to load the `multcompView` package. Like this: `library(multcompView); cld(lsmeans(ANC, list(pairwise ~ treatment * age)), alpha=0.05)` – DirtStats Nov 01 '16 at 14:19
2

Reparametrization is a possibility here:

treatAge <- interaction(treatment, age)
ANC1 <- aov(after_treatment ~ baseline + treatAge)
#fits are equivalent:
all.equal(logLik(ANC), logLik(ANC1))
#[1] TRUE

library(multcomp)
summary(glht(ANC1, linfct = mcp(treatAge="Tukey")))

#    Simultaneous Tests for General Linear Hypotheses
#
#Multiple Comparisons of Means: Tukey Contrasts
#
#
#Fit: aov(formula = after_treatment ~ baseline + treatAge)
#
#Linear Hypotheses:
#                                 Estimate Std. Error t value Pr(>|t|)
#Placebo.Old - Drug.Old == 0      -0.20635    0.21913  -0.942    0.783
#Drug.Young - Drug.Old == 0       -0.25431    0.20698  -1.229    0.617
#Placebo.Young - Drug.Old == 0    -0.06288    0.20647  -0.305    0.990
#Drug.Young - Placebo.Old == 0    -0.04796    0.22407  -0.214    0.996
#Placebo.Young - Placebo.Old == 0  0.14347    0.22269   0.644    0.916
#Placebo.Young - Drug.Young == 0   0.19143    0.20585   0.930    0.789
#(Adjusted p values reported -- single-step method)
Roland
  • 127,288
  • 10
  • 191
  • 288
  • There is a way to achieve the same results without reparametrization. See this [other answer](http://stackoverflow.com/questions/23628323/error-with-post-hoc-tests-in-ancova-with-factorial-design/31302451#31302451). – landroni Jul 08 '15 at 19:56
2

You need to use the which argument in TukeyHSD; "listing terms in the fitted model for which the intervals should be calculated". This is needed because you have a non-factor variable in the model ('baseline'). The variable causes trouble when included, which is default when which is not specified.

ANC = aov(after_treatment ~ baseline + treatment*age)
TukeyHSD(ANC, which = c("treatment:age"))

If you wish to use the more flexible glht, see section 3, page 8- here

Henrik
  • 65,555
  • 14
  • 143
  • 159