[R-sig-ME] Bug in mcmcsamp

[Douglas Bates]

On 8/12/07, Jarrod Hadfield <j.hadfield at ed.ac.uk> wrote:
> Hi,

> I'm using lme4 0.99875-6 and I think there may be a bug in mcmcsamp.
> When I fit 2 random effects, the posterior variances always seem to
> be approximately the same irrespective of the true variances (or the
> REML estimates of the true variances).  For example:

>     A<-gl(100,10)
>     B<-gl(100,1,1000)
>     y<-rnorm(100,0,sqrt(1))[A]+rnorm(100,0,sqrt(3))[B]+rnorm
> (1000,0,sqrt(3))
>     model1<-lmer(y~(1|A)+(1|B))
>     model1MCMC<-mcmcsamp(model1,n=5000,trans=FALSE)
>     plot(model1MCMC)

Thanks for reporting the problem, Jarrod.  I tried your script with
the development version of lme4 and got the expected results, which I
enclose Could you try the mcmcsamp.R script in lme4 0.99875-6 and tell
me if the results are similar?
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Hadfield.pdf
Type: application/pdf
Size: 338640 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20070812/dd518120/attachment.pdf>
-------------- next part --------------
## bug report by Jarrod Hadfield (2007-08-12)
library(lme4)
library(coda)
set.seed(1)
A<-gl(100,10)
B<-gl(100,1,1000)
y<-rnorm(100,0,sqrt(1))[A]+rnorm(100,0,sqrt(3))[B]+rnorm(1000,0,sqrt(3))
model1<-lmer(y~(1|A)+(1|B))
model1MCMC<-mcmcsamp(model1,n=5000,trans=FALSE)
head(as.data.frame(model1MCMC), n = 20)
HPDinterval(model1MCMC)
pdf(file = "Hadfield.pdf", height = 8, width = 6)
xyplot(model1MCMC, strip = FALSE, strip.left = TRUE,
       scales = list(x=list(axs='i')), type = c("g", "l"))
densityplot(model1MCMC, plot.points = FALSE)
dev.off()