How to perform field significance of a test in R

1
votes

I run a hypothesis test for 168 points to see if the values of the observations during A phase are significantly different than phase B and concluded that 30% of these points show a significant difference between two phases. How can I test the field significance of these results? Is it a bootstrapping method?In that case how to run bootstrap in this data. The dtt has p-value for each hypothesis testing and results.

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"))
2
rstatisticsbootstrap-modalhypothesis-test
It sounds like you are trying to determine whether there was a significant overall difference between phase A and phase B at 168 different points. The method you are using to analyse this doesn't seem right. It's difficult to say without knowing more about your original data's structure, but perhaps you should be doing a paired t-test or a mixed effects anova. – Allan Cameron
First of you need to correct for multiple testing, and if you do a Benjamini-Hochberg correction and you can see that with an FDR of 10%, table(p.adjust(dtt$pval,"BH")<0.1), about you get about 10 hits which is ok.. – StupidWolf
If you simulate data, by resampling your station label and doing the t.test, you will be able to see that first of all, you get less p < 0.05 and after FDR, you don't get so many hits.. This way you can show that you are getting significant results in a way – StupidWolf
You can do an anova for the overall effect of phase, but this will have to assume that the effect of phase B vs phase A is consistent, meaning same direction. Is this true? – StupidWolf

2 Answers

1
votes

I read that you have 168 tests of significance of the difference between "phase A" and "phase B". It's a little known fact that p-values under the null-hypothesis are distributed uniformly over the range [0-1]. So if there were no difference between phases, you would expect "significant results", i.e. p-values <= 0.05, exactly 5% of the time (give or take statistical inexactness). You are clearly seeing a "non-null" distribution of results.

1
votes

First I simulate something that might look like your original data to do the t.test, first 40 are simulated from normal distribution with different means between phase A and phase B, last 128 comes from the same distribution.

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))

Now we can do the t.test like you did, between phaseA and phaseB like what you did, separately for each station:

library(dplyr)
library(broom)
library(purrr)
obs_result = dat %>% group_by(stn) %>% do(tidy(t.test(value~type,data=.)))

If we see which station have 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

So most of the test for 1st 40 (truly different) show p < 0.05. But you do have some hits also from stations that don't show differences. Hence you need to control for multiple testing, by doing:

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

So you have 44 hits.. out of which 4 are actually false positives.This is ok as we allowed for an FDR of 0.1.

We can also simulate what you would expect under the null. In your case, we just swap the labels for the stations and perform the t.test again. This time, you p-values should look very different from what you observed:

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")

enter image description here

Lastly, do a qqplot to see how these distribution of p-values differ:

# 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)

enter image description here

So I suggest you do these checks and if indeed the plots show that the p-value distributions are different, then there's some claim of global difference..