2

我有一个逻辑回归模型,我用它来预测帝王蟹成熟时的大小,但是我在使用引导包设置引导代码时遇到了麻烦。这就是我所拥有的:

#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% CI 任何帮助将不胜感激!谢谢!

4

1 回答 1

3

在每次引导迭代中,您都想做类似的事情

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
于 2013-07-03T06:17:43.617 回答