我应该选择哪种自举回归模型？

机器算法验证 回归 物流 多重回归 引导程序
2022-03-13 07:37:38

我有一个二元逻辑回归模型，其中包含 DV（疾病：是/否）和 5 个预测变量（人口统计 [年龄、性别、吸烟（是/否）]、医学指数（序数）和一种随机治疗 [是/否]）。我还对所有双边交互项进行了建模。主要变量居中，没有多重共线性迹象（所有 VIF < 2.5）。

我有一些疑问：

  1. 自举比我的单一模型有优势吗？如果是这样的话，

  2. 我应该选择哪种自举模型？我只是想看看自举算法是否遵循随机方法来创建新样本，或者它们是否具有严格的算法。因此，我在每次尝试中重新采样 1000 次（所以我有几个自举模型，每个模型都有 1000 次试验）。然而，每次自举模型的系数都不同（尽管试验次数始终为 1000）。所以我想知道我应该为我的报告选择哪一个？有些变化很小，不会影响我的系数的显着性，但有些变化使我的一些系数不显着（例如，只有原始模型中 P 值接近 0.05 的那些变化为 0.06）。

  3. 我应该选择更高的数字，例如 10,000 吗？我如何确定这个限制？

  4. 我应该再次引导吗？如果它的结果​​每次都不同，我可以依赖它的结果吗？

  5. 你有什么其他想法可以帮助我处理我的案子吗？

非常感谢。

1个回答

Bootstrapping是一种重采样方法，用于估计回归系数的抽样分布，从而计算回归系数的标准误差/置信区间。这篇文章有一个很好的解释。有关您需要多少次复制的讨论，请参阅这篇文章。

  1. 非参数引导程序反复重新采样并随机抽取您的观测值并替换（即，有些观测值只绘制一次，有些则多次绘制，有些则根本不绘制），然后计算逻辑回归并存储系数。这是重复的n次。所以你最终会得到 10'000 个不同的回归系数。然后可以使用这 10'000 个系数来计算它们的置信度。由于使用了伪随机数生成器，您可以将种子设置为任意数字，以确保每次都获得完全相同的结果（参见下面的示例）。要真正获得稳定的估计，我建议进行 1000 多次重复，也许 10'000。您可以多次运行引导程序，看看无论您进行 1000 次还是 10'000 次复制，估计值是否会发生很大变化。简而言之：您应该进行复制，直到达到收敛。如果您的 bootstrap 估计值在您的估计值和观察到的单一模型之间存在差异，这可能表明观察到的模型不能恰当地反映您的样本结构。bootR，例如，提出“偏差”，即单个模型的回归系数与引导样本的平均值之间的差异。
  2. 执行 bootstrap 时，您对单个 bootstrap 样本不感兴趣，而是对 10'000 个 bootstrap 样本上的统计数据（例如回归系数）分布感兴趣。
  3. 我会说 10'000 比 1000 好。对于现代计算机，这不应该构成问题。在下面的示例中，我的 PC 花了大约 45 秒来绘制 10'000 个样本。当然，这取决于您的样本量。样本量越大，迭代次数就应该越高，以确保考虑到每个观察结果。
  4. 你是什​​么意思“结果每次都不同”？回想一下，在每个引导步骤中，观察都是新绘制的并带有替换。因此，由于您的观察结果不同，您最终可能会得到稍微不同的回归系数。但正如我所说：您对单个引导样本的结果并不真正感兴趣。当您的复制次数足够高时，引导程序每次都会产生非常相似的置信区间和点估计。

这是一个示例R：

library(boot)

mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")

head(mydata)

mydata$rank <- factor(mydata$rank)

my.mod <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial")

summary(my.mod)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.989979   1.139951  -3.500 0.000465 ***
gre          0.002264   0.001094   2.070 0.038465 *  
gpa          0.804038   0.331819   2.423 0.015388 *  
rank2       -0.675443   0.316490  -2.134 0.032829 *  
rank3       -1.340204   0.345306  -3.881 0.000104 ***
rank4       -1.551464   0.417832  -3.713 0.000205 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Set up the non-parametric bootstrap

logit.bootstrap <- function(data, indices) {

  d <- data[indices, ]
  fit <- glm(admit ~ gre + gpa + rank, data = d, family = "binomial")

  return(coef(fit))
}

set.seed(12345) # seed for the RNG to ensure that you get exactly the same results as here

logit.boot <- boot(data=mydata, statistic=logit.bootstrap, R=10000) # 10'000 samples

logit.boot

Bootstrap Statistics :
        original        bias    std. error
t1* -3.989979073 -7.217244e-02 1.165573039
t2*  0.002264426  4.054579e-05 0.001146039
t3*  0.804037549  1.440693e-02 0.354361032
t4* -0.675442928 -8.845389e-03 0.329099277
t5* -1.340203916 -1.977054e-02 0.359502576
t6* -1.551463677 -4.720579e-02 0.444998099

# Calculate confidence intervals (Bias corrected ="bca") for each coefficient

boot.ci(logit.boot, type="bca", index=1) # intercept
95%   (-6.292, -1.738 )  
boot.ci(logit.boot, type="bca", index=2) # gre
95%   ( 0.0000,  0.0045 ) 
boot.ci(logit.boot, type="bca", index=3) # gpa
95%   ( 0.1017,  1.4932 )
boot.ci(logit.boot, type="bca", index=4) # rank2
95%   (-1.3170, -0.0369 )
boot.ci(logit.boot, type="bca", index=5) # rank3
95%   (-2.040, -0.629 )
boot.ci(logit.boot, type="bca", index=6) # rank4
95%   (-2.425, -0.698 )

引导输出显示原始回归系数（“原始”）及其偏差，即原始系数与引导系数之间的差异。它还给出了标准错误。请注意，它们比原始标准误差大一点。

从置信区间来看，偏差校正（“bca”）通常是首选。它给出了原始尺度的置信区间。对于优势比的置信区间，只需将置信限取幂即可。

其它你可能感兴趣的问题
上一篇如何报告整个人口的数据？ 下一篇检测歌曲的各个部分