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:
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"
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0 0.0 1.0 2.3 3.0 30.0
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 0.000 1.000 2.322 3.000 30.000
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 0.000 1.000 2.315 3.000 30.000
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 4.00 7.00 11.99 15.00 108.00
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)
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)
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