Bootstrapping CI для модели логистической регрессии

Question

У меня есть модель логистической регрессии, которую я использую, чтобы предсказать размер в зрелости для короля-краба, но у меня возникли проблемы с настройкой кода для начальной загрузки с использованием загрузочного пакета. Это то, что у меня есть:

#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") 

Я хотел бы, чтобы грузиться статистика CL = 98,97, так как я заинтересован в размере, при котором 50% краба являются зрелыми, но я понятия не имею, как настроить свою функцию чтобы указать эту статистику и не говоря уже о функции бутстрапа вообще, чтобы получить мой 95% ДИ Любая помощь будет принята с благодарностью! Благодаря!

Answer 1

В каждой итерации начальной загрузки, вы хотите сделать что-то вроде

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 

где указать диапазон значений CL для любой точности вам нужно. Затем вычислите прогнозируемую вероятность зрелости на каждом из этих уровней. Затем найдите пороговое значение в диапазоне, где прогнозируемая вероятность равна 0,5. Это статистика, которая звучит так, как будто вы хотите загрузиться.

Для этого вам также не нужен boot. Если определить функцию, образцы и выходы, статистика как его результат, вы можете просто сделать replicate(1000, myfun), чтобы получить вашу самозагрузку распределения следующим образом:

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