R - boot.ci()返回引导样本的奇异置信区间

Question

我有两组正实数

> dput(group1)
c(2.10753, 2.57251, 2.61687, 4.62551, 7.13166, 6347.63, 4.22139, 
10.7373, 2.11568, 2.71866, 4.09376, 10.9046, 109807, 5.87156, 
3.17082, 3.4703, 2.47262, 9.24319, 34.6945, 5.72567, 12.0134, 
108.33, 6.60707, 6.24304, 3.59048, 10.3174, 48.0265, 5.32097, 
3.77157, 6.67401, 22.633, 34.8186, 21.5315, 9.42882, 7.10627, ...)
> dput(group2)
c(4.88474, 65.4318, 128.101, 24.1271, 5.44262, 54.8987, 2.85175, 
14.1089, 172.23, 66.8563, 6.74067, 2.19603, 2.12985, 4.12735,
16.401, 3.22688, 15.6943, 4.32861, 36.4752, 7.33769, 75.855, 
62.7653, 35.1786, 3.71099, 29.0186, 34.4472, 19.1061, 2.75174, ...)

Group1有~1000个值，group2有~ 30000个值。我对两组之间的中位数比率很感兴趣，并使用以下r函数来计算2000个引导样本中的每个样本的比率(关于boot()命令，请参阅下面的函数输出):

medianRatio <- function(x, i, noGroup1, noGroup2) {
    all <- x[i]
    currGroup1 <- all[1:noGroup1]
    currGroup2 <- all[c(noGroup1 + 1):length(all)]
    ratio <- median(currGroup1) / median(currGroup2)
    return(ratio)
}

boot()函数的调用如下所示

ORDINARY NONPARAMETRIC BOOTSTRAP

Call:
boot(data = c(group1, group2), statistic = medianRatio, R = noBs, 
    noGroup1 = length(group1), noGroup2 = length(group2))

Bootstrap Statistics :
    original      bias    std. error
t1*  1.08847 -0.08597889  0.05451763`

所得bootstrap样本分布的均值为1.002,sd为0.054(直方图的目视检查证实了1左右的正态分布)。

范围(Group1_Group2.BS $ t)[1]

然而，当我在引导对象上运行boot.ci()时，报告的置信区间是

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 2000 bootstrap replicates
CALL : 
boot.ci(boot.out = Group1_Group2.BS, type = "norm")
Intervals : 
Level      Normal        
95%   ( 1.068,  1.281 )  
Calculations and Intervals on Original Scale

我不明白这里发生了什么，因为报告的置信区间甚至没有覆盖bootstrap样本(对称)分布的平均值。我错过了什么?

Answer 1

了解正常置信区间是如何计算的可能有助于澄清问题。我在这个问题的答案中找到了一个很好的解释。

bootstrap正态置信区间围绕统计量的观测值建立，并进行偏差校正，以解决bootstrap统计量的中间值与观测值之间的差异。CI是通过使用原始样本中感兴趣的统计量的估计值，减去偏差，并加上1.96倍的自举标准误差来计算的。如果您"手工"处理boot对象的值，您将看到您获得的值与boot.ci相同。

1.08847 - -0.08597889 - 1.96*0.05451763
[1] 1.067594
1.08847 - -0.08597889 + 1.96*0.05451763
[1] 1.281303

您的偏差足够大，可以使百分位数CI和正常CI之间的差异明显。当我考虑自举置信区间时，我希望使用自举分布减去偏差的平均值，而不是观察到的统计量减去偏差。我没有花足够的时间来思考这个问题，我也没有带我的引导笔记，但是您可以更仔细地检查这个问题。

Answer 2

引导错误。您在函数medianRatio中假设x[i]的第一部分将保存组1的值，第二部分保存组2的值。这是不正确的，因为i只是一个从1替换为noGroup1+noGroup2的样本。您需要在引导命令中使用strata选项。我想下面的方法会有用。style ="f"选项指定引导将生成一个长度为noGroup1+noGroup2的向量，该向量表示为引导样本选择每个观测值的次数(0,1,2，…)。这段代码基于引导文档中的一个示例。

medRatio <- function(x, f, noGroup1){
 gp1 <- 1:noGroup1
 m1 <- median(rep(x[gp1],f[gp1]))
 m2 <- median(rep(x[-gp1],f[-gp1]))
 return(m1/m2)
}
n1 <- length(group1)
n2 <- length(group2)
boot(c(group1,group2),medRatio,R=2000,stype = "f",strata = rep(1:2,c(n1,n2)),noGroup1=n1)