逻辑回归中的重采样

Resampling in logistic regression

我有一个简单的数据集,其中包含一个 Y 和 10 个预测变量 (X1-X10),针对 100 个观测值编码为 0,1 或 2。

 n <- 100
 Y <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
 X1 <- sample(x=c(0,1,2), size=n, replace=TRUE, prob=c(0.1,0.4,0.5))
 X2 <- sample(x=c(0,1,2), size=n, replace=TRUE, prob=c(0.5,0.25,0.25))
 X3 <- sample(x=c(0,1,2), size=n, replace=TRUE, prob=c(0.3,0.4,0.4))
 X4 <- sample(x=c(0,1,2), size=n, replace=TRUE, prob=c(0.35,0.35,0.3))
 X5 <- sample(x=c(0,1,2), size=n, replace=TRUE, prob=c(0.1,0.2,0.7))
 X6 <- sample(x=c(0,1,2), size=n, replace=TRUE, prob=c(0.8,0.1,0.1))
 X7 <- sample(x=c(0,1,2), size=n, replace=TRUE, prob=c(0.1,0.1,0.8))
 X8 <- sample(x=c(0,1,2), size=n, replace=TRUE, prob=c(0.35,0.35,0.3))
 X9 <- sample(x=c(0,1,2), size=n, replace=TRUE, prob=c(0.35,0.35,0.3))
X10 <- c(0,2,2,2,2,2,2,2,0,2,0,2,2,0,0,0,0,0,2,0,0,2,2,0,0,2,2,2,0,2,0,2,0,2,1,2,1,1,1,1,1,1,1,1,1,1,1,0,1,2,2,2,2,2,2,2,2,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,1,0,0,0,0)

datasim <- data.frame(Y,X1,X2,X3,X4,X5,X6,X7,X8,X9,X10)

我正在尝试按如下方式进行 bootstrap 重采样,它可以为一个变量生成 100 组不同的样本。

 B <- 100
 n <- length(datasim$X1)
 boot.samples <- matrix(sample(datasim$X1, size=B*n, replace=TRUE),B,n)

现在,我正在尝试合并一个函数来使用 GLM 计算偏差差异。我的愿望是为每个 bootstrap 样本(100 个值)生成 dDeviance。我尝试了以下函数,但它只给了我 100 个相似的 dDeviance 值。

 xfunction <- function(x){
 glmfit <- glm(Y~X1, family="binomial", data=datasim)
 dDeviance <- glmfit$null.deviance-glmfit$deviance
 return(dDeviance)
 }

 boot.statistics <- apply(boot.samples,1,xfunction)

正如 Jeffrey 所说,数据应该 = x。

 xfunction <- function(x){
   glmfit <- glm(Y~X1, family="binomial", data=x)
   dDeviance <- glmfit$null.deviance-glmfit$deviance
   return(dDeviance)
 }

 boot.statistics <- apply(boot.samples,1,xfunction)

xfunction 的参数在这样的 apply 中使用时是矩阵中的一行。在您的原始代码中,该行未被使用,您每次都是 运行 相同数据的函数。解决此类问题的一种方法是按照建议(glmfit <- glm(Y~X1, family="binomial", data=x))每次将 glm 中的数据参数更改为新数据,但这假设 x 将是一个数据框,其中列名为 Y 和 X1 ,而您实际拥有的 x 是 X1 的值向量。最简单的解决方案是在每次拟合中更改 X1。

xfunction <- function(x){
  glmfit <- glm(Y~x, family="binomial")
  dDeviance <- glmfit$null.deviance-glmfit$deviance
  return(dDeviance)
}

boot.statistics <- apply(boot.samples,1,xfunction)