[R] Constructing predictions from HPDinterval() after lmer()

[Michael Kubovy]

Dear r-helpers,

Following up on http://finzi.psych.upenn.edu/R/Rhelp02a/archive/ 
81159.html where Douglas Bates gives a helpful application of lmer()  
to data(sleepstudy, package = 'lme4'), I need a bit more help in  
order to plot the correct confidence intervals of a designed  
experiment such as:

 > data(ratdrink, package = 'faraway')

I follow the steps Douglas took in his example:

 > summary( rd.er <- lmer(wt ~ weeks*treat + (1 | subject), data =  
ratdrink) )

Linear mixed-effects model fit by REML
Formula: wt ~ weeks * treat + (1 | subject)
    Data: ratdrink
    AIC   BIC logLik MLdeviance REMLdeviance
962.4 982.8 -474.2      964.3        948.4
Random effects:
Groups   Name        Variance Std.Dev.
subject  (Intercept) 71.206   8.4384
Residual             51.222   7.1570
number of obs: 135, groups: subject, 27

Fixed effects:
                       Estimate Std. Error t value
(Intercept)            52.8800     3.1928   16.56
weeks                  26.4800     0.7157   37.00
treatthiouracil         4.7800     4.5153    1.06
treatthyroxine         -0.7943     4.9756   -0.16
weeks:treatthiouracil  -9.3700     1.0121   -9.26
weeks:treatthyroxine    0.6629     1.1153    0.59

Correlation of Fixed Effects:
             (Intr) weeks  trtthr trtthy wks:trtthr
weeks       -0.448
treatthircl -0.707  0.317
treatthyrxn -0.642  0.288  0.454
wks:trtthrc  0.317 -0.707 -0.448 -0.203
wks:trtthyr  0.288 -0.642 -0.203 -0.448  0.454

 > rd.mc <- mcmcsamp(rd.er, 50000)
 > library(coda)
 > HPDinterval(rd.mc)

                            lower      upper
(Intercept)            46.420404  59.406398
weeks                  25.070131  27.930363
treatthiouracil        -4.420942  14.009291
treatthyroxine        -10.758369   9.435761
weeks:treatthiouracil -11.404620  -7.337025
weeks:treatthyroxine   -1.603858   2.842704
log(sigma^2)            3.683085   4.226153
log(sbjc.(In))          3.633853   4.955867
deviance              965.750351 980.825613
attr(,"Probability")
[1] 0.95


*******************To make sure my request is clear*******************
What I need is the analog of what is produced by all.effects() in the  
effects package:

 > summary(rd <- lm(wt ~ weeks*treat, data = ratdrink))

Call:
lm(formula = wt ~ weeks * treat, data = ratdrink)

Residuals:
     Min      1Q  Median      3Q     Max
-23.514  -6.660   0.230   6.914  28.343

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)
(Intercept)            52.8800     2.6547  19.919  < 2e-16 ***
weeks                  26.4800     1.0838  24.433  < 2e-16 ***
treatthiouracil         4.7800     3.7544   1.273    0.205
treatthyroxine         -0.7943     4.1371  -0.192    0.848
weeks:treatthiouracil  -9.3700     1.5327  -6.113 1.08e-08 ***
weeks:treatthyroxine    0.6629     1.6890   0.392    0.695
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.84 on 129 degrees of freedom
Multiple R-Squared: 0.9121,	Adjusted R-squared: 0.9087
F-statistic: 267.8 on 5 and 129 DF,  p-value: < 2.2e-16

 > rd.eff <- all.effects(rd)
 > rd.ci <- data.frame(y = rd.eff[[1]]$fit, Lower = rd.eff[[1]] 
$lower, Upper = rd.eff[[1]]$upper)

_____________________________
Professor Michael Kubovy
University of Virginia
Department of Psychology
USPS:     P.O.Box 400400    Charlottesville, VA 22904-4400
Parcels:    Room 102        Gilmer Hall
         McCormick Road    Charlottesville, VA 22903
Office:    B011    +1-434-982-4729
Lab:        B019    +1-434-982-4751
Fax:        +1-434-982-4766
WWW:    http://www.people.virginia.edu/~mk9y/