CODEWITHSUNDEEP
rglmconfidence-intervalstatistics-bootstrap
Jamie Murphy
Jamie Murphy

Reputation: 1

Creating confidence intervals for regression curve in GLMM using Bootstrapping

Have created a GLMM model and plotted the predicted probabilities of each factor. However, I cannot fathom how to create confidence intervals using the BootMer function. I keep getting the error message cannot simulate from non integer prior weights. I'm hoping someone would be able to help? Thanks in advance.

glmm1 <- glmer(cbind(Total_Dead, Total_Collected - Total_Dead) ~  
                 Species + timeseries + (1|Location), 
               data = dat, family= "binomial")
dat$timeseries <- dat$Study_Date - 1998

plot(predict(glmm1, data.frame(Species="An. Arab", timeseries= dat$timeseries), 
             type="response", re.form = NA) ~
       dat$timeseries, frame=FALSE, bty="n", pch="", xaxt="n", ylim=c(0, 0.5), 
     ylab="Predicted mortality", xlab="Year", 
     main = "Predicted mortality by species", 
     cex.lab=1.6, yaxt="n")
axis(1, at=c(1:17), labels=1999:2015, cex.axis=1.8)
axis(2, las=2, at=seq(0, 1, 0.2), labels=seq(0, 1, 0.2), cex.axis=1.8) 

COLS <- c("blue", "red", "purple", "aquamarine3", "orange")
PCH <- c(17, 15, 19, 20, 5)

for(i in 1:length(unique(levels(dat$Species)))){
  points((dat$Total_Dead[dat$Species == levels(dat$Species)[i]] / 
            dat$Total_Collected[dat$Species == levels(dat$Species)[i]]) ~ 
           dat$timeseries[dat$Species == levels(dat$Species)[i]], 
         pch=PCH[i], col=COLS[i])
  lines(predict(glmm1, data.frame(Species=levels(dat$Species)[i], 
                                  timeseries = dat$timeseries), type="response", 
                re.form = NA) ~ dat$timeseries, lty=4, col=COLS[i])
}

bootstrap <- bootMer(x=glmm1, FUN= fixef, nsim=200)

Upvotes: 0

Views: 761

Answers (1)

Derek Corcoran
Derek Corcoran

Reputation: 4102

for some reason Bootmer has problems with that, you have to use the mertools package

library(merTools)
preds <- predictInterval(glmm1, newdata = your.datarame, n.sims = 1000)

I would use then then the preds data.frame to plot, the resulting data.frame has the fit, upper and lower limit, then you can use geom_ribbon to plot it, if you need more help let me know.

now bear with me, you actually want to make a new standardized dataset for your graph. If you use this code it will work:

glmm1 <- glmer(cbind(Total_Dead, Total_Collected - Total_Dead) ~  
             Species + timeseries + (1|Location),
           data = dat,family= "binomial")

fit your model, then create your new data set, this will have your timeseries from 1 to 16 for each species, in your first location (Akron), note that you will have to do this for each location if you want the graph for each location, you can do that just by changing the number between [] from 1, to 2 up to your 17 locations

new.data <-data.frame(timeseries = rep(1:16, times = 5), Species = rep(unique(dat$Species), each = 16), Location = rep(unique(dat$Location)[1], each = 80))

Then predict the values and intervals for such dataset

preds <- predictInterval(glmm1, newdata = new.data, n.sims = 1000)

now join this prediction to your new.data

new.data <- cbind(new.data, preds)

and finally plot it with different colors for each species

ggplot(new.data, aes(x = timeseries, y = fit)) + geom_ribbon(aes(ymax=upr, ymin =lwr, fill=Species), alpha = 0.9)+ geom_line(aes(color=Species))

If you don't understand something don't hesitate to keep asking, currently your standard errors are quite big so first check to see if you like that better

ggplot(new.data, aes(x = timeseries, y = fit)) + geom_line(aes(color=Species))

Upvotes: 1

Related Questions