如何在 R 中执行测试的现场重要性
How to perform field significance of a test in R
I 运行 对 168 个点进行假设检验，以查看 A 阶段的观察值是否与 B 阶段有显着差异，并得出结论，其中 30% 的点显示两个阶段之间存在显着差异。如何检验这些结果的现场显着性？是bootstrapping方法吗？那么在这个数据中如何运行 bootstrap。 dtt 具有每个假设检验和结果的 p 值。
dtt<-dtt<-structure(list(stn = c("Stn_1", "Stn_2", "Stn_3", "Stn_4", "Stn_5",
"Stn_6", "Stn_7", "Stn_8", "Stn_9", "Stn_10", "Stn_11", "Stn_12",
"Stn_13", "Stn_14", "Stn_15", "Stn_16", "Stn_17", "Stn_18", "Stn_19",
"Stn_20", "Stn_21", "Stn_22", "Stn_23", "Stn_24", "Stn_25", "Stn_26",
"Stn_27", "Stn_28", "Stn_29", "Stn_30", "Stn_31", "Stn_32", "Stn_33",
"Stn_34", "Stn_35", "Stn_36", "Stn_37", "Stn_38", "Stn_39", "Stn_40",
"Stn_41", "Stn_42", "Stn_43", "Stn_44", "Stn_45", "Stn_46", "Stn_47",
"Stn_48", "Stn_49", "Stn_50", "Stn_51", "Stn_52", "Stn_53", "Stn_54",
"Stn_55", "Stn_56", "Stn_57", "Stn_58", "Stn_59", "Stn_60", "Stn_61",
"Stn_62", "Stn_63", "Stn_64", "Stn_65", "Stn_66", "Stn_67", "Stn_68",
"Stn_69", "Stn_70", "Stn_71", "Stn_72", "Stn_73", "Stn_74", "Stn_75",
"Stn_76", "Stn_77", "Stn_78", "Stn_79", "Stn_80", "Stn_81", "Stn_82",
"Stn_83", "Stn_84", "Stn_85", "Stn_86", "Stn_87", "Stn_88", "Stn_89",
"Stn_90", "Stn_91", "Stn_92", "Stn_93", "Stn_94", "Stn_95", "Stn_96",
"Stn_97", "Stn_98", "Stn_99", "Stn_100", "Stn_101", "Stn_102",
"Stn_103", "Stn_104", "Stn_105", "Stn_106", "Stn_107", "Stn_108",
"Stn_109", "Stn_110", "Stn_111", "Stn_112", "Stn_113", "Stn_114",
"Stn_115", "Stn_116", "Stn_117", "Stn_118", "Stn_119", "Stn_120",
"Stn_121", "Stn_122", "Stn_123", "Stn_124", "Stn_125", "Stn_126",
"Stn_127", "Stn_128", "Stn_129", "Stn_130", "Stn_131", "Stn_132",
"Stn_133", "Stn_134", "Stn_135", "Stn_136", "Stn_137", "Stn_138",
"Stn_139", "Stn_140", "Stn_141", "Stn_142", "Stn_143", "Stn_144",
"Stn_145", "Stn_146", "Stn_147", "Stn_148", "Stn_149", "Stn_150",
"Stn_151", "Stn_152", "Stn_153", "Stn_154", "Stn_155", "Stn_156",
"Stn_157", "Stn_158", "Stn_159", "Stn_160", "Stn_161", "Stn_162",
"Stn_163", "Stn_164", "Stn_165", "Stn_166", "Stn_167", "Stn_168"
), pval = c(0.205944631, 0.63991585, 0.473120067, 0.34875961,
0.292140039, 0.326105934, 0.529800338, 0.294475321, 0.141110971,
0.368350989, 0.552273175, 0.643845842, 0.07104491, 0.002432443,
0.003331365, 0.116333091, 0.585496713, 0.227960311, 0.172988608,
0.142913486, 0.001836251, 0.002553918, 0.066330084, 0.048866324,
0.507511564, 0.304430083, 0.367805688, 0.181954789, 0.318772861,
0.199522509, 0.002678304, 0.04779772, 0.017131339, 0.031137852,
NA, 0.26161318, 0.586668965, 0.0043644, 0.098939189, 0.028705313,
0.041562653, 0.09003053, 0.157823558, 0.161172547, 0.474951712,
0.136885745, NA, NA, NA, 0.050304544, NA, NA, NA, NA, 0.009360088,
0.126128118, 0.112494159, 0.220780636, 0.133918794, 0.547804304,
0.035639161, 0.099166469, 0.024599266, 0.063829305, 0.051450678,
0.094083816, 0.025413468, 0.041048015, 0.032694959, 0.017755539,
0.104045842, 0.005085752, 0.043865633, 0.254403589, 0.06702142,
0.750230985, 0.11802067, 0.086793641, 0.43275653, 0.249168613,
0.675590582, 0.146278867, 0.470232808, 0.560136445, 0.447567809,
0.790315815, 0.027195565, 0.304420281, 0.231429562, 0.444931845,
0.169611396, 0.964873224, 0.995751223, 0.921615572, 0.972023663,
0.779856105, 0.211587188, 0.021518622, 0.026315789, 0.07704875,
0.535102383, 0.342347191, 0.092329512, 0.566113195, 0.042437174,
0.038141264, 0.220853181, 0.048037069, 0.269139259, 0.056426699,
0.048866324, 0.552336171, 0.432409847, 0.248245841, 0.053496696,
0.294284245, 0.942909112, 0.519258899, 0.008882952, 0.055022687,
0.563962346, 0.159175889, 0.729845947, 0.610946838, 0.063839729,
0.542889465, 0.154845574, 0.396166658, 0.512275225, 0.636128255,
0.000265, 0.016868837, 0.041575921, 0.004828325, 0.170143423,
0.057359331, 0.463224782, 0.55416784, 0.286603865, 0.415871141,
0.473120067, 0.464303708, 0.011002312, 0.5, 0.348898333, 0.220983939,
0.009418115, 0.906935504, 0.601641288, 0.480741101, 0.79015601,
0.085766864, 0.016978584, NA, NA, 0.399770285, 0.927749655, 0.029916602,
0.119047989, 0.568590421, 0.540283117, 0.654617786, 0.917625002,
0.002563335, 0.010137604, 0.043544506, 0.085651037, 0.072239359
), hypo = c("no", "no", "no", "no", "no", "no", "no", "no", "no",
"no", "no", "no", "no", "yes", "yes", "no", "no", "no", "no",
"no", "yes", "yes", "no", "yes", "no", "no", "no", "no", "no",
"no", "yes", "yes", "yes", "yes", "yes", "no", "no", "yes", "no",
"yes", "yes", "no", "no", "no", "no", "no", "yes", "yes", "yes",
"no", "yes", "yes", "yes", "yes", "yes", "no", "no", "no", "no",
"no", "yes", "no", "yes", "no", "no", "no", "yes", "yes", "yes",
"yes", "no", "yes", "yes", "no", "no", "no", "no", "no", "no",
"no", "no", "no", "no", "no", "no", "no", "yes", "no", "no",
"no", "no", "no", "no", "no", "no", "no", "no", "yes", "yes",
"no", "no", "no", "no", "no", "yes", "yes", "no", "yes", "no",
"no", "yes", "no", "no", "no", "no", "no", "no", "no", "yes",
"no", "no", "no", "no", "no", "no", "no", "no", "no", "no", "no",
"yes", "yes", "yes", "yes", "no", "no", "no", "no", "no", "no",
"no", "no", "yes", "no", "no", "no", "yes", "no", "no", "no",
"no", "no", "yes", "yes", "yes", "no", "no", "yes", "no", "no",
"no", "no", "no", "yes", "yes", "yes", "no", "no")), row.names = c(NA,
-168L), spec = structure(list(cols = list(stnn = structure(list(), class = c("collector_character",
"collector")), wscore = structure(list(), class = c("collector_double",
"collector")), lat = structure(list(), class = c("collector_double",
"collector")), lon = structure(list(), class = c("collector_double",
"collector")), pval = structure(list(), class = c("collector_double",
"collector")), MWU_res = structure(list(), class = c("collector_character",
"collector")), hydro_zone = structure(list(), class = c("collector_character",
"collector"))), default = structure(list(), class = c("collector_guess",
"collector"))), class = "col_spec"), class = c("tbl_df", "tbl",
"data.frame"))
我读到您对 "phase A" 和 "phase B" 之间的差异进行了 168 次显着性检验。一个鲜为人知的事实是，零假设下的 p 值在 [0-1] 范围内均匀分布。因此，如果阶段之间没有差异，您会期望 "significant results"，即 p 值 <= 0.05，正好是 5% 的时间（给出或采用统计不精确性）。您清楚地看到 "non-null" 结果分布。
首先，我模拟了一些可能看起来像您的原始数据的东西来执行 t.test，前 40 个是从 A 阶段和 B 阶段之间具有不同均值的正态分布模拟的，最后 128 个来自相同的分布。
set.seed(100)
dat1 = lapply(1:40,function(i){
data.frame(
stn = paste0("Stn_",i),
value = c(rnorm(15,0,1),rnorm(15,2,1)),
type = rep(c("PhaseA","PhaseB"),each=15)
)
})
dat2 = lapply(41:168,function(i){
data.frame(
stn = paste0("Stn_",i),
value = rnorm(30),
type = rep(c("PhaseA","PhaseB"),each=15)
)
})
dat = do.call(rbind,c(dat1,dat2))
现在我们可以像您一样在 phaseA 和 phaseB 之间像您一样执行 t.test，分别针对每个站：
library(dplyr)
library(broom)
library(purrr)
obs_result = dat %>% group_by(stn) %>% do(tidy(t.test(value~type,data=.)))
如果我们看到哪个站的 p < 0.05:
which(obs_result$p.value<0.05)
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
[20] 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
[39] 39 40 52 57 67 90 113 124 166
所以前 40 个（真正不同的）的大部分测试显示 p < 0.05。但是您确实也有一些来自没有显示差异的电台的点击。因此，您需要通过执行以下操作来控制多次测试：
which(p.adjust(obs_result$p.value,"BH")<0.1)
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
[20] 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
[39] 39 40 57 90 124 166
所以你有 44 次命中.. 其中 4 次实际上是错误的 positives.This 没问题，因为我们允许 FDR 为 0.1。
我们还可以模拟您在 null 下的预期。在您的情况下，我们只是交换站点的标签并再次执行 t.test。这一次，您的 p 值应该与您观察到的非常不同：
permutated_pvalue = function(i,dat){
set.seed(i)
dat %>% group_by(stn) %>%
mutate(type=sample(type)) %>%
do(tidy(t.test(value~type,data=.))) %>%
pull(p.value)
}
par(mfrow=c(1,2))
hist(permutated_pvalue(0,dat),main="permutated pvalue")
hist(obs_result$p.value,main="observed pvalue")
最后，做一个 qqplot 来查看这些 p 值分布有何不同：
# 10 simulations to get null p-values
simulated_pvalues = 1:10 %>% map(perm,dat=dat) %>% unlist()
qqplot(-log10(simulated_pvalues),-log10(obs_result$p.value))
abline(a=0,b=1)
所以我建议你做这些检查，如果图表确实显示 p 值分布不同，那么就有一些全局差异的说法..
I 运行 对 168 个点进行假设检验，以查看 A 阶段的观察值是否与 B 阶段有显着差异，并得出结论，其中 30% 的点显示两个阶段之间存在显着差异。如何检验这些结果的现场显着性？是bootstrapping方法吗？那么在这个数据中如何运行 bootstrap。 dtt 具有每个假设检验和结果的 p 值。
dtt<-dtt<-structure(list(stn = c("Stn_1", "Stn_2", "Stn_3", "Stn_4", "Stn_5",
"Stn_6", "Stn_7", "Stn_8", "Stn_9", "Stn_10", "Stn_11", "Stn_12",
"Stn_13", "Stn_14", "Stn_15", "Stn_16", "Stn_17", "Stn_18", "Stn_19",
"Stn_20", "Stn_21", "Stn_22", "Stn_23", "Stn_24", "Stn_25", "Stn_26",
"Stn_27", "Stn_28", "Stn_29", "Stn_30", "Stn_31", "Stn_32", "Stn_33",
"Stn_34", "Stn_35", "Stn_36", "Stn_37", "Stn_38", "Stn_39", "Stn_40",
"Stn_41", "Stn_42", "Stn_43", "Stn_44", "Stn_45", "Stn_46", "Stn_47",
"Stn_48", "Stn_49", "Stn_50", "Stn_51", "Stn_52", "Stn_53", "Stn_54",
"Stn_55", "Stn_56", "Stn_57", "Stn_58", "Stn_59", "Stn_60", "Stn_61",
"Stn_62", "Stn_63", "Stn_64", "Stn_65", "Stn_66", "Stn_67", "Stn_68",
"Stn_69", "Stn_70", "Stn_71", "Stn_72", "Stn_73", "Stn_74", "Stn_75",
"Stn_76", "Stn_77", "Stn_78", "Stn_79", "Stn_80", "Stn_81", "Stn_82",
"Stn_83", "Stn_84", "Stn_85", "Stn_86", "Stn_87", "Stn_88", "Stn_89",
"Stn_90", "Stn_91", "Stn_92", "Stn_93", "Stn_94", "Stn_95", "Stn_96",
"Stn_97", "Stn_98", "Stn_99", "Stn_100", "Stn_101", "Stn_102",
"Stn_103", "Stn_104", "Stn_105", "Stn_106", "Stn_107", "Stn_108",
"Stn_109", "Stn_110", "Stn_111", "Stn_112", "Stn_113", "Stn_114",
"Stn_115", "Stn_116", "Stn_117", "Stn_118", "Stn_119", "Stn_120",
"Stn_121", "Stn_122", "Stn_123", "Stn_124", "Stn_125", "Stn_126",
"Stn_127", "Stn_128", "Stn_129", "Stn_130", "Stn_131", "Stn_132",
"Stn_133", "Stn_134", "Stn_135", "Stn_136", "Stn_137", "Stn_138",
"Stn_139", "Stn_140", "Stn_141", "Stn_142", "Stn_143", "Stn_144",
"Stn_145", "Stn_146", "Stn_147", "Stn_148", "Stn_149", "Stn_150",
"Stn_151", "Stn_152", "Stn_153", "Stn_154", "Stn_155", "Stn_156",
"Stn_157", "Stn_158", "Stn_159", "Stn_160", "Stn_161", "Stn_162",
"Stn_163", "Stn_164", "Stn_165", "Stn_166", "Stn_167", "Stn_168"
), pval = c(0.205944631, 0.63991585, 0.473120067, 0.34875961,
0.292140039, 0.326105934, 0.529800338, 0.294475321, 0.141110971,
0.368350989, 0.552273175, 0.643845842, 0.07104491, 0.002432443,
0.003331365, 0.116333091, 0.585496713, 0.227960311, 0.172988608,
0.142913486, 0.001836251, 0.002553918, 0.066330084, 0.048866324,
0.507511564, 0.304430083, 0.367805688, 0.181954789, 0.318772861,
0.199522509, 0.002678304, 0.04779772, 0.017131339, 0.031137852,
NA, 0.26161318, 0.586668965, 0.0043644, 0.098939189, 0.028705313,
0.041562653, 0.09003053, 0.157823558, 0.161172547, 0.474951712,
0.136885745, NA, NA, NA, 0.050304544, NA, NA, NA, NA, 0.009360088,
0.126128118, 0.112494159, 0.220780636, 0.133918794, 0.547804304,
0.035639161, 0.099166469, 0.024599266, 0.063829305, 0.051450678,
0.094083816, 0.025413468, 0.041048015, 0.032694959, 0.017755539,
0.104045842, 0.005085752, 0.043865633, 0.254403589, 0.06702142,
0.750230985, 0.11802067, 0.086793641, 0.43275653, 0.249168613,
0.675590582, 0.146278867, 0.470232808, 0.560136445, 0.447567809,
0.790315815, 0.027195565, 0.304420281, 0.231429562, 0.444931845,
0.169611396, 0.964873224, 0.995751223, 0.921615572, 0.972023663,
0.779856105, 0.211587188, 0.021518622, 0.026315789, 0.07704875,
0.535102383, 0.342347191, 0.092329512, 0.566113195, 0.042437174,
0.038141264, 0.220853181, 0.048037069, 0.269139259, 0.056426699,
0.048866324, 0.552336171, 0.432409847, 0.248245841, 0.053496696,
0.294284245, 0.942909112, 0.519258899, 0.008882952, 0.055022687,
0.563962346, 0.159175889, 0.729845947, 0.610946838, 0.063839729,
0.542889465, 0.154845574, 0.396166658, 0.512275225, 0.636128255,
0.000265, 0.016868837, 0.041575921, 0.004828325, 0.170143423,
0.057359331, 0.463224782, 0.55416784, 0.286603865, 0.415871141,
0.473120067, 0.464303708, 0.011002312, 0.5, 0.348898333, 0.220983939,
0.009418115, 0.906935504, 0.601641288, 0.480741101, 0.79015601,
0.085766864, 0.016978584, NA, NA, 0.399770285, 0.927749655, 0.029916602,
0.119047989, 0.568590421, 0.540283117, 0.654617786, 0.917625002,
0.002563335, 0.010137604, 0.043544506, 0.085651037, 0.072239359
), hypo = c("no", "no", "no", "no", "no", "no", "no", "no", "no",
"no", "no", "no", "no", "yes", "yes", "no", "no", "no", "no",
"no", "yes", "yes", "no", "yes", "no", "no", "no", "no", "no",
"no", "yes", "yes", "yes", "yes", "yes", "no", "no", "yes", "no",
"yes", "yes", "no", "no", "no", "no", "no", "yes", "yes", "yes",
"no", "yes", "yes", "yes", "yes", "yes", "no", "no", "no", "no",
"no", "yes", "no", "yes", "no", "no", "no", "yes", "yes", "yes",
"yes", "no", "yes", "yes", "no", "no", "no", "no", "no", "no",
"no", "no", "no", "no", "no", "no", "no", "yes", "no", "no",
"no", "no", "no", "no", "no", "no", "no", "no", "yes", "yes",
"no", "no", "no", "no", "no", "yes", "yes", "no", "yes", "no",
"no", "yes", "no", "no", "no", "no", "no", "no", "no", "yes",
"no", "no", "no", "no", "no", "no", "no", "no", "no", "no", "no",
"yes", "yes", "yes", "yes", "no", "no", "no", "no", "no", "no",
"no", "no", "yes", "no", "no", "no", "yes", "no", "no", "no",
"no", "no", "yes", "yes", "yes", "no", "no", "yes", "no", "no",
"no", "no", "no", "yes", "yes", "yes", "no", "no")), row.names = c(NA,
-168L), spec = structure(list(cols = list(stnn = structure(list(), class = c("collector_character",
"collector")), wscore = structure(list(), class = c("collector_double",
"collector")), lat = structure(list(), class = c("collector_double",
"collector")), lon = structure(list(), class = c("collector_double",
"collector")), pval = structure(list(), class = c("collector_double",
"collector")), MWU_res = structure(list(), class = c("collector_character",
"collector")), hydro_zone = structure(list(), class = c("collector_character",
"collector"))), default = structure(list(), class = c("collector_guess",
"collector"))), class = "col_spec"), class = c("tbl_df", "tbl",
"data.frame"))
我读到您对 "phase A" 和 "phase B" 之间的差异进行了 168 次显着性检验。一个鲜为人知的事实是，零假设下的 p 值在 [0-1] 范围内均匀分布。因此，如果阶段之间没有差异，您会期望 "significant results"，即 p 值 <= 0.05，正好是 5% 的时间（给出或采用统计不精确性）。您清楚地看到 "non-null" 结果分布。
首先，我模拟了一些可能看起来像您的原始数据的东西来执行 t.test，前 40 个是从 A 阶段和 B 阶段之间具有不同均值的正态分布模拟的，最后 128 个来自相同的分布。
set.seed(100)
dat1 = lapply(1:40,function(i){
data.frame(
stn = paste0("Stn_",i),
value = c(rnorm(15,0,1),rnorm(15,2,1)),
type = rep(c("PhaseA","PhaseB"),each=15)
)
})
dat2 = lapply(41:168,function(i){
data.frame(
stn = paste0("Stn_",i),
value = rnorm(30),
type = rep(c("PhaseA","PhaseB"),each=15)
)
})
dat = do.call(rbind,c(dat1,dat2))
现在我们可以像您一样在 phaseA 和 phaseB 之间像您一样执行 t.test，分别针对每个站：
library(dplyr)
library(broom)
library(purrr)
obs_result = dat %>% group_by(stn) %>% do(tidy(t.test(value~type,data=.)))
如果我们看到哪个站的 p < 0.05:
which(obs_result$p.value<0.05)
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
[20] 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
[39] 39 40 52 57 67 90 113 124 166
所以前 40 个（真正不同的）的大部分测试显示 p < 0.05。但是您确实也有一些来自没有显示差异的电台的点击。因此，您需要通过执行以下操作来控制多次测试：
which(p.adjust(obs_result$p.value,"BH")<0.1)
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
[20] 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
[39] 39 40 57 90 124 166
所以你有 44 次命中.. 其中 4 次实际上是错误的 positives.This 没问题，因为我们允许 FDR 为 0.1。
我们还可以模拟您在 null 下的预期。在您的情况下，我们只是交换站点的标签并再次执行 t.test。这一次，您的 p 值应该与您观察到的非常不同：
permutated_pvalue = function(i,dat){
set.seed(i)
dat %>% group_by(stn) %>%
mutate(type=sample(type)) %>%
do(tidy(t.test(value~type,data=.))) %>%
pull(p.value)
}
par(mfrow=c(1,2))
hist(permutated_pvalue(0,dat),main="permutated pvalue")
hist(obs_result$p.value,main="observed pvalue")
最后，做一个 qqplot 来查看这些 p 值分布有何不同：
# 10 simulations to get null p-values
simulated_pvalues = 1:10 %>% map(perm,dat=dat) %>% unlist()
qqplot(-log10(simulated_pvalues),-log10(obs_result$p.value))
abline(a=0,b=1)
所以我建议你做这些检查，如果图表确实显示 p 值分布不同，那么就有一些全局差异的说法..