r - R：在插入符号模型和 glm 模型上使用 rms:bootcov 计算引导估计

Question

如何使用bootcov包中的函数计算回归系数的引导估计rms？我使用示例数据集尝试了以下操作，但出现错误：

library(mlbench)
data(PimaIndiansDiabetes)

library(caret)
trControl <- trainControl(method = "repeatedcv",
                          repeats = 3,
                          classProbs = TRUE,
                          number = 10, 
                          savePredictions = TRUE,
                          summaryFunction = twoClassSummary)

caret_model <- train(diabetes~., 
                     data=PimaIndiansDiabetes, 
                     method="glm", 
                     trControl=trControl)

library(rms)
set.seed(1234)
reduced_model_bootcov <- bootcov(caret_model$finalModel, B=100)

错误是：

bootcov(caret_model$finalModel, B = 100) 中的错误：您没有在 fit 中指定 x=TRUE 和 y=TRUE

如果我使用该函数glm来构建模型，我会这样做：

model <- glm(diabetes~., 
             data=PimaIndiansDiabetes, 
             family=binomial,
             x=TRUE, y=TRUE)
model_bootcov <- bootcov(model, B=100)

但同样，我得到了一个不同的错误：

bootcov 中的错误（模型，B = 100）：钳工无效

标签： rlogistic-regressionstatistics-bootstrap

Answer 1

原来在 rms 中有一个称为 glm 的拟合函数，它是 glm 的包装器，但如果您对使用 bootcov 感兴趣，也可以使用它。所以要让 bootcov 工作：

library(mlbench)
library(rms)
data(PimaIndiansDiabetes)

model <- rms::Glm(diabetes~., 
             data=PimaIndiansDiabetes, 
             family=binomial,
             x=TRUE, y=TRUE)
model_bootcov <- bootcov(model, B=1000)

要使用引导：

library(boot)
glm.fun <- function(dat, inds){
  fit <- glm(diabetes~.,family=binomial,data=dat[inds,])
      coef(fit)
     }
model_boot <- boot(PimaIndiansDiabetes, glm.fun, R = 1000)

我们可以比较两个不同的模型是如何引导的，当然种子是不同的，很可能你需要先设置相似的种子：

library(tidyr)
library(dplyr)
library(ggplot2)

melt_matrix = function(mat,NAMES,X){
colnames(mat) = NAMES
data.frame(mat) %>% 
tibble::rownames_to_column("B") %>% 
pivot_longer(-B) %>%
mutate(type=X)
}

VAR = names(coef(model))

plotdf = rbind(
melt_matrix(model_boot$t,VAR,"boot"),
melt_matrix(model_bootcov$boot.Coef,VAR,"bootcov")
)

ggplot(plotdf,aes(x=type,y=value))+ geom_violin() + facet_wrap(~name,scale="free_y")