[R-sig-ME] confidence intervals for random slopes

[Ben Bolker]

On 14-08-10 10:14 AM, Bokony Veronika wrote:
> Dear all,
> 
> let me ask your advice about calculating confidence intervals for
> random slopes. I have a random intercept & random slope model like
> this:
> 
> lme(Y ~ X * fixfactor, random=~X|randomfactor)
> 
> where randomfactor has 9 levels. For each of these 9 levels, I can
> calculate the slope of X from the fixed and random effect estimates.
> I would like to add some measure of uncertainty to each of these 9
> estimates, i.e. an SE or CI for each random slope. The random effects
> SD which is given by the lme summary output is not what I'm
> interested in. I got a very general tip that I could calculate
> credible intervals from posterior distributions using a bayesian
> approach, but I found no evident way of extracting these from
> MCMCglmm either.
> 
> I would be very grateful for any working example that can achieve
> this.


library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
rslopes <- coef(fm1)$Subject[,"Days"]
varfix <- vcov(fm1)["Days","Days"]
re <- ranef(fm1,condVar=TRUE)
varcm <- attr(re$Subject,"postVar")[2,2,]
vartot <- varfix+varcm
library(plotrix)
plotCI(1:length(rslopes),rslopes,2*sqrt(vartot))

  Note that this approach ignores correlation between the conditional
mean/BLUP and the fixed-effect estimate, a topic which has been
previously debated on this list ...

Now with MCMCglmm, for comparison:

library(MCMCglmm)
## turns out we have to specify a slightly more informative prior ...
priorb <- list(R = list(V = diag(1), nu = 0.002),
      G = list(G1 = list(V = diag(2), nu = 0.002)))
fm2 <- MCMCglmm(Reaction~Days,random=~us(1+Days):Subject,
                data=sleepstudy,prior=priorb, pr=TRUE, verbose=FALSE)

fslope <- fm2$Sol[,"Days"]  ## fixed slope
## slope random effects
rslopes2 <- fm2$Sol[,grep("Subject\\.Days\\.Subject",colnames(fm2$Sol))]
dim(rslopes2)
allslopes <- sweep(rslopes2,1,fslope,"+")  ## fixed + random
allslopemean <- colMeans(allslopes)
allsd <- apply(allslopes,2,sd)

Compare:

plotCI(1:length(rslopes),rslopes,2*sqrt(vartot))
plotCI(1:length(rslopes),allslopemean,2*allsd,col=2,add=TRUE)