thr3ads.net - R help - [R] HPDinterval problem [Apr 2007]

If this information is useful, please help other people find it:
Share via:

Seyed Reza Jafarzadeh

2007-Apr-03 00:26 UTC

[R] HPDinterval problem

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.16Generalized 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

Seyed Reza Jafarzadeh

2007-Apr-04 05:00 UTC

head link

[R] HPDinterval problem

Hi,

I am providing more examples where HPDinterval failed. It seems to be
working OK for (generalized linear mixed) models without crossed
random-effects (m1.17, m1.18, m1.19, m1.20, m1.21, m1.22, and m1.24
below).

Thank you,
Reza

> m1.1 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (1 | prov) +
(1 | pm), data = mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.1, n = 10000))                      lower        upper
(Intercept)    -0.207561922 -0.207561922
pv1o            0.056574609  0.056574609
pv2o            0.023042057  0.023042057
pv1toa          0.026497315  0.026497315
pv2toa         -0.001074887 -0.001074887
sesblf2         0.307805373  0.307805373
sesblf3         0.067866694  0.067866694
sesblf4         0.232652035  0.232652035
log(prov.(In)) -1.948540913 -0.815550437
log(pm.(In))   -4.609549269 -3.008113214
attr(,"Probability")
[1] 0.95

> m1.3 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1toa |
prov) + (1 | pm), data = mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.3, n = 10000))                           lower         upper
(Intercept)        -0.3854582560 -0.3854582560
pv1o                0.0545945359  0.0545945359
pv2o                0.0266911717  0.0266911717
pv1toa              0.0369314516  0.0369314516
pv2toa             -0.0008906397 -0.0008906397
sesblf2             0.3326814534  0.3326814534
sesblf3             0.1012759194  0.1012759194
sesblf4             0.1968001587  0.1968001587
log(prov.(In))     -1.2423994216 -0.0463231047
log(prov.pv1t)     -8.5013756480 -7.3008649434
atanh(prv.(I).pv1) -1.3266358579 -0.4613822430
log(pm.(In))       -4.5813293741 -2.9416249086
attr(,"Probability")
[1] 0.95

> m1.5 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (1 | prov) +
(pv1toa | pm), data = mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.5, n = 10000))                         lower        upper
(Intercept)       -0.298634101 -0.298634101
pv1o               0.056017516  0.056017516
pv2o               0.021658991  0.021658991
pv1toa             0.028086682  0.028086682
pv2toa             0.003447681  0.003447681
sesblf2            0.413727463  0.413727463
sesblf3            0.046766676  0.046766676
sesblf4            0.255977008  0.255977008
log(prov.(In))    -1.875689638 -0.751072995
log(pm.(In))      -3.556592560 -1.722182602
log(pm.pv1t)      -9.709527247 -7.885488338
atanh(pm.(I).pv1) -1.901663364 -0.616765080
attr(,"Probability")
[1] 0.95

> m1.7 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1toa |
prov) + (pv1o | pm), data = mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.7, n = 10000))                          lower        upper
(Intercept)        -0.389923132 -0.389923132
pv1o                0.063098026  0.063098026
pv2o                0.034944761  0.034944761
pv1toa              0.032622126  0.032622126
pv2toa              0.003154919  0.003154919
sesblf2             0.300371141  0.300371141
sesblf3             0.020146759  0.020146759
sesblf4             0.187532056  0.187532056
log(prov.(In))     -1.322210629 -0.131769983
log(prov.pv1t)     -8.411977395 -7.219326685
atanh(prv.(I).pv1) -1.257375358 -0.396660253
log(pm.(In))       -3.671581481 -1.835388518
log(pm.pv1o)       -7.107693068 -5.275740794
atanh(pm.(I).pv1)  -1.679642476 -0.382000422
attr(,"Probability")
[1] 0.95

> m1.8 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o | prov)
+ (pv1toa | pm), data = mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.8, n = 10000))                          lower        upper
(Intercept)        -0.457054707 -0.457054707
pv1o                0.156145534  0.156145534
pv2o                0.024773645  0.024773645
pv1toa              0.024579764  0.024579764
pv2toa              0.001907060  0.001907060
sesblf2             0.394565315  0.394565315
sesblf3             0.061645816  0.061645816
sesblf4             0.259691274  0.259691274
log(prov.(In))     -1.301168272 -0.105499439
log(prov.pv1o)     -5.255142692 -4.057937301
atanh(prv.(I).pv1) -1.851880731 -0.999139040
log(pm.(In))       -3.798743689 -1.968443546
log(pm.pv1t)       -9.788997900 -7.962938034
atanh(pm.(I).pv1)  -1.761670102 -0.480121774
attr(,"Probability")
[1] 0.95

> m1.9 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1toa |
prov) + (pv1toa | pm), data = mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.9, n = 10000))                          lower        upper
(Intercept)        -0.485035838 -0.485035838
pv1o                0.053927806  0.053927806
pv2o                0.027072754  0.027072754
pv1toa              0.036558311  0.036558311
pv2toa              0.004437478  0.004437478
sesblf2             0.470407622  0.470407622
sesblf3             0.094079640  0.094079640
sesblf4             0.240104293  0.240104293
log(prov.(In))     -1.195701809  0.010070569
log(prov.pv1t)     -8.449706678 -7.242412741
atanh(prv.(I).pv1) -1.288497599 -0.436028040
log(pm.(In))       -3.343337299 -1.524187369
log(pm.pv1t)       -9.575747795 -7.757473396
atanh(pm.(I).pv1)  -2.020459923 -0.729784182
attr(,"Probability")
[1] 0.95

> m1.10 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o | prov)
+ (pv1toa | prov) + (1 | pm), data = mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.10, n = 10000))                           lower         upper
(Intercept)        -0.4456071894 -0.4456071894
pv1o                0.1464835742  0.1464835742
pv2o                0.0229765752  0.0229765752
pv1toa              0.0278685330  0.0278685330
pv2toa             -0.0003038598 -0.0003038598
sesblf2             0.3443156778  0.3443156778
sesblf3             0.0944738799  0.0944738799
sesblf4             0.2254927178  0.2254927178
log(prov.(In))     -1.6491223578 -0.4377761632
log(prov.pv1o)     -5.6182267338 -4.4035977132
atanh(prv.(I).pv1) -2.5384179212 -1.6957391764
log(prov.(In))     -2.5659880481 -1.3799410576
log(prov.pv1t)     -9.0668805185 -7.8753819065
atanh(prv.(I).pv1) -2.0180897439 -1.1751174649
log(pm.(In))       -4.4612011753 -2.8135021093
attr(,"Probability")
[1] 0.95

> m1.11 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (1 | prov) +
(pv1o | pm) + (pv1toa | pm), data = mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.11, n = 10000))                        lower       upper
(Intercept)       -0.29985554 -0.29985554
pv1o               0.07375569  0.07375569
pv2o               0.02568464  0.02568464
pv1toa             0.02426277  0.02426277
pv2toa             0.00551527  0.00551527
sesblf2            0.35060701  0.35060701
sesblf3            0.01707800  0.01707800
sesblf4            0.23239228  0.23239228
log(prov.(In))    -2.01160823 -0.87660101
log(pm.(In))      -4.16132557 -2.35116260
log(pm.pv1o)      -6.98272087 -5.17284124
atanh(pm.(I).pv1) -7.98465439 -6.71075618
log(pm.(In))      -3.74016652 -1.93693950
log(pm.pv1t)      -9.66383022 -7.88527146
atanh(pm.(I).pv1) -1.73988647 -0.46569668
attr(,"Probability")
[1] 0.95

> m1.12 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o | prov)
+ (pv1toa | prov) + (pv1o | pm), data = mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.12, n = 10000))                          lower        upper
(Intercept)        -0.443061575 -0.443061575
pv1o                0.145784210  0.145784210
pv2o                0.032899568  0.032899568
pv1toa              0.026326962  0.026326962
pv2toa              0.002251301  0.002251301
sesblf2             0.321052133  0.321052133
sesblf3             0.016162887  0.016162887
sesblf4             0.213160215  0.213160215
log(prov.(In))     -1.558074949 -0.356582122
log(prov.pv1o)     -5.675655643 -4.488932586
atanh(prv.(I).pv1) -2.680508143 -1.837073359
log(prov.(In))     -2.879667257 -1.665153073
log(prov.pv1t)     -8.939603075 -7.728957475
atanh(prv.(I).pv1) -1.979096998 -1.125821675
log(pm.(In))       -3.715644215 -1.907578884
log(pm.pv1o)       -7.149593475 -5.332610229
atanh(pm.(I).pv1)  -1.662717526 -0.365192495
attr(,"Probability")
[1] 0.95

> m1.13 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o | prov)
+ (pv1toa | prov) + (pv1toa | pm), data = mydata[1:1392, ], family =
quasipoisson)
> HPDinterval(mcmcsamp(m1.13, n = 10000))                          lower        upper
(Intercept)        -0.520671363 -0.520671363
pv1o                0.143483258  0.143483258
pv2o                0.023551563  0.023551563
pv1toa              0.028365257  0.028365257
pv2toa              0.004494341  0.004494341
sesblf2             0.400505707  0.400505707
sesblf3             0.072760610  0.072760610
sesblf4             0.252793971  0.252793971
log(prov.(In))     -1.664999620 -0.457571686
log(prov.pv1o)     -5.676457356 -4.482666134
atanh(prv.(I).pv1) -2.422760733 -1.569174953
log(prov.(In))     -2.466638158 -1.254008900
log(prov.pv1t)     -8.913571866 -7.705391788
atanh(prv.(I).pv1) -1.971407906 -1.105821462
log(pm.(In))       -3.708401541 -1.878728047
log(pm.pv1t)       -9.633734032 -7.816379065
atanh(pm.(I).pv1)  -1.760071642 -0.488676259
attr(,"Probability")
[1] 0.95

> m1.14 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o | prov)
+ (pv1o | pm) + (pv1toa | pm), data = mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.14, n = 10000))                          lower        upper
(Intercept)        -0.443656150 -0.443656150
pv1o                0.162193691  0.162193691
pv2o                0.031649149  0.031649149
pv1toa              0.022429405  0.022429405
pv2toa              0.002443997  0.002443997
sesblf2             0.360466084  0.360466084
sesblf3             0.016536949  0.016536949
sesblf4             0.231307640  0.231307640
log(prov.(In))     -1.339206184 -0.149150781
log(prov.pv1o)     -5.335011824 -4.152221345
atanh(prv.(I).pv1) -1.943320599 -1.087748000
log(pm.(In))       -4.267623575 -2.442556980
log(pm.pv1o)       -7.190055027 -5.371607393
atanh(pm.(I).pv1)  -3.831335751 -2.554657428
log(pm.(In))       -4.045019561 -2.242185600
log(pm.pv1t)       -9.952372879 -8.157867980
atanh(pm.(I).pv1)  -1.606404138 -0.325092337
attr(,"Probability")
[1] 0.95

> m1.15 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1toa |
prov) + (pv1o | pm) + (pv1toa | pm), data = mydata[1:1392, ], family =
quasipoisson)
> HPDinterval(mcmcsamp(m1.15, n = 10000))                          lower        upper
(Intercept)        -0.492628035 -0.492628035
pv1o                0.066932526  0.066932526
pv2o                0.032558697  0.032558697
pv1toa              0.031964223  0.031964223
pv2toa              0.006872102  0.006872102
sesblf2             0.462985660  0.462985660
sesblf3             0.059477673  0.059477673
sesblf4             0.228330298  0.228330298
log(prov.(In))     -1.284303557 -0.103539433
log(prov.pv1t)     -8.417909708 -7.204367435
atanh(prv.(I).pv1) -1.282230449 -0.436899596
log(pm.(In))       -3.964190389 -2.154068579
log(pm.pv1o)       -7.288337695 -5.491624377
atanh(pm.(I).pv1)  -2.373590093 -1.106568121
log(pm.(In))       -3.689622102 -1.867271754
log(pm.pv1t)       -9.611438322 -7.798522903
atanh(pm.(I).pv1)  -2.542602654 -1.259583953
attr(,"Probability")
[1] 0.95




> m1.17 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (1 | prov),
data =  mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.17, n = 10000))                       lower        upper
(Intercept)    -0.4556673441  0.013126565
pv1o            0.0453480105  0.064682544
pv2o            0.0151423969  0.034835504
pv1toa          0.0183202799  0.026733287
pv2toa         -0.0020464788  0.006746995
sesblf2         0.2216241893  0.441465720
sesblf3        -0.0005681114  0.209583281
sesblf4         0.1468518151  0.356648966
log(prov.(In)) -1.9257817242 -0.683522096
attr(,"Probability")
[1] 0.95

> m1.18 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (1 | pm),
data =  mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.18, n = 10000))                   lower        upper
(Intercept)  -0.30685638  0.228099238
pv1o          0.07232589  0.090089888
pv2o          0.03628791  0.055208902
pv1toa        0.01722480  0.026059233
pv2toa       -0.01216203 -0.002886872
sesblf2      -0.06651472  0.679455774
sesblf3      -0.32185743  0.425081143
sesblf4      -0.19286322  0.550831225
log(pm.(In)) -4.34513788 -2.000295046
attr(,"Probability")
[1] 0.95

> m1.19 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o |
prov), data =  mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.19, n = 10000))                          lower        upper
(Intercept)        -0.732466301 -0.095762021
pv1o                0.109101825  0.205227082
pv2o                0.014971677  0.035573998
pv1toa              0.015082087  0.024262828
pv2toa             -0.003380130  0.005867948
sesblf2             0.252479830  0.479557058
sesblf3             0.023607887  0.241040795
sesblf4             0.140224883  0.358214042
log(prov.(In))     -1.188516571  0.105022406
log(prov.pv1o)     -5.250761065 -3.673800387
atanh(prv.(I).pv1) -1.895520991 -0.894677411
attr(,"Probability")
[1] 0.95

> m1.20 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1toa |
prov), data =  mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.20, n = 10000))                          lower        upper
(Intercept)        -0.654205909 -0.062846199
pv1o                0.044454141  0.063085006
pv2o                0.018996380  0.038600912
pv1toa              0.023692650  0.041228419
pv2toa             -0.001622876  0.006471938
sesblf2             0.253318503  0.460385395
sesblf3             0.041646544  0.232981172
sesblf4             0.114749796  0.309772745
log(prov.(In))     -1.179070095  0.115172146
log(prov.pv1t)     -8.502052140 -7.093668453
atanh(prv.(I).pv1) -1.406595407 -0.442011756
attr(,"Probability")
[1] 0.95

> m1.21 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o | prov)
+ (pv1toa | prov), data =  mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.21, n = 10000))                          lower        upper
(Intercept)        -0.809795120 -0.092818083
pv1o                0.102320721  0.195798084
pv2o                0.013958996  0.035777516
pv1toa              0.015285309  0.032155176
pv2toa             -0.001270502  0.008386427
sesblf2             0.265843858  0.488668321
sesblf3             0.022096224  0.246423841
sesblf4             0.146323603  0.368723416
log(prov.(In))     -2.053426653 -0.003235248
log(prov.pv1o)     -5.494391785 -3.713515236
atanh(prv.(I).pv1) -6.542210732 -0.849163554
log(prov.(In))     -2.693832187 -0.335722676
log(prov.pv1t)     -9.432610814 -7.470030780
atanh(prv.(I).pv1) -4.334208695 -0.366701490
attr(,"Probability")
[1] 0.95

> m1.22 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o | pm),
data =  mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.22, n = 10000))                         lower        upper
(Intercept)       -0.464495795  0.251100466
pv1o               0.058394653  0.150951592
pv2o               0.044325579  0.064938826
pv1toa             0.012090092  0.022003334
pv2toa            -0.007826689  0.002510138
sesblf2           -0.198332573  0.732862207
sesblf3           -0.451325552  0.428932303
sesblf4           -0.282425435  0.600581823
log(pm.(In))      -3.340172245 -0.892699307
log(pm.pv1o)      -6.211385523 -4.186685594
atanh(pm.(I).pv1) -1.804149382 -0.051858193
attr(,"Probability")
[1] 0.95

> m1.24 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o | pm) +
(pv1toa | pm), data =  mydata[1:1392, ], family = quasipoisson)
> HPDinterval(mcmcsamp(m1.24, n = 10000))                          lower        upper
(Intercept)        -0.618999654  0.371571883
pv1o                0.058341562  0.153168044
pv2o                0.041622796  0.062613888
pv1toa              0.006833982  0.025585021
pv2toa             -0.006074951  0.004975857
sesblf2            -0.703034251  1.379272724
sesblf3            -0.494385351  0.459530213
sesblf4            -0.296029582  0.674519473
log(pm.(In))       -3.918382813 -0.627437580
log(pm.pv1o)       -6.248238619 -4.142139977
atanh(pm.(I).pv1)  -7.304481224 -0.036363295
log(pm.(In))       -6.083394919 -0.025413685
log(pm.pv1t)      -10.436085207 -7.435829249
atanh(pm.(I).pv1)  -6.066684183  3.950585563
attr(,"Probability")
[1] 0.95






On 4/1/07, Seyed Reza Jafarzadeh <srjafarzadeh at gmail.com>
wrote:> 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
>
>
>
>
> Thanks,
> Reza
>

Seyed Reza Jafarzadeh

2007-Apr-04 06:59 UTC

head link

[R] HPDinterval problem

Hi,

Sorry for multiple postings. Please consider the following reproducible example:

> grp1 <- rep(1:24, times = 24, each = 1)
> grp2 <- rep(1:12, times = 2, each = 24)
> set.seed(1)
> out <- as.integer(abs(rnorm(576))*10)
> x1 <- out [25:552]
> x2 <- out [49:576]
> mydf <- data.frame(cbind(out[1:528], x1, x2, grp1[49:576],
grp2[49:576]))
> names(mydf) <- c("out", "x1", "x2",
"grp1", "grp2")
> me1 <- lmer(out ~ x1 + x2 + (x1 | grp1) + (x2 | grp2), data = mydf,
family = quasipoisson)
> me1Generalized linear mixed model fit using Laplace
Formula: out ~ x1 + x2 + (x1 | grp1) + (x2 | grp2)
   Data: mydf
 Family: quasipoisson(log link)
  AIC  BIC logLik deviance
 2559 2598  -1271     2541
Random effects:
 Groups   Name        Variance  Std.Dev. Corr
 grp1     (Intercept) 0.3250162 0.570102
          x1          0.0021856 0.046750 -0.813
 grp2     (Intercept) 0.3183619 0.564236
          x2          0.0038276 0.061868 -0.940
 Residual             4.3279637 2.080376
number of obs: 528, groups: grp1, 24; grp2, 12

Fixed effects:
             Estimate Std. Error t value
(Intercept)  2.036473   0.212687   9.575
x1          -0.001159   0.011234  -0.103
x2          -0.006830   0.018860  -0.362

Correlation of Fixed Effects:
   (Intr) x1
x1 -0.487
x2 -0.751  0.007
> HPDinterval(mcmcsamp(me1, 20000))                         lower        upper
(Intercept)        2.036472792  2.036472792
x1                -0.001158922 -0.001158922
x2                -0.006830165 -0.006830165
log(grp1.(In))    -3.254548062 -2.052213965
log(grp1.x1)      -8.351051505 -7.149742494
atanh(gr1.(I).x1) -1.630982865 -0.783008788
log(grp2.(In))    -3.312492769 -1.509030356
log(grp2.x2)      -7.739344906 -5.926156515
atanh(gr2.(I).x2) -2.564511232 -1.291401604
attr(,"Probability")
[1] 0.95


Thanks,
Reza



On 4/3/07, Seyed Reza Jafarzadeh <srjafarzadeh at gmail.com>
wrote:> Hi,
>
> I am providing more examples where HPDinterval failed. It seems to be
> working OK for (generalized linear mixed) models without crossed
> random-effects (m1.17, m1.18, m1.19, m1.20, m1.21, m1.22, and m1.24
> below).
>
> Thank you,
> Reza
>
>
> > m1.1 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (1 |
prov) + (1 | pm), data = mydata[1:1392, ], family = quasipoisson)
> > HPDinterval(mcmcsamp(m1.1, n = 10000))
>                       lower        upper
> (Intercept)    -0.207561922 -0.207561922
> pv1o            0.056574609  0.056574609
> pv2o            0.023042057  0.023042057
> pv1toa          0.026497315  0.026497315
> pv2toa         -0.001074887 -0.001074887
> sesblf2         0.307805373  0.307805373
> sesblf3         0.067866694  0.067866694
> sesblf4         0.232652035  0.232652035
> log(prov.(In)) -1.948540913 -0.815550437
> log(pm.(In))   -4.609549269 -3.008113214
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.3 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1toa |
prov) + (1 | pm), data = mydata[1:1392, ], family = quasipoisson)
> > HPDinterval(mcmcsamp(m1.3, n = 10000))
>                            lower         upper
> (Intercept)        -0.3854582560 -0.3854582560
> pv1o                0.0545945359  0.0545945359
> pv2o                0.0266911717  0.0266911717
> pv1toa              0.0369314516  0.0369314516
> pv2toa             -0.0008906397 -0.0008906397
> sesblf2             0.3326814534  0.3326814534
> sesblf3             0.1012759194  0.1012759194
> sesblf4             0.1968001587  0.1968001587
> log(prov.(In))     -1.2423994216 -0.0463231047
> log(prov.pv1t)     -8.5013756480 -7.3008649434
> atanh(prv.(I).pv1) -1.3266358579 -0.4613822430
> log(pm.(In))       -4.5813293741 -2.9416249086
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.5 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (1 |
prov) + (pv1toa | pm), data = mydata[1:1392, ], family = quasipoisson)
> > HPDinterval(mcmcsamp(m1.5, n = 10000))
>                          lower        upper
> (Intercept)       -0.298634101 -0.298634101
> pv1o               0.056017516  0.056017516
> pv2o               0.021658991  0.021658991
> pv1toa             0.028086682  0.028086682
> pv2toa             0.003447681  0.003447681
> sesblf2            0.413727463  0.413727463
> sesblf3            0.046766676  0.046766676
> sesblf4            0.255977008  0.255977008
> log(prov.(In))    -1.875689638 -0.751072995
> log(pm.(In))      -3.556592560 -1.722182602
> log(pm.pv1t)      -9.709527247 -7.885488338
> atanh(pm.(I).pv1) -1.901663364 -0.616765080
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.7 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1toa |
prov) + (pv1o | pm), data = mydata[1:1392, ], family = quasipoisson)
> > HPDinterval(mcmcsamp(m1.7, n = 10000))
>                           lower        upper
> (Intercept)        -0.389923132 -0.389923132
> pv1o                0.063098026  0.063098026
> pv2o                0.034944761  0.034944761
> pv1toa              0.032622126  0.032622126
> pv2toa              0.003154919  0.003154919
> sesblf2             0.300371141  0.300371141
> sesblf3             0.020146759  0.020146759
> sesblf4             0.187532056  0.187532056
> log(prov.(In))     -1.322210629 -0.131769983
> log(prov.pv1t)     -8.411977395 -7.219326685
> atanh(prv.(I).pv1) -1.257375358 -0.396660253
> log(pm.(In))       -3.671581481 -1.835388518
> log(pm.pv1o)       -7.107693068 -5.275740794
> atanh(pm.(I).pv1)  -1.679642476 -0.382000422
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.8 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o |
prov) + (pv1toa | pm), data = mydata[1:1392, ], family = quasipoisson)
> > HPDinterval(mcmcsamp(m1.8, n = 10000))
>                           lower        upper
> (Intercept)        -0.457054707 -0.457054707
> pv1o                0.156145534  0.156145534
> pv2o                0.024773645  0.024773645
> pv1toa              0.024579764  0.024579764
> pv2toa              0.001907060  0.001907060
> sesblf2             0.394565315  0.394565315
> sesblf3             0.061645816  0.061645816
> sesblf4             0.259691274  0.259691274
> log(prov.(In))     -1.301168272 -0.105499439
> log(prov.pv1o)     -5.255142692 -4.057937301
> atanh(prv.(I).pv1) -1.851880731 -0.999139040
> log(pm.(In))       -3.798743689 -1.968443546
> log(pm.pv1t)       -9.788997900 -7.962938034
> atanh(pm.(I).pv1)  -1.761670102 -0.480121774
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.9 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1toa |
prov) + (pv1toa | pm), data = mydata[1:1392, ], family = quasipoisson)
> > HPDinterval(mcmcsamp(m1.9, n = 10000))
>                           lower        upper
> (Intercept)        -0.485035838 -0.485035838
> pv1o                0.053927806  0.053927806
> pv2o                0.027072754  0.027072754
> pv1toa              0.036558311  0.036558311
> pv2toa              0.004437478  0.004437478
> sesblf2             0.470407622  0.470407622
> sesblf3             0.094079640  0.094079640
> sesblf4             0.240104293  0.240104293
> log(prov.(In))     -1.195701809  0.010070569
> log(prov.pv1t)     -8.449706678 -7.242412741
> atanh(prv.(I).pv1) -1.288497599 -0.436028040
> log(pm.(In))       -3.343337299 -1.524187369
> log(pm.pv1t)       -9.575747795 -7.757473396
> atanh(pm.(I).pv1)  -2.020459923 -0.729784182
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.10 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o |
prov) + (pv1toa | prov) + (1 | pm), data = mydata[1:1392, ], family =
quasipoisson)
> > HPDinterval(mcmcsamp(m1.10, n = 10000))
>                            lower         upper
> (Intercept)        -0.4456071894 -0.4456071894
> pv1o                0.1464835742  0.1464835742
> pv2o                0.0229765752  0.0229765752
> pv1toa              0.0278685330  0.0278685330
> pv2toa             -0.0003038598 -0.0003038598
> sesblf2             0.3443156778  0.3443156778
> sesblf3             0.0944738799  0.0944738799
> sesblf4             0.2254927178  0.2254927178
> log(prov.(In))     -1.6491223578 -0.4377761632
> log(prov.pv1o)     -5.6182267338 -4.4035977132
> atanh(prv.(I).pv1) -2.5384179212 -1.6957391764
> log(prov.(In))     -2.5659880481 -1.3799410576
> log(prov.pv1t)     -9.0668805185 -7.8753819065
> atanh(prv.(I).pv1) -2.0180897439 -1.1751174649
> log(pm.(In))       -4.4612011753 -2.8135021093
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.11 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (1 |
prov) + (pv1o | pm) + (pv1toa | pm), data = mydata[1:1392, ], family =
quasipoisson)
> > HPDinterval(mcmcsamp(m1.11, n = 10000))
>                         lower       upper
> (Intercept)       -0.29985554 -0.29985554
> pv1o               0.07375569  0.07375569
> pv2o               0.02568464  0.02568464
> pv1toa             0.02426277  0.02426277
> pv2toa             0.00551527  0.00551527
> sesblf2            0.35060701  0.35060701
> sesblf3            0.01707800  0.01707800
> sesblf4            0.23239228  0.23239228
> log(prov.(In))    -2.01160823 -0.87660101
> log(pm.(In))      -4.16132557 -2.35116260
> log(pm.pv1o)      -6.98272087 -5.17284124
> atanh(pm.(I).pv1) -7.98465439 -6.71075618
> log(pm.(In))      -3.74016652 -1.93693950
> log(pm.pv1t)      -9.66383022 -7.88527146
> atanh(pm.(I).pv1) -1.73988647 -0.46569668
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.12 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o |
prov) + (pv1toa | prov) + (pv1o | pm), data = mydata[1:1392, ], family =
quasipoisson)
> > HPDinterval(mcmcsamp(m1.12, n = 10000))
>                           lower        upper
> (Intercept)        -0.443061575 -0.443061575
> pv1o                0.145784210  0.145784210
> pv2o                0.032899568  0.032899568
> pv1toa              0.026326962  0.026326962
> pv2toa              0.002251301  0.002251301
> sesblf2             0.321052133  0.321052133
> sesblf3             0.016162887  0.016162887
> sesblf4             0.213160215  0.213160215
> log(prov.(In))     -1.558074949 -0.356582122
> log(prov.pv1o)     -5.675655643 -4.488932586
> atanh(prv.(I).pv1) -2.680508143 -1.837073359
> log(prov.(In))     -2.879667257 -1.665153073
> log(prov.pv1t)     -8.939603075 -7.728957475
> atanh(prv.(I).pv1) -1.979096998 -1.125821675
> log(pm.(In))       -3.715644215 -1.907578884
> log(pm.pv1o)       -7.149593475 -5.332610229
> atanh(pm.(I).pv1)  -1.662717526 -0.365192495
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.13 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o |
prov) + (pv1toa | prov) + (pv1toa | pm), data = mydata[1:1392, ], family =
quasipoisson)
> > HPDinterval(mcmcsamp(m1.13, n = 10000))
>                           lower        upper
> (Intercept)        -0.520671363 -0.520671363
> pv1o                0.143483258  0.143483258
> pv2o                0.023551563  0.023551563
> pv1toa              0.028365257  0.028365257
> pv2toa              0.004494341  0.004494341
> sesblf2             0.400505707  0.400505707
> sesblf3             0.072760610  0.072760610
> sesblf4             0.252793971  0.252793971
> log(prov.(In))     -1.664999620 -0.457571686
> log(prov.pv1o)     -5.676457356 -4.482666134
> atanh(prv.(I).pv1) -2.422760733 -1.569174953
> log(prov.(In))     -2.466638158 -1.254008900
> log(prov.pv1t)     -8.913571866 -7.705391788
> atanh(prv.(I).pv1) -1.971407906 -1.105821462
> log(pm.(In))       -3.708401541 -1.878728047
> log(pm.pv1t)       -9.633734032 -7.816379065
> atanh(pm.(I).pv1)  -1.760071642 -0.488676259
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.14 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o |
prov) + (pv1o | pm) + (pv1toa | pm), data = mydata[1:1392, ], family =
quasipoisson)
> > HPDinterval(mcmcsamp(m1.14, n = 10000))
>                           lower        upper
> (Intercept)        -0.443656150 -0.443656150
> pv1o                0.162193691  0.162193691
> pv2o                0.031649149  0.031649149
> pv1toa              0.022429405  0.022429405
> pv2toa              0.002443997  0.002443997
> sesblf2             0.360466084  0.360466084
> sesblf3             0.016536949  0.016536949
> sesblf4             0.231307640  0.231307640
> log(prov.(In))     -1.339206184 -0.149150781
> log(prov.pv1o)     -5.335011824 -4.152221345
> atanh(prv.(I).pv1) -1.943320599 -1.087748000
> log(pm.(In))       -4.267623575 -2.442556980
> log(pm.pv1o)       -7.190055027 -5.371607393
> atanh(pm.(I).pv1)  -3.831335751 -2.554657428
> log(pm.(In))       -4.045019561 -2.242185600
> log(pm.pv1t)       -9.952372879 -8.157867980
> atanh(pm.(I).pv1)  -1.606404138 -0.325092337
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.15 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1toa
| prov) + (pv1o | pm) + (pv1toa | pm), data = mydata[1:1392, ], family =
quasipoisson)
> > HPDinterval(mcmcsamp(m1.15, n = 10000))
>                           lower        upper
> (Intercept)        -0.492628035 -0.492628035
> pv1o                0.066932526  0.066932526
> pv2o                0.032558697  0.032558697
> pv1toa              0.031964223  0.031964223
> pv2toa              0.006872102  0.006872102
> sesblf2             0.462985660  0.462985660
> sesblf3             0.059477673  0.059477673
> sesblf4             0.228330298  0.228330298
> log(prov.(In))     -1.284303557 -0.103539433
> log(prov.pv1t)     -8.417909708 -7.204367435
> atanh(prv.(I).pv1) -1.282230449 -0.436899596
> log(pm.(In))       -3.964190389 -2.154068579
> log(pm.pv1o)       -7.288337695 -5.491624377
> atanh(pm.(I).pv1)  -2.373590093 -1.106568121
> log(pm.(In))       -3.689622102 -1.867271754
> log(pm.pv1t)       -9.611438322 -7.798522903
> atanh(pm.(I).pv1)  -2.542602654 -1.259583953
> attr(,"Probability")
> [1] 0.95
>
>
>
>
>
> > m1.17 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (1 |
prov), data =  mydata[1:1392, ], family = quasipoisson)
> > HPDinterval(mcmcsamp(m1.17, n = 10000))
>                        lower        upper
> (Intercept)    -0.4556673441  0.013126565
> pv1o            0.0453480105  0.064682544
> pv2o            0.0151423969  0.034835504
> pv1toa          0.0183202799  0.026733287
> pv2toa         -0.0020464788  0.006746995
> sesblf2         0.2216241893  0.441465720
> sesblf3        -0.0005681114  0.209583281
> sesblf4         0.1468518151  0.356648966
> log(prov.(In)) -1.9257817242 -0.683522096
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.18 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (1 |
pm), data =  mydata[1:1392, ], family = quasipoisson)
> > HPDinterval(mcmcsamp(m1.18, n = 10000))
>                    lower        upper
> (Intercept)  -0.30685638  0.228099238
> pv1o          0.07232589  0.090089888
> pv2o          0.03628791  0.055208902
> pv1toa        0.01722480  0.026059233
> pv2toa       -0.01216203 -0.002886872
> sesblf2      -0.06651472  0.679455774
> sesblf3      -0.32185743  0.425081143
> sesblf4      -0.19286322  0.550831225
> log(pm.(In)) -4.34513788 -2.000295046
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.19 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o |
prov), data =  mydata[1:1392, ], family = quasipoisson)
> > HPDinterval(mcmcsamp(m1.19, n = 10000))
>                           lower        upper
> (Intercept)        -0.732466301 -0.095762021
> pv1o                0.109101825  0.205227082
> pv2o                0.014971677  0.035573998
> pv1toa              0.015082087  0.024262828
> pv2toa             -0.003380130  0.005867948
> sesblf2             0.252479830  0.479557058
> sesblf3             0.023607887  0.241040795
> sesblf4             0.140224883  0.358214042
> log(prov.(In))     -1.188516571  0.105022406
> log(prov.pv1o)     -5.250761065 -3.673800387
> atanh(prv.(I).pv1) -1.895520991 -0.894677411
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.20 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1toa
| prov), data =  mydata[1:1392, ], family = quasipoisson)
> > HPDinterval(mcmcsamp(m1.20, n = 10000))
>                           lower        upper
> (Intercept)        -0.654205909 -0.062846199
> pv1o                0.044454141  0.063085006
> pv2o                0.018996380  0.038600912
> pv1toa              0.023692650  0.041228419
> pv2toa             -0.001622876  0.006471938
> sesblf2             0.253318503  0.460385395
> sesblf3             0.041646544  0.232981172
> sesblf4             0.114749796  0.309772745
> log(prov.(In))     -1.179070095  0.115172146
> log(prov.pv1t)     -8.502052140 -7.093668453
> atanh(prv.(I).pv1) -1.406595407 -0.442011756
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.21 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o |
prov) + (pv1toa | prov), data =  mydata[1:1392, ], family = quasipoisson)
> > HPDinterval(mcmcsamp(m1.21, n = 10000))
>                           lower        upper
> (Intercept)        -0.809795120 -0.092818083
> pv1o                0.102320721  0.195798084
> pv2o                0.013958996  0.035777516
> pv1toa              0.015285309  0.032155176
> pv2toa             -0.001270502  0.008386427
> sesblf2             0.265843858  0.488668321
> sesblf3             0.022096224  0.246423841
> sesblf4             0.146323603  0.368723416
> log(prov.(In))     -2.053426653 -0.003235248
> log(prov.pv1o)     -5.494391785 -3.713515236
> atanh(prv.(I).pv1) -6.542210732 -0.849163554
> log(prov.(In))     -2.693832187 -0.335722676
> log(prov.pv1t)     -9.432610814 -7.470030780
> atanh(prv.(I).pv1) -4.334208695 -0.366701490
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.22 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o |
pm), data =  mydata[1:1392, ], family = quasipoisson)
> > HPDinterval(mcmcsamp(m1.22, n = 10000))
>                          lower        upper
> (Intercept)       -0.464495795  0.251100466
> pv1o               0.058394653  0.150951592
> pv2o               0.044325579  0.064938826
> pv1toa             0.012090092  0.022003334
> pv2toa            -0.007826689  0.002510138
> sesblf2           -0.198332573  0.732862207
> sesblf3           -0.451325552  0.428932303
> sesblf4           -0.282425435  0.600581823
> log(pm.(In))      -3.340172245 -0.892699307
> log(pm.pv1o)      -6.211385523 -4.186685594
> atanh(pm.(I).pv1) -1.804149382 -0.051858193
> attr(,"Probability")
> [1] 0.95
>
>
> > m1.24 <- lmer(o ~ pv1o + pv2o + pv1toa + pv2toa + sesblf + (pv1o |
pm) + (pv1toa | pm), data =  mydata[1:1392, ], family = quasipoisson)
> > HPDinterval(mcmcsamp(m1.24, n = 10000))
>                           lower        upper
> (Intercept)        -0.618999654  0.371571883
> pv1o                0.058341562  0.153168044
> pv2o                0.041622796  0.062613888
> pv1toa              0.006833982  0.025585021
> pv2toa             -0.006074951  0.004975857
> sesblf2            -0.703034251  1.379272724
> sesblf3            -0.494385351  0.459530213
> sesblf4            -0.296029582  0.674519473
> log(pm.(In))       -3.918382813 -0.627437580
> log(pm.pv1o)       -6.248238619 -4.142139977
> atanh(pm.(I).pv1)  -7.304481224 -0.036363295
> log(pm.(In))       -6.083394919 -0.025413685
> log(pm.pv1t)      -10.436085207 -7.435829249
> atanh(pm.(I).pv1)  -6.066684183  3.950585563
> attr(,"Probability")
> [1] 0.95
>
>
>
>
>
>
> On 4/1/07, Seyed Reza Jafarzadeh <srjafarzadeh at gmail.com> wrote:
> > 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
> >
> >
> >
> >
> > Thanks,
> > Reza
> >
>

Possibly Parallel Threads

Search for more apparently analagous threads

R help - Apr 2007 - HPDinterval problem

[R] HPDinterval problem

[R] HPDinterval problem

[R] HPDinterval problem

Possibly Parallel Threads