使用带有 R 的引导程序计算 p 值

机器算法验证 r 假设检验 p 值 引导程序 置换检验
2022-01-20 18:59:11

我使用“引导”包来计算近似的2 面引导 p 值,但结果与使用 t.test 的 p 值相差太远。我无法弄清楚我在 R 代码中做错了什么。有人可以给我一个提示吗

time = c(14,18,11,13,18,17,21,9,16,17,14,15,
         12,12,14,13,6,18,14,16,10,7,15,10)
group=c(rep(1:2, each=12))
sleep = data.frame(time, group)

require(boot)
diff = function(d1,i){
    d = d1[i,]
    Mean= tapply(X=d$time, INDEX=d$group, mean)
    Diff = Mean[1]-Mean[2]
    Diff
}

set.seed(1234)
b3 = boot(data = sleep, statistic = diff, R = 5000, strata=sleep$group)

pvalue = mean(abs(b3$t) > abs(b3$t0))
pvalue 

2 面自举 p 值 (pvalue) = 0.4804,但 t.test 的 2 面 p 值为 0.04342。两个 p 值大约相差 11 倍。这怎么可能发生?

3个回答

您正在使用 bootstrap 在观察到的数据的经验分布下生成数据。这对于给出两个均值之间差异的置信区间很有用:

> quantile(b3$t,c(0.025,0.975))
     2.5%     97.5% 
0.4166667 5.5833333 

得到一个p-value,您需要在原假设下生成排列。这可以这样做,例如:

diff2 = function(d1,i){
    d = d1; 
    d$group <- d$group[i];  # randomly re-assign groups
    Mean= tapply(X=d$time, INDEX=d$group, mean)
    Diff = Mean[1]-Mean[2]
    Diff
}

> set.seed(1234)
> b4 = boot(data = sleep, statistic = diff2, R = 5000)
> mean(abs(b4$t) > abs(b4$t0))
[1] 0.046

在此解决方案中,组的大小不是固定的,您通过从初始组集中引导来随机重新分配一个组给每个人。对我来说这似乎是合法的,但是更经典的解决方案是固定每个组的个体数量,因此您只需置换组而不是引导(这通常是由实验设计驱动的,其中组大小是事先固定的):

> R <- 10000; d <- sleep
> b5 <- numeric(R); for(i in 1:R) { 
+    d$group <- sample(d$group, length(d$group)); 
+    b5[i] <- mean(d$time[d$group==1])-mean(d$time[d$group==2]); 
+ }
> mean(abs(b5) > 3)
[1] 0.0372

猫王的答案依赖于排列,但在我看来,它并不清楚原始引导方法有什么问题。让我讨论一个完全基于引导程序的解决方案。

原始模拟的关键问题是引导程序始终为您提供测试统计量的真实分布。但是,在计算 p 值时,您必须将获得的检验统计值与其在 H0 下的分布进行比较,即不与真实分布进行比较!

[让我们说清楚。例如,已知经典 t 检验的检验统计量 T 在 H0 下具有经典的“中心”t 分布和一般的非中心分布。然而,每个人都熟悉这样一个事实,即 T 的观测值与经典的“中心”t 分布进行比较,即人们不会试图获得真正的 [非中心] t 分布来与 T 进行比较。]

您的 p 值 0.4804 太大了,因为检验统计量 Mean[1]-Mean[2] 的观察值“t0”非常靠近自举样本“t”的中心。这是自然的并且通常总是如此[即,不管 H0 的有效性如何],因为自举样本“t”模拟了 Mean[1]-Mean[2] 的实际分布。但是,如上所述 [以及 Elvis],您真正需要的是 H0 下的 Mean[1]-Mean[2] 分布。很明显,

1) 在 H0 下,Mean[1]-Mean[2] 的分布将以 0 为中心,

2)它的形状不依赖于H0的有效性。

这两点意味着 H0 下 Mean[1]-Mean[2] 的分布可以通过自举样本“t”SHIFTED 来模拟,使其以 0 为中心。在 R 中:

b3.under.H0 <- b3$t - mean(b3$t)

相应的 p 值将是:

mean(abs(b3.under.H0) > abs(b3$t0))

这给了你一个“非常好的”值 0.0232。:-)

让我注意,上面提到的点“2)”被称为检验统计量的“平移等方差”,它通常不必成立!即对于某些测试统计,自举“t”的移动不会为您提供对 HO 下测试统计分布的有效估计!看看这个讨论,尤其是 P. Dalgaard 的回复:http: //tolstoy.newcastle.edu.au/R/e6/help/09/04/11096.html

您的测试问题确实会产生测试统计量的完美对称分布,但请记住,在测试统计量的自举分布偏斜的情况下,获取 TWO-SIDED p 值存在一些问题。再次,阅读上面的链接。

[最后,我会在你的情况下使用“纯”置换测试;即猫王回答的后半部分。:-)]

有多种计算引导 CI 和 p 值的方法。主要问题是引导程序不可能在零假设下生成数据。置换测试是一种可行的基于重采样的替代方法。要使用适当的引导程序,您必须对检验统计量的抽样分布做出一些假设。

关于缺乏测试不变性的评论:完全有可能找到 95% 的 CI,不包括 null 但 ap > 0.05,反之亦然。为了获得更好的一致性,null 下 bootstrap 样本的计算必须β0=β^β^而不是β0=β^β^. 也就是说,如果密度在 bootstrap 样本中向右倾斜,则密度必须在空值中向左倾斜。使用诸如此类的非分析(例如重采样)解决方案来反转 CI 的测试是不可能的。

正常引导

一种方法是正态自举,您获取自举分布的均值和标准差,通过移动分布并使用原始自举样本中估计点处的零分布的正态百分位数来计算零下的采样分布. 当 bootstrap 分布正常时,这是一种合理的方法,目视检查通常就足够了。使用这种方法的结果通常非常接近稳健的或基于三明治的误差估计,这种估计对异方差和/或有限样本方差假设具有稳健性。正态检验统计量的假设是我将讨论的下一个自举检验中假设的更强条件。

百分位引导

另一种方法是百分位自举法,我认为我们大多数人在谈到自举法时都会考虑这一点。在这里,参数的自举分布估计了备择假设下样本的经验分布。这种分布可能是非正态的。采用经验分位数很容易计算出 95% CI。但是一个重要的假设是这样的分布是关键的。这意味着如果基础参数发生变化,分布的形状只会移动一个常数,而尺度不一定会改变。这是一个强有力的假设!如果这成立,您可以生成“在零假设下的统计分布”(DSNH 或F0) 通过从估计值中减去自举分布,然后计算 DSNH 的百分比比您的估计值“更极端”2×min(F0(β^),1F0(β^))

学生化引导

计算值的最简单的引导解决方案是使用学生化引导。在每次 bootstrap 迭代中,计算统计量及其标准误差并返回学生统计量。这给出了假设的自举学生分布,可以很容易地用于计算 cis 和 p 值。这也是偏置校正加速引导背后的直觉的基础。t 分布在零值下更容易移动,因为异常结果被其相应的高方差降低了权重。p

编程示例

例如,我将使用citybootstrap 包中的数据。使用以下代码计算引导置信区间:

ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)
city.boot <- boot(city, ratio, R = 999, stype = "w", sim = "ordinary")
boot.ci(city.boot, conf = c(0.90, 0.95),
        type = c("norm", "basic", "perc", "bca"))

并产生这个输出:

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates

CALL : 
boot.ci(boot.out = city.boot, conf = c(0.9, 0.95), type = c("norm", 
    "basic", "perc", "bca"))

Intervals : 
Level      Normal              Basic         
90%   ( 1.111,  1.837 )   ( 1.030,  1.750 )   
95%   ( 1.042,  1.906 )   ( 0.895,  1.790 )  

Level     Percentile            BCa          
90%   ( 1.291,  2.011 )   ( 1.292,  2.023 )   
95%   ( 1.251,  2.146 )   ( 1.255,  2.155 )  
Calculations and Intervals on Original Scale

正常 bootstrap 的 95% CI 是通过计算获得的:

with(city.boot, 2*t0 - mean(t) + qnorm(c(0.025, 0.975)) %o% sqrt(var(t)[1,1]))

因此得到 p 值:

> with(city.boot, pnorm(abs((2*t0 - mean(t) - 1) / sqrt(var(t)[1,1])), lower.tail=F)*2)
[1] 0.0315

这同意 95% 的正常 CI 不包括 1 的空值比值。

获得百分位数 CI(由于关系的方法而存在一些差异):

quantile(city.boot$t, c(0.025, 0.975))

百分位引导的 p 值为:

cvs <- quantile(city.boot$t0 - city.boot$t + 1, c(0.025, 0.975))
mean(city.boot$t > cvs[1] & city.boot$t < cvs[2])

给出 0.035 的 ap,这也与从值中排除 1 的置信区间一致。我们一般不能观察到,虽然百分位数 CI 的宽度几乎与正常 CI 一样宽,并且百分位数 CI 离零值更远,百分位数 CI 应该提供较低的 p 值。这是因为百分位数方法的 CI 下的抽样分布形状是非正态的。