Skip to content
Prev 4799 / 7420 Next

ANCOVA with random effects for slope and intercept

peterhouk1 . <peterhouk at ...> writes:
There are a few things going on here.

* Plotting your data (see below), it looks pretty clear that
you've already mean-corrected it (all groups have mean value
of zero at the center point of the data) -- you probably want
to take this into account in your model.

* Even with this correction, you still get an estimate of zero
for the variance in slopes among groups.  This is not terribly
surprising for small-to-moderate sized, fairly
noisy (large residual variance) data sets.  This does indeed
mean that you will get essentially the same answer from a
plain old linear model.

* I did the stuff below with lmer, but I'd expect similar answers
from lme  . I used lmerTest for easy access to denominator df
approximates.

* The overall effect is not super-strong (R^2 ~ 0.12), but
very clearly statistically significant (depending on how you
do the calculation, p-value ranges from 0.001 to 0.009 ...)

## get data, draw pictures
mydata <- read.table(header=TRUE,"houk.txt")
library("ggplot2"); theme_set(theme_bw())

## make grouping variable into a factor; not totally necessary
## but doesn't hurt/probably a good idea
mydata <- transform(mydata,group_factor=factor(group_factor))
meanpred <- mean(mydata$predictor)
ggplot(mydata,aes(predictor,response,colour=group_factor))+geom_point()+
    geom_smooth(method="lm",se=FALSE)+
        geom_vline(xintercept=meanpred)

## centered version of predictor
mydata <- transform(mydata,cpred=predictor-meanpred)
library("lme4")

## original random-slopes model
(m1 <- lmer(response~predictor+(predictor|group_factor), data=mydata))

## model with centered predictor; suppress among-group variation in
## intercept
(m2 <- lmer(response~cpred+(cpred-1|group_factor), data=mydata))

summary(m2B <- lm(response~cpred, data=mydata))

## get p-value, using common-sense/classical value of 10 for denominator df
(cc1 <- coef(summary(m2)))
tval <- cc1["cpred","t value"]
pt(abs(tval),df=10,lower.tail=FALSE)

## refit model with lmerTest for convenient access to df estimates
detach("package:lme4")
library("lmerTest")
(m3 <- lmer(response~cpred+(cpred-1|group_factor), data=mydata))
anova(m3,ddf="Kenward-Roger")  ## gives denDF=8.67, p-value=0.0097
anova(m3,ddf="Satterthwaite")  ## gives denDF=69.9 (?), p-value=0.001
## essentially equivalent to lm() because denDF are estimated to be large

detach("package:lmerTest")