[R-sig-ME] Discrepancy between mcmcsamp value and variance component esitamate
Alan Bergland
Alan_Bergland at brown.edu
Thu Jan 31 21:28:15 CET 2008
Hi all,
I am getting a funny result using the mcmcsamp function, perhaps it
is a bug or perhaps there are convergence errors. I get this result
using the latest stable version of lme4 on CRAN and the alpha version
on R-forge, on both a Mac and linux machines. I apologize for not
being able to construct a self contained example, so perhaps my
output will suffice.
My data is structured in the following way:
> str(ov)
'data.frame': 7958 obs. of 7 variables:
$ ovn : num 27 16 29 22 13 18 17 12 16 18 ...
$ geno : Factor w/ 214 levels "2","3","10","12",..: 52 99 25 25 41
32 24 32 99 207 ...
$ food : Factor w/ 4 levels "0.2","0.4","0.6",..: 3 3 3 3 3 3 3 3 3
3 ...
$ block: Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
$ gf : Factor w/ 856 levels "2:0.8","2:0.6",..: 206 394 98 98 162
126 94 126 394 826 ...
$ gb : Factor w/ 642 levels "2:1","2:2","2:3",..: 154 295 73 73
121 94 70 94 295 619 ...
$ bf : Factor w/ 12 levels "1:0.8","1:0.6",..: 2 2 2 2 2 2 2 2 2
2 ...
where gf is the interaction of 'geno' and 'food'; gb is the
interaction of 'geno' and 'block'; and 'bf' is the interaction of
'block' and 'food'. There is no geno:food:block interaction.
I fit the following model. I should also note that lmer2 produces
the same results.
> full.ovn<-lmer(ovn~1+(1|block)+(1|gb)+(1|gf)+(1|bf)+(1|geno)+(1|
food), ov)
> summary(full.ovn)
Linear mixed-effects model fit by REML
Formula: ovn ~ 1 + (1 | block) + (1 | gb) + (1 | gf) + (1 | bf) + (1
| geno) + (1 | food)
Data: ov
AIC BIC logLik MLdeviance REMLdeviance
49097 49146 -24542 49087 49083
Random effects:
Groups Name Variance Std.Dev.
gf (Intercept) 6.43894 2.53751
gb (Intercept) 6.00979 2.45149
geno (Intercept) 3.37417 1.83689
bf (Intercept) 0.84951 0.92169
food (Intercept) 19.30673 4.39394
block (Intercept) 0.35833 0.59861
Residual 22.38789 4.73158
number of obs: 7958, groups: gf, 787; gb, 527; geno, 214; bf, 12;
food, 4; block, 3
Fixed effects:
Estimate Std. Error t value
(Intercept) 21.009 2.249 9.342
And, then run mcmcsamp:
> full.mcmc<-mcmcsamp(full.ovn, n=10000)
> summary(full.mcmc)
Iterations = 1:10000
Thinning interval = 1
Number of chains = 1
Sample size per chain = 10000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
(Intercept) 21.0016 3.82131 0.0382131 0.0350391
log(sigma^2) 3.1106 0.01734 0.0001734 0.0003096
log(gf.(In)) 1.7228 0.08343 0.0008343 0.0019563
log(gb.(In)) 1.8292 0.09193 0.0009193 0.0019939
log(geno.(In)) 1.6577 0.13624 0.0013624 0.0015901
log(bf.(In)) 0.5115 0.62757 0.0062757 0.0119208
log(food.(In)) 3.2965 1.01289 0.0101289 0.0154601
log(blck.(In)) -50.0119 43.12205 0.4312205 5.3344899
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
(Intercept) 13.7807 19.34548 21.0349 22.6964 27.9320
log(sigma^2) 3.0763 3.09890 3.1105 3.1222 3.1446
log(gf.(In)) 1.5619 1.66684 1.7226 1.7782 1.8891
log(gb.(In)) 1.6491 1.76648 1.8300 1.8926 2.0060
log(geno.(In)) 1.3956 1.56445 1.6588 1.7480 1.9267
log(bf.(In)) -0.6863 0.08243 0.5012 0.9256 1.7489
log(food.(In)) 1.6572 2.58508 3.1775 3.8774 5.6251
log(blck.(In)) -152.0792 -61.58789 -38.9542 -19.1747 -0.4419
Warning message:
In glm.fit(x = X, y = Y, weights = weights, start = start, etastart =
etastart, :
algorithm did not converge
and then, the HPD interval:
> ovn.hpd
lower upper
(Intercept) 13.744164 27.856893
log(sigma^2) 3.075533 3.143706
log(gf.(In)) 1.564024 1.890378
log(gb.(In)) 1.650405 2.006905
log(geno.(In)) 1.399007 1.929362
log(bf.(In)) -0.721077 1.706577
log(food.(In)) 1.518141 5.352884
log(blck.(In)) -142.946749 1.565227
attr(,"Probability")
[1] 0.95
>
As you can see, the HPD interval for the 'geno' variance component do
no enclose the estimate for the variance component itself. Further,
the mean value of 'geno' from the MCMC run does not equal the
variance component for 'geno'. For other variance components there
is no problem.
Is the warning message at the end of the summary(full.mcmc) a clue
that there is just poor convergence? Or, is this a bug?
Thanks for the help,
Alan
More information about the R-sig-mixed-models
mailing list