Hi, Can anyone tell me why I am not getting the correct intervals for fixed effect terms for the following generalized linear mixed model from HPDinterval:
> sessionInfo() R version 2.4.1 (2006-12-18) i386-pc-mingw32 locale: LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252 attached base packages: [1] "stats" "graphics" "grDevices" "utils" "datasets" "methods" "base" other attached packages: coda lme4 Matrix lattice "0.10-7" "0.9975-13" "0.9975-11" "0.14-17" > summary(o[1:1392]) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.0 0.0 1.0 2.3 3.0 30.0 > summary(pv1o[1:1392]) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.000 0.000 1.000 2.322 3.000 30.000 > summary(pv2o[1:1392]) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.000 0.000 1.000 2.315 3.000 30.000 > summary(pv1toa[1:1392]) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00 4.00 7.00 11.99 15.00 108.00 > summary(pv2toa[1:1392]) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00 4.00 7.00 11.94 15.00 108.00 > m1.16 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o | prov) + > (pv1toa | prov) + (pv1o | pm) + (pv1toa | pm), data = mydata[1:1392,], family > = quasipoisson) > m1.16 Generalized linear mixed model fit using Laplace Formula: o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o | prov) + (pv1toa | prov) + (pv1o | pm) + (pv1toa | pm) Data: mydata[1:1392, ] Family: quasipoisson(log link) AIC BIC logLik deviance 2285 2390 -1123 2245 Random effects: Groups Name Variance Std.Dev. Corr prov (Intercept) 0.68110734 0.825292 pv1o 0.01182079 0.108723 -0.927 prov (Intercept) 0.09896363 0.314585 pv1toa 0.00029002 0.017030 -0.182 pm (Intercept) 0.05023998 0.224143 pv1o 0.00234292 0.048404 -0.789 pm (Intercept) 0.01918719 0.138518 pv1toa 0.00011984 0.010947 -0.086 Residual 1.54376281 1.242483 number of obs: 1392, groups: prov, 24; prov, 24; pm, 12; pm, 12 Fixed effects: Estimate Std. Error t value (Intercept) -0.372258 0.233326 -1.595 pv1o 0.151591 0.028778 5.268 pv2o 0.029524 0.007247 4.074 pv1toa 0.025669 0.006221 4.126 pv2toa 0.004344 0.003755 1.157 sesblf2 0.074507 0.186258 0.400 sesblf3 -0.037145 0.188021 -0.198 sesblf4 0.155999 0.187175 0.833 Correlation of Fixed Effects: (Intr) pv1o pv2o pv1toa pv2toa ssblf2 ssblf3 pv1o -0.638 pv2o -0.036 -0.099 pv1toa -0.073 -0.050 -0.029 pv2toa -0.043 -0.035 -0.156 -0.458 sesblf2 -0.411 -0.009 0.040 0.002 0.012 sesblf3 -0.412 -0.005 0.039 -0.002 0.037 0.516 sesblf4 -0.413 -0.006 0.044 0.003 0.028 0.513 0.514 > s1.16 <- mcmcsamp(m1.16, n = 100000) > HPDinterval(s1.16) lower upper (Intercept) -0.372258256 -0.372258256 pv1o 0.151590854 0.151590854 pv2o 0.029524315 0.029524315 pv1toa 0.025668727 0.025668727 pv2toa 0.004343653 0.004343653 sesblf2 0.074507427 0.074507427 sesblf3 -0.037144631 -0.037144631 sesblf4 0.155998825 0.155998825 log(prov.(In)) -1.547675380 -0.345723770 log(prov.pv1o) -5.610048117 -4.407086692 atanh(prv.(I).pv1) -2.509960360 -1.663905782 log(prov.(In)) -4.030294678 -2.823797787 log(prov.pv1t) -9.370781684 -8.165302813 atanh(prv.(I).pv1) -1.146944941 -0.289800204 log(pm.(In)) -4.420270387 -2.597929912 log(pm.pv1o) -7.227500164 -5.401277510 atanh(pm.(I).pv1) -2.172644329 -0.886969199 log(pm.(In)) -6.071675906 -4.252728431 log(pm.pv1t) -10.230334351 -8.403082501 atanh(pm.(I).pv1) -0.810182999 0.503799332 attr(,"Probability") [1] 0.95 If I use only one grouping factor (either "prov" or "pm") in the model, it looks to be ok, but it doesn't seem right with two. Thanks, Reza ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.