[R-sig-ME] confidence intervals with mvrnorm - upper value equal to inf
Alessandra Bielli
b|e|||@@|e@@@ndr@ @end|ng |rom gm@||@com
Sat Oct 10 13:28:30 CEST 2020
Dear list,
I am trying to predict a value and CI for two different treatments from a
glmmTMB fitted model using posterior predictive simulations (mvrnorm
function in the MASS package) as in the Salamander example, Brooks 2017
appendix B
<https://www.biorxiv.org/content/biorxiv/suppl/2017/05/01/132753.DC1/132753-2.pdf>.
The dependent variable is a count, majority of values are zeros but some
positive values appear.
m1 <- glmmTMB(Dolphins.TOT ~ Used_piNgers + offset(log(Effort)) +
(1|Trip_Code), ziformula =~1,
data=x, family = "truncated_poisson")
newdata0 = with(x,
expand.grid(
Used_piNgers = c("Y","N"),
Effort=1))
X.cond = model.matrix(lme4::nobars(formula(m1)[-2]), newdata0)
beta.cond = fixef(m1)$cond
pred.cond = X.cond %*% beta.cond
ziformula = m1$modelInfo$allForm$ziformula
X.zi = model.matrix(lme4::nobars(ziformula), newdata0)
beta.zi = fixef(m1)$zi
pred.zi = X.zi %*% beta.zi
pred.ucount = exp(pred.cond)*(1-plogis(pred.zi))
set.seed(101)
pred.condpar.psim = mvrnorm(1000,mu=beta.cond,Sigma=vcov(m1)$cond)
pred.cond.psim = X.cond %*% t(pred.condpar.psim)
pred.zipar.psim = mvrnorm(1000,mu=beta.zi,Sigma=vcov(m1)$zi)
pred.zi.psim = X.zi %*% t(pred.zipar.psim)
pred.ucount.psim = exp(pred.cond.psim)*(1-plogis(pred.zi.psim))
ci.ucount = t(apply(pred.ucount.psim,1,quantile,c(0.025,0.975)))
ci.ucount = data.frame(ci.ucount)
names(ci.ucount) = c("ucount.low","ucount.high")
pred.ucount = data.frame(newdata0, pred.ucount, ci.ucount)
For my upper CI, I get a value equal to Inf:
Used_piNgers Effort pred.ucount ucount.low ucount.high
1 Y 1 6.758889e-11 0.00000000 Inf
2 N 1 1.575418e-02 0.00223033 0.1096139
Is the Inf caused by the very low variability of values in my dataset? I
tried to lower the upper bound of the CI ci.ucount =
t(apply(pred.ucount.psim,1,quantile,c(0.025,0.975))) and only when reaching
0.475 ci.ucount = t(apply(pred.ucount.psim,1,quantile,c(0.025,0.475))) I
obtained:
Used_piNgers Effort pred.ucount ucount.low ucount.high
1 Y 1 6.758889e-11 0.00000000 7.117465e+12
2 N 1 1.575418e-02 0.00223033 1.454579e-02
I found a related post
<https://stackoverflow.com/questions/38272798/bootstrap-confidence-interval-with-inf-in-final-estimates-boot-dplyr-package>
but
the explanation is not clear to me. I would like to publish these results
and I would like to know:
1. is this a sign that something is wrong? if yes, what is it?
2. if nothing is wrong, what does the Inf mean and what's the best way
to report it and plot it in a publication?
I also posted this question on Cross validated
https://stats.stackexchange.com/questions/491196/bootstrap-confidence-interval-with-mvrnorm-upper-value-equal-to-inf
Thanks,
Alessandra
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list