Supplementary R code for reproducing the examples in the paper ``Comparing three groups’’.
One-way ANOVA example
This example is taken from Dobson’s (1983) book, Section 7.
Genetically similar seeds are randomly assigned to be raised either under standard conditions (control) or in two different nutritionally enriched environments (treatment A and treatment B). After a predetermined period all plants are harvested, dried and weighed. The results, expressed as dried weight in grams, for samples of \(10\) plants from each group are given in following Table.
|
Control
|
TreatmentA
|
TreatmentB
|
|
4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14
|
4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69
|
6.31, 5.12, 5.54, 5.50, 5.37, 5.29, 4.92, 6.15, 5.80, 5.26
|
We assume the linear model \[Y_{u,i} = \mu_i + \varepsilon_{u,i}, \qquad i=1,2,3, \quad u=1,\ldots,10,\] where \(Y_{u,i}\) is the response (dried weight) for unit \(u\) in group \(i\) (\(i=1\) for control, \(i=2\) for treatment A and \(i=3\) for treatment B), and \(\varepsilon_{u,i} \stackrel{\mathrm{iid}}{\sim} N(0,\sigma^2)\). We are interested in testing the global null hypothesis that all three group means are equal, i.e. \(H_{123}: \mu_1 = \mu_2=\mu_3\), and the three null hypotheses about pairwise comparisons between groups, i.e. \(H_{12}: \mu_1=\mu_2\), \(H_{13}: \mu_1=\mu_3\) and \(H_{23}: \mu_2=\mu_3\).
To get estimates \(\hat{\mu}_1\), \(\hat{\mu}_2\), \(\hat{\mu}_3\) and \(\hat{\sigma}^2\):
fit <- lm(weight ~ group -1, plant.dried)
fit$coefficients
#> groupControl groupTreatmentA groupTreatmentB
#> 5.032 4.661 5.526
summary(fit)$sigma^2
#> [1] 0.3885959
Partial \(F\)-test unadjusted \(p\)-values for \(H_{12}\), \(H_{13}\) and \(H_{23}\):
fit_12 <- lm(weight ~ I(group=="TreatmentB"), plant.dried)
p_12 <- anova(fit_12,fit)[2,"Pr(>F)"]
p_12
#> [1] 0.1943879
fit_13 <- lm(weight ~ I(group=="TreatmentA"), plant.dried)
p_13 <- anova(fit_13,fit)[2,"Pr(>F)"]
p_13
#> [1] 0.08768168
fit_23 <- lm(weight ~ I(group=="Control"), plant.dried)
p_23 <- anova(fit_23,fit)[2,"Pr(>F)"]
p_23
#> [1] 0.004459236
Tukey’s HSD adjusted \(p\)-values for \(H_{12}\), \(H_{13}\) and \(H_{23}\):
p_tuk <- TukeyHSD(aov(weight ~ group, data = plant.dried))$group[,'p adj']
p_tuk
#> TreatmentA-Control TreatmentB-Control TreatmentB-TreatmentA
#> 0.39087114 0.19799599 0.01200642
Dunnett adjusted \(p\)-values for \(H_{12}\) and \(H_{13}\)
p_dun <- DunnettTest(weight ~ factor(group), plant.dried, control="Control")$Control[,'pval']
p_dun
#> TreatmentA-Control TreatmentB-Control
#> 0.3226957 0.1534859
ANOVA F-test (method A), closed Tukey (method B), closed Dunnett (method C) and Gatekeeping (method D) \(p\)-values for \(H_{123}\):
fit_123 <- lm(weight ~ 1, plant.dried)
p_A <- anova(fit_123,fit)[2,"Pr(>F)"]
p_A
#> [1] 0.01590996
p_B <- min(p_tuk)
p_B
#> [1] 0.01200642
p_C <- min(p_dun)
p_C
#> [1] 0.1534859
p_D <- p_12
p_D
#> [1] 0.1943879
Adjusted \(p\)-values for the four hypotheses and the four methods:
adjp <- matrix(NA,nrow=4,ncol=4,
dimnames = list(c("A", "B", "C", "D"),
c("H12", "H13", "H23", "H123")))
adjp["A",] <- pmax(c(p_12,p_13,p_23,p_A),p_A)
adjp["B",] <- pmax(c(p_12,p_13,p_23,p_B),p_B)
adjp["C",] <- pmax(c(p_12,p_13,p_23,p_C),p_C)
adjp["D",] <- pmax(c(p_12,p_13,p_23,p_D),p_D)
adjp
#> H12 H13 H23 H123
#> A 0.1943879 0.08768168 0.01590996 0.01590996
#> B 0.1943879 0.08768168 0.01200642 0.01200642
#> C 0.1943879 0.15348586 0.15348586 0.15348586
#> D 0.1943879 0.19438788 0.19438788 0.19438788
We see that at the significance level of \(\alpha=5\%\), methods C and D do not reject any hypothesis, while methods A and B reject \(H_{123}\) and \(H_{13}\).
For method A, add adjusted \(p\)-values to box-plots:
bp <- ggboxplot(plant.dried, x = "group", y = "weight")
stats <- compare_means(weight ~ group, data = plant.dried, method = "t.test")
stats$p.adj <- round(adjp["A",-4],3)
bp + stat_compare_means(method = "anova", label.y = 10) +
stat_pvalue_manual(stats, label = "p = {p.adj}", y.position = c(9, 8, 7))
Alternative analyses
The methods A, B, C and D are not tied to \(F\)-tests, and naturally generalize to other tests. For example, if we believe that variances may differ between groups when the means do, we would prefer a two-sample \(t\)-test over the \(F\)-test for the pairwise hypotheses. By replacing partial F-tests by two-samples \(t\) tests, we obtain
p_12 <- stats$p[1]
p_13 <- stats$p[2]
p_23 <- stats$p[3]
adjp["A",] <- pmax(c(p_12,p_13,p_23,p_A),p_A)
adjp["B",] <- pmax(c(p_12,p_13,p_23,p_B),p_B)
adjp["C",] <- pmax(c(p_12,p_13,p_23,p_C),p_C)
adjp["D",] <- pmax(c(p_12,p_13,p_23,p_D),p_D)
adjp
#> H12 H13 H23 H123
#> A 0.2503825 0.04789926 0.01590996 0.01590996
#> B 0.2503825 0.04789926 0.01200642 0.01200642
#> C 0.2503825 0.15348586 0.15348586 0.15348586
#> D 0.2503825 0.19438788 0.19438788 0.19438788
resulting in the rejection of \(H_{123}\), \(H_{13}\) and \(H_{23}\) at \(\alpha=5\%\) for methods A and B.Here method A is Fisher’s LSD.
We could also use methods A, B, C and D in a permutation framework. A standard choice is to consider non-standardized test statistics \(T_{ij} = (\hat{\mu}_i - \hat{\mu}_j)^2\) for testing \(H_{ij}\), and \(T_{123} = T_{12}+T_{23}+ T_{13}\), \(\tilde{T} = \max(T_{12}, T_{13},T_{23})\), \(\tilde{T}_{1}=\max(T_{12}, T_{13})\) and \(T_{12}\) for testing \(H_{123}\) in methods A, B, C and D, respectively.
The construction of the permutation null distribution for each test statistic proceeds as follows. The observations of the groups are pooled, and the test statistic is recalculated for every permutation of the group labels.
Permutation tests of \(H_{123}\) uses a global permutation distribution, constructed by permuting the observations of all three groups.
labels <- factor(plant.dried$group, labels=0:2)
labels <- as.numeric(levels(labels))[labels]
B = 10^5
labels.perm <- mt.sample.label(labels, test="f", B=B)
c12 <- ( (labels.perm + 2) %% 3 ) -1
c13 <- labels.perm - 1
c23 <- ( (labels.perm + 1) %% 3 ) -1
y <- plant.dried$weight
n <- 10
stats.perm <- rbind(
( y %*% t(c12) / n )^2,
( y %*% t(c13) / n )^2,
( y %*% t(c23) / n )^2
)
hist(colSums(stats.perm), main="", xlab=expression(T[123]), xatx="n")
points(sum(stats.perm[,1]),0, pch=19)
hist(apply(stats.perm,2,max), main="", xlab=expression( tilde(T) ) )
points(max(stats.perm[,1]),0, pch=19)
The permutation \(p\)-value is calculated as the proportion of permutations where the test statistic is greater than or equal to the value computed on the original data.
p_A <- mean( colSums(stats.perm) >= sum(stats.perm[,1]) )
p_A
#> [1] 0.01668
p_B <- mean( apply(stats.perm,2,max) >= max(stats.perm[,1]) )
p_B
#> [1] 0.01136
p_C <- mean( apply(stats.perm[-3,],2,max) >= max(stats.perm[-3,1]) )
p_C
#> [1] 0.20692
The permutation version of the ANOVA F-test can also be obtained by
t_123 <- mt.sample.teststat(y, labels, B=B, test="f")
p_A <- mean( t_123 >= t_123[1] )
p_A
#> [1] 0.01668
Permutation tests of \(H_{ij}\) use a local permutation distribution, constructed by permuting the observations of groups \(i\) and \(j\).
t_12 <- (mt.sample.teststat(y[labels!=2],labels[labels!=2], B=0, test="t.equalvar") ) ^2
#>
#> We're doing 184756 complete permutations
p_12 <- mean(t_12 >= t_12[1])
p_12
#> [1] 0.2471801
t_13 <- (mt.sample.teststat(y[labels!=1],labels[labels!=1]=="2", B=0, test="t.equalvar") ) ^2
#>
#> We're doing 184756 complete permutations
p_13 <- mean(t_13 >= t_13[1])
p_13
#> [1] 0.04797679
t_23 <- (mt.sample.teststat(y[labels!=0],labels[labels!=0]=="2", B=0, test="t.equalvar") ) ^2
#>
#> We're doing 184756 complete permutations
p_23 <- mean(t_23 >= t_23[1])
p_23
#> [1] 0.00853017
Adjusted \(p\)-values with permutation tests for the four hypotheses and the four methods:
p_D <- p_12
adjp["A",] <- pmax(c(p_12,p_13,p_23,p_A),p_A)
adjp["B",] <- pmax(c(p_12,p_13,p_23,p_B),p_B)
adjp["C",] <- pmax(c(p_12,p_13,p_23,p_C),p_C)
adjp["D",] <- pmax(c(p_12,p_13,p_23,p_D),p_D)
adjp
#> H12 H13 H23 H123
#> A 0.2471801 0.04797679 0.0166800 0.0166800
#> B 0.2471801 0.04797679 0.0113600 0.0113600
#> C 0.2471801 0.20692000 0.2069200 0.2069200
#> D 0.2471801 0.24718006 0.2471801 0.2471801
Rank tests can be obtained as a special case of permutation tests by replacing the observations with their ranks.
r <- rank(y, ties.method = "average")
ranks.perm <- rbind(
( r %*% t(c12) / n )^2,
( r %*% t(c13) / n )^2,
( r %*% t(c23) / n )^2
)
p_A <- mean( colSums(ranks.perm) >= sum(ranks.perm[,1]) )
p_B <- mean( apply(ranks.perm,2,max) >= max(ranks.perm[,1]) )
p_C <- mean( apply(ranks.perm[-3,],2,max) >= max(ranks.perm[-3,1]) )
Kruskal-Wallis test can also be obtained by
w_123 <- mt.sample.teststat(y, labels, B=B, test="f", nonpara = "y")
p_A <- mean( w_123 >= w_123[1] )
p_A
#> [1] 0.01442
Wilcoxon-Mann-Whitney tests are obtained by
w_12 <- (mt.sample.teststat(y[labels!=2],labels[labels!=2], B=0, test="wilcoxon"))^2
#>
#> We're doing 184756 complete permutations
p_12 <- mean(w_12 >= w_12[1])
w_13 <- (mt.sample.teststat(y[labels!=1],labels[labels!=1]=="2", B=0, test="wilcoxon"))^2
#>
#> We're doing 184756 complete permutations
p_13 <- mean(w_13 >= w_13[1])
w_23 <- (mt.sample.teststat(y[labels!=0],labels[labels!=0]=="2", B=0, test="wilcoxon"))^2
#>
#> We're doing 184756 complete permutations
p_23 <- mean(w_23 >= w_23[1])
p_D <- p_12
Adjusted \(p\)-values with rank tests for the four hypotheses and the four methods:
adjp["A",] <- pmax(c(p_12,p_13,p_23,p_A),p_A)
adjp["B",] <- pmax(c(p_12,p_13,p_23,p_B),p_B)
adjp["C",] <- pmax(c(p_12,p_13,p_23,p_C),p_C)
adjp["D",] <- pmax(c(p_12,p_13,p_23,p_D),p_D)
adjp
#> H12 H13 H23 H123
#> A 0.1967568 0.06301284 0.0144200 0.0144200
#> B 0.1967568 0.06301284 0.0106300 0.0106300
#> C 0.1967568 0.16581000 0.1658100 0.1658100
#> D 0.1967568 0.19675680 0.1967568 0.1967568
ANCOVA example
As with ANOVA, we are interested in comparing means for groups defined by a factor while controlling for the effects of other covariates that are not of primary interest. The following table displays data from Winer (1971), discussed in Dobson (1983).
|
Group A
|
Group B
|
Group C
|
|
y
|
x
|
y
|
x
|
y
|
x
|
|
6
|
3
|
8
|
4
|
6
|
3
|
|
4
|
1
|
9
|
5
|
7
|
2
|
|
5
|
3
|
7
|
5
|
7
|
2
|
|
3
|
1
|
9
|
4
|
7
|
3
|
|
4
|
2
|
8
|
3
|
8
|
4
|
|
3
|
1
|
5
|
1
|
5
|
1
|
|
6
|
4
|
7
|
2
|
7
|
4
|
The response \(y\) is the achievement score, the levels A, B and C of the group factor represent three different training methods, and the covariate \(x\) is the aptitude score measured before training commenced. We want to compare the training methods, taking into account differences in initial aptitude between the three groups of subjects.
achievement$method <- factor(achievement$method)
ggscatter(achievement, x="x", y="y", color = "method")
We assume that the response in group \(i\) is normally distributed with mean \(\mu_i(x)\) and variance \(\sigma^2\), with \[\mu_i(x) = \gamma + \tau_i + \beta (x - \bar{x})\] where \(\gamma\) is the common mean, \(\tau_i\) is the \(i\)th group effect such that \(\sum_i \tau_i =0\), \(\beta\) is the regression slope and \(\bar{x}\) is the average covariate value, with \(i=1,2,3\) for group A,B,C, respectively.
Analysis of covariance compares the adjusted means \(\hat{\mu}_i(\bar{x}_i)\), i.e. the estimated group means adjusted for the group average covariate values.
Estimates \(\hat\mu_1(\bar{x}_1)\), \(\hat\mu_2(\bar{x}_2)\) and \(\hat\mu_3(\bar{x}_3)\):
fit <- lm(y ~ x + method, data = achievement)
effect("method", fit)
method effect
method
A B C
4.888435 7.076190 6.749660
Graph of the fitted model:
plot(predictorEffects(fit))
Let \(Y\) and \(X\) denote the response vector and the design matrix, respectively, and let \(\hat{\theta} = (X'X)^{-1}XY\) be the least square estimator.
Least square estimates \(\hat{\theta}\) and \(\hat{\sigma}^2\):
X <- model.matrix(fit)
df <- summary(fit)$df[2]
theta <- coef(fit)
theta
#> (Intercept) x methodB methodC
#> 2.8367347 0.7428571 2.1877551 1.8612245
s2 <- summary(fit)$sigma^2
s2
#> [1] 0.6060024
We can use the multcomp package to calculate the test statistics and find their distribution.
Matrix \(K\) specifying pairwise contrasts of factor levels:
tuk <- glht(fit,linfct=mcp(method="Tukey"))
K <- tuk$linfct
K
(Intercept) x methodB methodC
B - A 0 0 1 0
C - A 0 0 0 1
C - B 0 0 -1 1
attr(,"type")
[1] "Tukey"
Estimated differences \(K\hat{\theta} = \left(\begin{array}{c} \hat{\mu}_2(\bar{x}_2) - \hat{\mu}_1(\bar{x}_1)\\ \hat{\mu}_3(\bar{x}_3) - \hat{\mu}_1(\bar{x}_1)\\ \hat{\mu}_3(\bar{x}_3) - \hat{\mu}_2(\bar{x}_2)\\ \end{array}\right)\)
K %*% theta
#> [,1]
#> B - A 2.1877551
#> C - A 1.8612245
#> C - B -0.3265306
and estimated covariance matrix \(\hat{\Sigma} = \hat{\sigma}^2 K(X'X)^{-1}K'\)
S <- s2 * K %*% solve(crossprod(X)) %*% t(K)
S
#> B - A C - A C - B
#> B - A 0.2065355 0.10141265 -0.10512287
#> C - A 0.1014126 0.17973949 0.07832684
#> C - B -0.1051229 0.07832684 0.18344971
Unadjusted \(p\)-values \(p_{12}\), \(p_{13}\) and \(p_{23}\)
p_raw <- pf( (K %*% theta)^2 / diag(S), df1 = 1, df2 = df, lower.tail = F)
p_12 <- p_raw[1]
p_12
#> [1] 0.0001620169
p_13 <- p_raw[2]
p_13
#> [1] 0.0003995625
p_23 <- p_raw[3]
p_23
#> [1] 0.456288
Tukey’s HSD adjusted \(p\)-values \(\tilde{p}^{Tuk}_{12}\), \(\tilde{p}^{Tuk}_{13}\) and \(\tilde{p}^{Tuk}_{23}\) :
p_tuk <- summary(tuk)$test$pvalues[1:3]
p_tuk
#> [1] 0.0004274782 0.0011832849 0.7302152807
Dunnett adjusted \(p\)-values \(\tilde{p}^{Dun}_{12}\) and \(\tilde{p}^{Dun}_{13}\)
dun <- glht(fit,linfct=mcp(method="Dunnett"))
p_dun <- summary(dun)$test$pvalues[1:2]
p_dun
#> [1] 0.0003111268 0.0007633035
ANCOVA F-test (method A), closed Tukey (method B), closed Dunnett (method C) and Gatekeeping (method D) \(p\)-values for \(H_{123}\):
fit_123 <- lm(y ~ x, achievement)
p_A <- anova(fit_123,fit)[2,"Pr(>F)"]
p_B <- min(p_tuk)
p_C <- min(p_dun)
p_D <- p_12
Adjusted \(p\)-values for the four hypotheses and the four methods:
adjp <- matrix(NA,nrow=4,ncol=4,
dimnames = list(c("A", "B", "C", "D"),
c("H12", "H13", "H23", "H123")))
adjp["A",] <- pmax(c(p_12,p_13,p_23,p_A),p_A)
adjp["B",] <- pmax(c(p_12,p_13,p_23,p_B),p_B)
adjp["C",] <- pmax(c(p_12,p_13,p_23,p_C),p_C)
adjp["D",] <- pmax(c(p_12,p_13,p_23,p_D),p_D)
adjp
#> H12 H13 H23 H123
#> A 0.0002578664 0.0003995625 0.456288 0.0002578664
#> B 0.0004274782 0.0004274782 0.456288 0.0004274782
#> C 0.0003111268 0.0003995625 0.456288 0.0003111268
#> D 0.0001620169 0.0003995625 0.456288 0.0001620169
Four-group example
Jakob et al. (2019) collected a dataset from 847 subjects to study the relation between Harry Potter Houses and personality traits. The raw data can be found at https://osf.io/rtf74/.
In this analysis, the outcome variable is the Machiavellianism score, as measured by the Dark Triad Questionnaire (Jones and Paulhus, 2014), and the affiliation to four Houses of Hogwarts - Gryffindor, Hufflepuff, Ravenclaw, or Slytherin - is determined by the Sorting Hat Quiz, available on https://www.pottermore.com.
rm(list=ls())
set.seed(123)
source_data("https://github.com/aldosolari/comparing3groups/blob/main/Harry4.RData?raw=true")
#> Downloading data from: https://github.com/aldosolari/comparing3groups/blob/main/Harry4.RData?raw=true
#> SHA-1 hash of the downloaded data file is:
#> 748c8a18d72cdbc1b94960145e0257bc23e56b2f
#> [1] "Harry4"
table(Harry4$group)
#>
#> Gryffindor Hufflepuff Ravenclaw Slytherin
#> 43 36 42 24
ggboxplot(Harry4, x = "group", y = "score", col="group", add = "jitter")
To get estimates \(\hat{\mu}_1\), \(\hat{\mu}_2\), \(\hat{\mu}_3\) and \(\hat{\mu}_4\):
levels(Harry4$group) <- 1:4
fit <- lm(score ~ group -1, Harry4)
fit$coefficients
#> group1 group2 group3 group4
#> 27.48837 26.91667 27.28571 32.04167
To compute unadjusted \(p\)-values for \(H_{12}\), \(H_{13}\), \(H_{14}\), \(H_{23}\), \(H_{24}\) and \(H_{34}\):
tuk <- glht(fit,linfct=mcp(group="Tukey"))
p_pair = summary(tuk, test=adjusted("none"))$test$pvalues
round(p_pair,4)
#> 2 - 1 3 - 1 4 - 1 3 - 2 4 - 2 4 - 3
#> 0.7015 0.8875 0.0075 0.8056 0.0037 0.0055
Define the 14 hypotheses and the inclusion relationship between each elementary hypothesis and its ancestors:
hyp = list(`1234` = 1:4,
`123` = 1:3,
`124` = c(1:2,4),
`134` = c(1,3:4),
`234` = 2:4,
`12_34` = list(1:2,3:4),
`13_24` = list(c(1,3),c(2,4)),
`14_23` = list(c(1,4),2:3),
`12` = 1:2,
`13` = c(1,3),
`14` = c(1,4),
`23` = 2:3,
`24` = c(2,4),
`34` = 3:4)
subhyp <- sapply(hyp[c("12","13","14","23","24","34")],
function(y)
sapply(hyp, function(x) ifelse(any(lengths(x)<2),
all(y %in% x),
list(y) %in% x)))
To compute adjusted \(p\)-values for \(H_{12}\), \(H_{13}\), \(H_{14}\), \(H_{23}\), \(H_{24}\) and \(H_{34}\) for classic closed testing with partial \(F\)-tests:
rawp_A <- vector("numeric", length = length(hyp))
names(rawp_A) <- names(hyp)
rawp_A["1234"] <- anova(lm(score ~ 1, Harry4),fit)[2,"Pr(>F)"]
rawp_A[c("123","124","134","234")] <- sapply(4:1, function(i)
anova(lm(score ~ I(group==i), Harry4),fit)[2,"Pr(>F)"] )
anova_12_34 <- glht(fit, linfct=c("group1 - group2 = 0",
"group3 - group4 = 0"))
rawp_A["12_34"] <-summary(anova_12_34, test=Ftest())$test$pvalue
anova_13_24 <- glht(fit, linfct=c("group1 - group3 = 0",
"group2 - group4 = 0"))
rawp_A["13_24"] <-summary(anova_13_24, test=Ftest())$test$pvalue
anova_14_23 <- glht(fit, linfct=c("group1 - group4 = 0",
"group2 - group3 = 0"))
rawp_A["14_23"] <-summary(anova_14_23, test=Ftest())$test$pvalue
rawp_A[c("12","13","14","23","24","34")] <- p_pair
adjp_A <- apply(subhyp, 2, function(x) max(rawp_A[x]) )
Closed Tukey adjusted \(p\)-values:
rawp_B <- vector("numeric", length = length(hyp))
names(rawp_B) <- names(hyp)
rawp_B["1234"] <- min(summary(tuk)$test$pvalues)
tuk_123 <- glht(fit, linfct=c("group1 - group2 = 0",
"group1 - group3 = 0",
"group2 - group3 = 0"))
rawp_B["123"] <- min(summary(tuk_123)$test$pvalues)
tuk_124 <- glht(fit, linfct=c("group1 - group2 = 0",
"group1 - group4 = 0",
"group2 - group4 = 0"))
rawp_B["124"] <- min(summary(tuk_124)$test$pvalues)
tuk_134 <- glht(fit, linfct=c("group1 - group3 = 0",
"group1 - group4 = 0",
"group3 - group4 = 0"))
rawp_B["134"] <- min(summary(tuk_134)$test$pvalues)
tuk_234 <- glht(fit, linfct=c("group2 - group3 = 0",
"group2 - group4 = 0",
"group3 - group4 = 0"))
rawp_B["234"] <- min(summary(tuk_234)$test$pvalues)
tuk_12_34 <- glht(fit, linfct=c("group1 - group2 = 0",
"group3 - group4 = 0"))
rawp_B["12_34"] <- min(summary(tuk_12_34)$test$pvalues)
tuk_13_24 <- glht(fit, linfct=c("group1 - group3 = 0",
"group2 - group4 = 0"))
rawp_B["13_24"] <- min(summary(tuk_13_24)$test$pvalues)
tuk_14_23 <- glht(fit, linfct=c("group1 - group4 = 0",
"group2 - group3 = 0"))
rawp_B["14_23"] <- min(summary(tuk_14_23)$test$pvalues)
rawp_B[c("12","13","14","23","24","34")] <- p_pair
adjp_B <- apply(subhyp, 2, function(x) max( rawp_B[x]) )
Closed Dunnett adjusted \(p\)-values:
rawp_C <- vector("numeric", length = length(hyp))
names(rawp_C) <- names(hyp)
dun <- glht(fit, linfct=c("group4 - group1 = 0",
"group4 - group2 = 0",
"group4 - group3 = 0"))
rawp_C["1234"] <- min(summary(dun)$test$pvalues)
rawp_C["123"] <- rawp_A["123"]
dun_124 <- glht(fit, linfct=c("group4 - group1 = 0",
"group4 - group2 = 0"))
rawp_C["124"] <- min(summary(dun_124)$test$pvalues)
dun_134 <- glht(fit, linfct=c("group4 - group1 = 0",
"group4 - group3 = 0"))
rawp_C["134"] <- min(summary(dun_134)$test$pvalues)
dun_234 <- glht(fit, linfct=c("group4 - group3 = 0",
"group4 - group2 = 0"))
rawp_C["234"] <- min(summary(dun_234)$test$pvalues)
rawp_C[c("14_23","13_24","12_34")] <- rawp_A[c("14","24","34")]
rawp_C[c("12","13","14","23","24","34")] <- p_pair
adjp_C <- apply(subhyp, 2, function(x) max( rawp_C[x]) )
Gatekeeping adjusted \(p\)-values:
rawp_D <- vector("numeric", length = length(hyp))
names(rawp_D) <- names(hyp)
rawp_D[c("1234",
"123","124","134","234",
"12_34","13_24","14_23")] <- rawp_A[c("24",
"12","24","14","24",
"34","24","14")]
rawp_D[c("12","13","14","23","24","34")] <- p_pair
adjp_D <- apply(subhyp, 2, function(x) max(rawp_D[x]) )
Adjusted \(p\)-values for the six hypotheses and four methods
adjp <- matrix(rbind(adjp_A, adjp_B, adjp_C, adjp_D),nrow=4,ncol=6,
dimnames = list(c("A", "B", "C","D"),names(p_pair)))
round(adjp,4)
#> 2 - 1 3 - 1 4 - 1 3 - 2 4 - 2 4 - 3
#> A 0.9279 0.9279 0.0269 0.9279 0.0155 0.0194
#> B 0.9219 0.9219 0.0190 0.9219 0.0190 0.0190
#> C 0.9279 0.9279 0.0101 0.9279 0.0101 0.0101
#> D 0.7015 0.8875 0.0075 0.8056 0.0037 0.0075
Compare to Tukey’s HSD adjusted \(p\)-values:
adjp_tuk = c(summary(tuk)$test$pvalues)
names(adjp_tuk) <- rownames(summary(tuk)$linfct)
round(adjp_tuk,4)
#> 2 - 1 3 - 1 4 - 1 3 - 2 4 - 2 4 - 3
#> 0.9805 0.9990 0.0369 0.9947 0.0194 0.0276
Compare to Dunnett’s adjusted \(p\)-values:
adjp_dun = c(summary(dun)$test$pvalues)
names(adjp_dun) <- rownames(summary(dun)$linfct)
round(adjp_dun,4)
#> group4 - group1 group4 - group2 group4 - group3
#> 0.0196 0.0097 0.0142