kornl / gMCP

gMCP - Graph Based Multiple Comparison Procedures (an R package with a Java GUI)
http://gsrmtp.r-forge.r-project.org/
10 stars 1 forks source link

[Question] Could you, please, cross-validate the calculation gMCP vs. Mediana in this scenario? #18

Closed Generalized closed 1 year ago

Generalized commented 1 year ago

For this 2-family parallel gatekeeper with Hochberg truncated at: 0.5 in the 1st family and 1 (classic) in the 2nd family:

obraz

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)

obraz

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

Let's filter it: obraz

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 1 year ago

OK, I now understand. This is a different procedure. The weighted Simes, although more powerful than Hochberg, provides only weak control of the FWER.