Bootstrapping CI for a Logistic Regression Model

Question

I have a logistic regression model that I am using to predict size at maturity for king crab, but I am having trouble setting up the code for bootstrapping using the boot package. This is what I have:

#FEMALE GKC SAM#
LowerChatham<-read.table(file=file.choose(),header=TRUE)

#LOGISTIC REGRESSION FIT#
glm.out<-glm(Mature~CL,family=binomial(link=logit),data=LowerChatham)
plot(Mature~CL,data=LowerChatham)
lines(LowerChatham$CL,glm.out$fitted,col="red")
title(main="Lower Chatham")
summary(glm.out)
segments(98.9,0,98.9,0.5,col=1,lty=3,lwd=3)

SAM<-data.frame(CL=98.97)
predict(glm.out,SAM,type="response")

I would like to to bootstrap the statistic CL=98.97 since I am interested in the size at which 50% of crab are mature, but I have no idea how to setup my function to specify the that statistic and let alone the bootstrap function in general to get my 95% C.I. Any help would be greatly appreciated! Thanks!

Answer 1

In each bootstrap iteration, you want to do something like

range <- 1:100 # this could be any substantively meaningful range
p <- predict(glm.out, newdata = data.frame(CL=range), "response")
range[match(TRUE,p>.5)] # predicted probability of 50% maturity

where you specify a range of values of CL to whatever precision you need. Then calculate the predicted probability of maturity at each of those levels. Then find the threshold value in the range where the predicted probability cross 0.5. This is the statistic it sounds like you want to bootstrap.

You also don't need the boot to do this. If you define a function that samples and outputs that statistic as its result, you can just do replicate(1000, myfun) to get your bootstrap distribution, as follows:

myfun <- function(){
    srows <- sample(1:nrow(LowerChatham),nrow(LowerChatham),TRUE)
    glm.out <- (Mature ~ CL, family=binomial(link=logit), data=LowerChatham[srows,])
    range <- 1:100 # this could be any substantively meaningful range
    p <- predict(glm.out, newdata = data.frame(CL=range), "response")
    return(range[match(TRUE,p>.5)]) # predicted probability of 50% maturity
}
bootdist <- replicate(1000, myfun()) # your distribution
quantile(unlist(bootdist),c(.025,.975)) # 95% CI