[R] HPDinterval problem

[Seyed Reza Jafarzadeh]

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