[R] HPDinterval problem

From: Seyed Reza Jafarzadeh <srjafarzadeh_at_gmail.com>
Date: Tue 03 Apr 2007 - 00:26:50 GMT

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. Received on Tue Apr 03 10:33:39 2007

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Tue 03 Apr 2007 - 01:30:37 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.