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
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 >
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 > > >