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

[jos matejus]

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)
$fixef
               lower    upper
(Intercept) 54.49488 78.99753
attr(,"Probability")
[1] 0.95

$ST
         lower    upper
[1,] 0.2582595 1.352534
attr(,"Probability")
[1] 0.95

$sigma
        lower    upper
[1,] 8.503474 23.07903
attr(,"Probability")
[1] 0.95

$ranef
           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
attr(,"Probability")
[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
Jos