[R-sig-ME] HPDinterval function with Rails example (P&B)

jos matejus matejus106 at googlemail.com
Wed Sep 10 16:31:29 CEST 2008

Dear lmers

I am trying to get to grips using the functions mcmcsamp() and
HPDintervals() to generate HPD confidence intervals for my parameter
estimates for some simple mixed effects models. In order to gently
introduce myself to this subject I have been working through examples
given in Pinheiro and Bates- 'mixed-effects models in S and S-Plus
using lme4. I have a couple of questions that I have been unable to
find the answers to.

The example I am following at the moment is using the Rails dataset
(page 9, P&B). I have fitted the following model:

> library(lme4)
> data(Rail, package="nlme")
> fmRail.lme <- lmer(travel~1+(1|Rail), data=Rail)
> summary(fmRail.lme)
Linear mixed model fit by REML

Formula: travel ~ 1 + (1 | Rail)
   Data: Rail
   AIC   BIC logLik deviance REMLdev
 128.2 130.8 -61.09    128.6   122.2
Random effects:
 Groups   Name        Variance Std.Dev.
 Rail     (Intercept) 615.266  24.8046
 Residual              16.167   4.0208
Number of obs: 18, groups: Rail, 6

Fixed effects:
            Estimate Std. Error t value
(Intercept)    66.50      10.17   6.538

In order to obtain HPD intervals for the estimated parameters

> Railmcmc <- mcmcsamp(fmRail.lme, n=10000, saveb=TRUE)
> HPDinterval(Railmcmc)
               lower    upper
(Intercept) 54.49488 78.99753
[1] 0.95

         lower    upper
[1,] 0.2582595 1.352534
[1] 0.95

        lower    upper
[1,] 8.503474 23.07903
[1] 0.95

           lower       upper
[1,] -38.8745199  0.09389242
[2,] -26.8906507  5.99597084
[3,] -24.6996428  7.63476719
[4,]  -5.3189656 27.55519181
[5,]  -4.7411818 28.25215767
[6,]  -0.3872951 35.66205005
[1] 0.95

  I understand the HPD intervals for the fixed effects and for sigma
(although they are quite a bit different than those obtained using the
intervals() function in nlme), however I am unclear as to what the
intervals in the  ST represent. Also, is there a way of calculating
(or accessing)  the HPD interval for the random effect associated with
rail (sigma b ) as I cant seem to see it (I remember you got it using
the function HPDinterval in the coda package). Apologies in advance if
this is a rather trivial question, but its just got me stumped and I
would really like to get to grips with this.

Best wishes

More information about the R-sig-mixed-models mailing list