CODEWITHSUNDEEP
rlme4mixed-models
Leprechault
Leprechault

Reputation: 1823

Confidence Intervals for negative binomial GLMM: uncertainty in the variance parameters

I'd like to calculate the 95 confidence Intervals for prediction in negative binomial GLMMs, using lme4 package. But, a cause of no option for computing standard errors of predictions by default I try to:

#Package
library(lme4)

# Create a data set
set.seed(101)
dd <- expand.grid(f1 = factor(1:3),
                  f2 = LETTERS[1:2], g=1:9, rep=1:15,
                  f3 = rnorm(10,25),
          KEEP.OUT.ATTRS=FALSE)
mu <- 5*(-4 + with(dd, as.integer(f1) + 4*as.numeric(f2)))          
dd$y <- rnbinom(nrow(dd), mu = mu, size = 0.5)
str(dd)

# Create a negative binomial generalized mixed model
m.nb <- glmer.nb(y ~ f1 + f2 + poly(f3,2) + (1|g), data=dd, verbose=TRUE)
m.nb

#Compute the confidence interval for model predictions
mm<-model.matrix(~f1+f2+poly(f3,2),dd)
dd$y_est<-mm%*%fixef(m.nb) #predicted values
pvar <- diag(mm %*% tcrossprod(vcov(m.nb),mm))
CIupper <- dd$y_est+sqrt(pvar) * 1.96 
CIlower <- dd$y_est-sqrt(pvar) * 1.96

But I am not so sure if this is an efficient method for incorporating uncertainty in the variance parameters, especially using a poly() variable. Please, any help with this?

Upvotes: 0

Views: 515

Answers (1)

Leprechault
Leprechault

Reputation: 1823

Compute the confidence interval for model predictions in the code above is OK, but it needs to apply the inverse link function (exp()) to the dd$y_estand pvar + sigma(m.nb)^2:

mm<-model.matrix(~f1+f2+poly(f3,2),dd)
dd$y_est<-exp(mm%*%fixef(m.nb)) #predicted values
pvar <- diag(mm %*% tcrossprod(vcov(m.nb),mm))
CIupper <- dd$y_est+sqrt(pvar+sigma(m.nb)^2) * 1.96 
CIlower <- dd$y_est-sqrt(pvar+sigma(m.nb)^2) * 1.96

Upvotes: 1

Related Questions