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

[Robert Kushler]

"ST" is the ratio of "sigma-Rail" to "sigma".

Try the code below.  This approach handles simple models correctly,
but not more complex ones.  I'm sure I'll be scolded for accessing
the slots directly (and perhaps other infelicities), but others may
provide improved code.

Regards,   Rob Kushler


# code from original message:
library(lme4)
data(Rail, package="nlme")
fmRail.lme <- lmer(travel~1+(1|Rail), data=Rail)
summary(fmRail.lme)
Railmcmc <- mcmcsamp(fmRail.lme, n=10000, saveb=TRUE)
HPDinterval(Railmcmc)


# new stuff:
plot(t(Railmcmc at ST))  # inspect the chain (could be "stuck at zero")
sigmaRail <- as.vector(Railmcmc at ST*Railmcmc at sigma)
hist(sigmaRail)
plot(density(sigmaRail,adjust=2),main='Posterior distribution of "sigma-Rail"')
sort(sigmaRail)[c(50,950)]  # frequentist "equal tails" CI

# ad hoc HPD interval
hpdint <- function(x,prob=0.95) {
  xx <- sort(x)
  nn <- length(xx)
  len <- function(a)  xx[floor((a+prob)*nn)] - xx[floor(a*nn)+1]
  tail <- optim((1-prob)/2,len,method="L-BFGS-B",lower=0,upper=1-prob)$par
  c(xx[floor(tail*nn)+1],xx[floor((tail+prob)*nn)])
}
hpdint(sigmaRail,0.9)





jos matejus wrote:
> 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
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>