mcmcsamp(lme4): What is contained in $ST and $sigma?

Dear List,

I am new to lme4. I believe it will be very useful for analysis of
pharmaceutical stability data in which batches are considered to have
random slopes and intercepts. We now use SAS Proc Mixed for this, but I
like the lme4 approach and the other capabilities of R better.

What attracts me to lme4 is the possibility of obtaining a posterior
sample using mcmcsamp(). How else can one obtain meaningful interval
estimates or make predictions about future data?

However, I do not know how to interpret the contents of the mcmcsamp
output slots $ST and $sigma. I find that the samples contained in those
slots do not agree with the estimates provided by lmer(). I hope someone
more knowledgeable will be willing to explain this to me.

To make things concrete, I provide below an example analysis of simulated
data.

The following code generates 50 batches worth of stability tests. Batch
intercepts and slopes are drawn from N(100,,4^2) and N(-1,2^2). Analytical
results at each time point (0,3,...,36 months) are drawn from
N(intercept+slope*month, 1^2). So there are 3 random paramters (variance
for slope, intercept, and error) and 2 fixed effects (mean slope and
intercept). With such a large balanced data set, I expect rapid
convergence of the Gibbs sampler.

B<-50 # number of batches
Months<-c(0,3,6,9,12,24,30,36) # pull points
n.Months<-length(Months)
Init.mean<-100 # initial level for the process
Init.SD<-4 # SD of the process initial
Slope.mean<- -1
Slope.SD<- 2
Error.SD<- 1
# Prepare Data Vectors
Batch<-rep(NA,B*n.Months)
Month<-rep(NA,B*n.Months)
Y<-rep(NA,B*n.Months)
for(i in 1:B){
 slope<-rnorm(1,Slope.mean,Slope.SD)
 intercept<-rnorm(1,Init.mean,Init.SD)
 for(j in 1:n.Months){
   k<-(i-1)*n.Months+j
   Batch[k]<-i
   Month[k]<-Months[j]
   Y[k]<-intercept+Month[k]*slope+rnorm(1,0,Error.SD)
 }
}
Stab<-cbind(Batch=factor(Batch),Month=Month,Y=Y);Stab

I use the following to analyze these data:

library(lme4)
anal<-lmer(Y~1+Month+(1|Batch)+(-1+Month|Batch))
anal

The following output is obtained:

Linear mixed model fit by REML
Formula: Y ~ 1 + Month + (1 | Batch) + (-1 + Month | Batch)
 AIC  BIC logLik deviance REMLdev
 1821 1841 -905.4     1811    1811
Random effects:
 Groups   Name        Variance Std.Dev.
 Batch    (Intercept) 19.5551  4.42212
 Batch    Month        4.4818  2.11704
 Residual              0.9780  0.98894
Number of obs: 400, groups: Batch, 50

Fixed effects:
           Estimate Std. Error t value
(Intercept) 100.0254     0.6302  158.73
Month        -1.3988     0.2994   -4.67

Correlation of Fixed Effects:
     (Intr)
Month -0.001

This is terrific. The estimates agree well with theory. Armed with
confidence (or rather credibility?) I obtain interval estimate as follows:

HPDinterval(mcmcsamp(anal,n=10000))

The following output is obtained:

$fixef
               lower       upper
(Intercept) 99.447112 100.5605457
Month       -1.905921  -0.8496636
attr(,"Probability")
[1] 0.95

$ST
        lower     upper
[1,] 0.7546098 0.9719147
[2,] 0.6913345 1.7690498
attr(,"Probability")
[1] 0.95

$sigma
       lower    upper
[1,] 1.710004 2.169471
attr(,"Probability")
[1] 0.95

While the fixed effect posteriors compare very well with truth and with
lmer point estimates, the $sigma is strange. I expect it to be close to 1
since this is the error SD used to generate the data and also the value
reported by lmer. Instead $sigma contains a values between 1.7 and 2.2.
Why?

mcmcsamp(lme4): What is contained in $ST and $sigma?

Thread (4 messages)