[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