[Question] Could you, please, cross-validate the calculation gMCP vs. Mediana in this scenario?
Generalized opened this issue · 1 comments
Generalized commented
For this 2-family parallel gatekeeper with Hochberg truncated at: 0.5 in the 1st family and 1 (classic) in the 2nd family:
I get a discrepancy between gMCP and Mediana.
Mediana:
rawp <- c(0.01, 0.013, 0.01, 0.01)
family <- families(family1 = c(1, 2), family2 = c(3, 4))
gamma <- families(family1 = 0.5, family2 = 1)
AdjustPvalues(rawp,
proc = "ParallelGatekeepingAdj",
par = parameters(family = family,
proc = families(family1 ="HochbergAdj", family2 = "HochbergAdj"),
gamma = gamma))
[1] 0.01733333 0.01733333 0.01733333 0.01733333
gMCP with Simes test:
m <- rbind(H1=c(0, 0.5, 0.25, 0.25),
H2=c(0.5, 0, 0.25, 0.25),
H3=c(0, 0, 0, 1),
H4=c(0, 0, 1, 0))
weights <- c(0.5, 0.5, 0, 0)
graph <- new("graphMCP", m=m, weights=weights)
pvalues <- c(0.01, 0.013, 0.01, 0.01)
gMCP(graph, pvalues, test="Simes", alpha=0.05)
I get from the CTP:
Remaining hypotheses (new numbering):
1: H1
2: H2
3: H3
4: H4
Subset {4} (padj: 0.01): p_4=0.01<=a*(w_4) =0.05*(1)=0.05
Subset {3} (padj: 0.01): p_3=0.01<=a*(w_3) =0.05*(1)=0.05
Subset {3,4} (padj: 0.01): p_4=0.01<=a*(w_3+w_4) =0.05*(0.5+0.5)=0.05
Subset {2} (padj: 0.0173333333333333): p_2=0.013<=a*(w_2) =0.05*(0.75)=0.0375
Subset {2,4} (padj: 0.013): p_4=0.01<=a*(w_4) =0.05*(0.25)=0.0125
Subset {2,3} (padj: 0.013): p_3=0.01<=a*(w_3) =0.05*(0.25)=0.0125
Subset {2,3,4} (padj: 0.013): p_4=0.01<=a*(w_3+w_4) =0.05*(0.125+0.125)=0.0125
Subset {1} (padj: 0.0133333333333333): p_1=0.01<=a*(w_1) =0.05*(0.75)=0.0375
Subset {1,4} (padj: 0.01): p_4=0.01<=a*(w_1+w_4) =0.05*(0.75+0.25)=0.05
Subset {1,3} (padj: 0.01): p_3=0.01<=a*(w_1+w_3) =0.05*(0.75+0.25)=0.05
Subset {1,3,4} (padj: 0.01): p_4=0.01<=a*(w_1+w_3+w_4) =0.05*(0.75+0.125+0.125)=0.05
Subset {1,2} (padj: 0.013): p_2=0.013<=a*(w_1+w_2) =0.05*(0.5+0.5)=0.05
Subset {1,2,4} (padj: 0.013): p_4=0.01<=a*(w_1+w_4) =0.05*(0.5+0)=0.025
Subset {1,2,3} (padj: 0.013): p_3=0.01<=a*(w_1+w_3) =0.05*(0.5+0)=0.025
Subset {1,2,3,4} (padj: 0.013): p_4=0.01<=a*(w_1+w_3+w_4) =0.05*(0.5+0+0)=0.025
For the same parameters and Holm (Bonferroni's approach) it agrees.
EDIT: Checked also with multxpert:
F1 <- list(label = "F1", rawp=c(0.01, 0.013), proc = "Hochberg", procpar = 0.5)
F2 <- list(label = "F2", rawp=c(0.01, 0.01), proc = "Hochberg", procpar = 1)
multxpert::pargateadjp(gateproc = list(F1, F2), independence = TRUE, printDecisionRules=TRUE)
Hypothesis testing problem
Global familywise error rate=0.05
Independence condition is imposed (the families are tested from first to last)
Family 1 (F1) is tested using Hochberg procedure (truncation parameter=0.5) at alpha1=0.05.
Null hypothesis 1 (raw p-value=0.01) is rejected.
Null hypothesis 2 (raw p-value=0.013) is rejected.
Details on the decision rule for this family can be obtained by running the PValAdjP function for Hochberg procedure with gamma=0.5 and alpha=0.05.
One or more null hypotheses are rejected in Family 1 and the parallel gatekeeping procedure passes this family. Based on the error rate function of Hochberg procedure (truncation parameter=0.5), alpha2=0.05 is carried over to Family 2.
Family 2 (F2) is tested using Hochberg procedure (truncation parameter=1) at alpha2=0.05.
Null hypothesis 3 (raw p-value=0.01) is rejected.
Null hypothesis 4 (raw p-value=0.01) is rejected.
Details on the decision rule for this family can be obtained by running the PValAdjP function for Hochberg procedure with gamma=1 and alpha=0.05.
Family Procedure Parameter Raw.pvalue Adj.pvalue
1 F1 Hochberg 0.5 0.010 0.0173
2 F1 Hochberg 0.5 0.013 0.0173
3 F2 Hochberg 1.0 0.010 0.0173
4 F2 Hochberg 1.0 0.010 0.0173
Generalized commented
OK, I now understand. This is a different procedure. The weighted Simes, although more powerful than Hochberg, provides only weak control of the FWER.