Simultaneous Directional Inference
2023-08-18
Supplementary R code for reproducing the examples in the paper
``Simultaneous Directional Inference’’.
Set up
Install the latest version of the nplus R
package:
devtools::install_github("aldosolari/nplus")
Introduction
We consider the problem of inference on the signs of \(n>1\) parameters. Let \(\theta = (\theta_1,\ldots,\theta_n)\) be a
vector of \(n\) unknown real-valued
parameters. For any \(I\subseteq
\{1,\ldots,n\}\), let \(n^+(I)\)
and \(n^-(I)\) be the number of
parameters \(\theta_i\), \(i \in I\) with positive values and negative
values, respectively: \[n^+(I) = |\{ i \in I
: \theta_i > 0 \}|, \quad n^-(I) = |\{ i \in I : \theta_i < 0 \}|.
\] For simplicity of notation, we use \(n^+\) and \(n^-\) instead of \(n^+(\{1,\ldots,n\})\) and \(n^-(\{1,\ldots,n\})\).
Our first goal is to provide post-hoc lower bounds \(\ell_{\alpha}^{+}(I)\) and \(\ell_{\alpha}^{-}(I)\) for \(n^+(I)\) and \(n^-(I)\) such that \[
\mathrm{pr}_{\theta}\Big(
n^{+}(I)
\geq
\ell_{\alpha}^{+}(I),\,\,
n^{-}(I)
\geq
\ell_{\alpha}^{-}(I),
\mathrm{\,\,for\,\,all\,\,}I \Big) \geq 1- \alpha.
\] It is worth noting that once we have obtained the confidence
lower bounds on the number of positive and negative parameters, we can
automatically derive the corresponding upper bounds, i.e. if \(n^{+}(I) \geq \ell_\alpha^+(I), n^{-}(I) \geq
\ell_\alpha^-(I)\) holds, then also \(n^{+}(I) \leq |I| - \ell^-_\alpha(I),
n^{-}_\alpha(I) \leq |I| - \ell^{+}_\alpha(I)\) holds.
For finding \(\ell_{\alpha}^{+}(I)\)
and \(\ell_{\alpha}^{-}(I)\), we use
the directional closed testing (DCT) procedure: first, select
from each pair \[
H_i^-: \theta_i \leq 0, \qquad
H_i^+: \theta_i \geq 0.
\]
the hypothesis to test. Let \(p_i\)
and \(q_i=1-p_i\) be the \(p\)-values for \(H_i^-\) and \(H_i^+\), respectively. We select \(H_i^-\) if \(p_i
\leq 1/2\) and \(H_i^+\) if
\(p_i > 1/2\).
Second, on the selected \(n\)
one-sided hypotheses, apply the closed testing procedure at level \(\alpha\). Of course, the second step has to
take care to adjust for the first step of selection from the same
data.
In some settings it is enough to infer on positive and non-positive
findings, i.e., the interest is only in lower bounds on \(n^+(I)\) and on \(n^-(I)+n^0( I)\), where \(n^0( I) = | I|-n^+(I)-n^-(I) = |i\in I: \theta_i =
0|\). Thus we would like to provide lower bounds \(\bar{\ell}_{\alpha}^{+}(I)\) and \(\bar{\ell}_{\alpha}^{-}(I)\) for \(n^+(I)\) and \(n^-(I)+n^0(I)\) such that \[
\mathrm{pr}_{\theta}\Big(
n^{+}(I)
\geq
\bar{\ell}_{\alpha}^{+}(I),\,\,
n^{-}(I)+n^{0}(I)
\geq
\bar{\ell}_{\alpha}^{-}(I),
\mathrm{\,\,for\,\,all\,\,}I \Big) \geq 1- \alpha.\] So the
hypotheses considered for \(i=1,\ldots,n\) are: \[
H_i^-: \theta_i \leq 0, \quad K_i: \theta_i > 0, \quad i=1,\ldots,n.
\] Exactly one of the hypotheses pairs is true. For finding \(\bar{\ell}_{\alpha}^+(I)\) and \(\bar{\ell}_{\alpha}^-(I)\) we use the
partitioning principle to the \(2^n\) orthants defined by \(\theta_i \in \{(-\infty,0],(0,\infty)\},
i=1,\ldots, n\). Heller & Solari (2023) showed that the
bounds obtained through the adaptive local tests with
partitioning principle (AP) uniformly improves upon the DCT bounds:
\(\bar{\ell}_{\alpha}^{+}(I)\geq
\ell_{\alpha}^{+}(I)\), \(\bar{\ell}_{\alpha}^{-}(I)\geq
\ell_{\alpha}^{-}(I)\).
Heller, R. and Solari, A. (2022). Simultaneous directional inference.
arXiv preprint arXiv:2301.01653
Subgroup Analysis
The data analysed in Gail and Simon (1985) consists of the difference
in disease-free survival probabilities at 3 years for breast cancer
patients who underwent a PFT treatment compared to those who received a
PF treatment. The study includes a total of 1260 patients divided in
\(n=4\) subgroups, classified based on
age and progesterone receptor levels. Table 2 in Gail and Simon (1985)
reports the Kaplan-Meier estimate \(\hat{\pi}_{ij}\) of the disease-free
survival probability at 3 years for patients who underwent treatment
\(j\) and belong to subgroup \(i\), along with the corresponding standard
error \(\mathrm{SE}(\hat{\pi}_{ij})\),
calculated using Greenwood’s formula.
We consider the arcsine-square root transformation \(g(\hat{\pi}_{ij}) = \arcsin{ \sqrt{ \hat{\pi}_{ij}
} }\) and compute the standardized difference \(\hat{\theta}_i= \{ g(\hat{\pi}_{i1}) -
g(\hat{\pi}_{i2}) \} /\sqrt{ 1/(4 \tilde{m}_{i1}) + 1/(4 \tilde{m}_{i2})
}\), where \(\tilde{m}_{ij}=\hat{\pi}_{ij}
(1-\hat{\pi}_{ij})\) is the effective sample size. The \(p\)-value \(p_i\) for \(H_i^-\) is calculated as \(1-\Phi(\hat{\theta}_i)\), where \(\Phi(\cdot)\) is the \(N(0,1)\) CDF.
The Kaplan-Meier estimate of survival probability at 3 years (and
standard error SE) for patients in subgroup \(i\) receiving treatment PFT or PF.
|
Subgroup i
|
PFT
|
SE(PFT)
|
PF
|
SE(PF)
|
|
1
|
0.599
|
0.0542
|
0.436
|
0.0572
|
|
2
|
0.526
|
0.0510
|
0.639
|
0.0463
|
|
3
|
0.651
|
0.0431
|
0.698
|
0.0438
|
|
4
|
0.639
|
0.0386
|
0.790
|
0.0387
|
pi_1 = c(.599, .526, .651, .639)
pi_2 = c(.436, .639, .698, .790)
SE_1 = c(.0542, .0510, .0431, .0386 )
SE_2 = c(.0572, .0463, .0438, .0387)
m1_tilde = (pi_1*(1-pi_1))/(SE_1^2)
m2_tilde = (pi_2*(1-pi_2))/(SE_2^2)
thetaHat = (asin(sqrt(pi_1)) - asin(sqrt(pi_2)))/sqrt( 1/(4*m1_tilde) + 1/(4*m2_tilde) )
round(thetaHat,3)
#> [1] 2.051 -1.635 -0.764 -2.708
p_right_tail = 1-pnorm(thetaHat)
round(p_right_tail,4)
#> [1] 0.0202 0.9490 0.7774 0.9966
Based on the \(p\)-values \(p_i\), we select the hypotheses \(H_1^-\), \(H_2^+\), \(H_3^+\), \(H_4^+\) and we adjust for the selection the
corresponding \(p\)-values:
n = length(p_right_tail)
S = rep("+",n)
S[p_right_tail <= 0.5] <- "-"
hyp_selected <- paste0("H",paste0(1:n,S))
p_selected <- 2*pmin(p_right_tail, 1-p_right_tail)
names(p_selected) <- hyp_selected
round(p_selected,4)
#> H1- H2+ H3+ H4+
#> 0.0403 0.1020 0.4451 0.0068
Function to compute the power set of a set:
powset <- function(set) {
n <- length(set)
masks <- 2^(1:n-1)
lapply( 1:2^n-1, function(u) set[ bitwAnd(u, masks) != 0 ] )
}
Function to compute the \(p\)-value
of the Simes local test:
p_combi_test <- function(x) { min(sort(x)*length(x)/(1:length(x))) }
Computing the \(p\)-value for each
intersection hypothesis by using the Simes local test:
res = unlist(
lapply(powset(1:n), function(x) p_combi_test(p_selected[x]))
)
names(res) <- powset(hyp_selected)
round(res,4)
character(0) H1-
Inf 0.0403
H2+ c("H1-", "H2+")
0.1020 0.0806
H3+ c("H1-", "H3+")
0.4451 0.0806
c("H2+", "H3+") c("H1-", "H2+", "H3+")
0.2039 0.1209
H4+ c("H1-", "H4+")
0.0068 0.0135
c("H2+", "H4+") c("H1-", "H2+", "H4+")
0.0135 0.0203
c("H3+", "H4+") c("H1-", "H3+", "H4+")
0.0135 0.0203
c("H2+", "H3+", "H4+") c("H1-", "H2+", "H3+", "H4+")
0.0203 0.0271
Adjusted \(p\)-values by closed
testing:
res_ct = res
closure = powset(1:n)
for (k in 1:length(closure))
res_ct[[k]] = max(res[sapply(1:length(closure), function(i) all(closure[[k]] %in% closure[[i]]))])
round(res_ct,4)
#> character(0) H1-
#> Inf 0.1209
#> H2+ c("H1-", "H2+")
#> 0.2039 0.1209
#> H3+ c("H1-", "H3+")
#> 0.4451 0.1209
#> c("H2+", "H3+") c("H1-", "H2+", "H3+")
#> 0.2039 0.1209
#> H4+ c("H1-", "H4+")
#> 0.0271 0.0271
#> c("H2+", "H4+") c("H1-", "H2+", "H4+")
#> 0.0271 0.0271
#> c("H3+", "H4+") c("H1-", "H3+", "H4+")
#> 0.0271 0.0271
#> c("H2+", "H3+", "H4+") c("H1-", "H2+", "H3+", "H4+")
#> 0.0271 0.0271
Partitioning procedure with adaptive tests:
Sc = rep("-",n)
Sc[p_right_tail <= 0.5] <- "+"
hyp_part <- powset(1:n)
for (i in 1:length(hyp_part)){
hyp_part[[i]] <- Sc
hyp_part[[i]][ powset(1:n)[[i]] ] <- S[ powset(1:n)[[i]] ]
}
res_part <- res
names(res_part) <- hyp_part
round(res_part,4)
#> c("+", "-", "-", "-") c("-", "-", "-", "-") c("+", "+", "-", "-")
#> Inf 0.0403 0.1020
#> c("-", "+", "-", "-") c("+", "-", "+", "-") c("-", "-", "+", "-")
#> 0.0806 0.4451 0.0806
#> c("+", "+", "+", "-") c("-", "+", "+", "-") c("+", "-", "-", "+")
#> 0.2039 0.1209 0.0068
#> c("-", "-", "-", "+") c("+", "+", "-", "+") c("-", "+", "-", "+")
#> 0.0135 0.0135 0.0203
#> c("+", "-", "+", "+") c("-", "-", "+", "+") c("+", "+", "+", "+")
#> 0.0135 0.0203 0.0203
#> c("-", "+", "+", "+")
#> 0.0271
\(95\%\) confidence lower and upper
bounds \(\bar{\ell}_{\alpha}^{+}\) and
\(n - \bar{\ell}_{\alpha}^{-}\) for
\(n^+\) by using the partitioning
procedure with adaptive Simes’ tests:
library(nplus)
nplus_bound(p_right_tail, alpha=0.05, method="Simes")
#> $lo
#> [1] 1
#>
#> $up
#> [1] 3
Index set of positive discoveries (i.e. \(i:\theta_i >0\)) and non-positive
discoveries (\(i:\theta_i \leq 0\))
with familywise error rate control at level \(\alpha\):
nplus_fwer(p_right_tail, alpha=0.05, method="Simes")
#> $positive
#> integer(0)
#>
#> $nonpositive
#> [1] 4
Gail, M. and Simon, R. (1985). Testing for qualitative interactions
between treatment effects and patient subsets. Biometrics, pages
361–372.
Testing for qualitative interactions with follow-up inference
If testing for mixed signs of the parameter vector \(\theta\) is of primary concern, DCT is
applied only if the hypothesis of no qualitative interaction, \(H_0:\{n^+=0\}\cup\{ n^−=0\}\), is
rejected.
res_qi = res
powset_ix = powset(1:n)
p_Sminus =res[sapply(1:length(powset_ix), function(i) identical(which(S=="-") , unlist(powset_ix[i])))]
p_Splus =res[sapply(1:length(powset_ix), function(i) identical(which(S=="+") , unlist(powset_ix[i])))]
p_H0 = max(p_Sminus,p_Splus)
round(p_H0,4)
#> [1] 0.0403
res_qi[sapply(1:length(powset_ix), function(i) all(which(S=="+") %in% unlist(powset_ix[i])))] <- p_H0
res_qi[sapply(1:length(powset_ix), function(i) all(which(S=="-") %in% unlist(powset_ix[i])))] <- p_H0
round(res_qi,4)
#> character(0) H1-
#> Inf 0.0403
#> H2+ c("H1-", "H2+")
#> 0.1020 0.0403
#> H3+ c("H1-", "H3+")
#> 0.4451 0.0403
#> c("H2+", "H3+") c("H1-", "H2+", "H3+")
#> 0.2039 0.0403
#> H4+ c("H1-", "H4+")
#> 0.0068 0.0403
#> c("H2+", "H4+") c("H1-", "H2+", "H4+")
#> 0.0135 0.0403
#> c("H3+", "H4+") c("H1-", "H3+", "H4+")
#> 0.0135 0.0403
#> c("H2+", "H3+", "H4+") c("H1-", "H2+", "H3+", "H4+")
#> 0.0403 0.0403
res_qi_ct <- res_qi
for (k in 1:length(closure))
res_qi_ct[[k]] = max(res_qi_ct[sapply(1:length(closure), function(i) all(closure[[k]] %in% closure[[i]]))])
round(res_qi_ct,4)
#> character(0) H1-
#> Inf 0.0403
#> H2+ c("H1-", "H2+")
#> 0.2039 0.0403
#> H3+ c("H1-", "H3+")
#> 0.4451 0.0403
#> c("H2+", "H3+") c("H1-", "H2+", "H3+")
#> 0.2039 0.0403
#> H4+ c("H1-", "H4+")
#> 0.0403 0.0403
#> c("H2+", "H4+") c("H1-", "H2+", "H4+")
#> 0.0403 0.0403
#> c("H3+", "H4+") c("H1-", "H3+", "H4+")
#> 0.0403 0.0403
#> c("H2+", "H3+", "H4+") c("H1-", "H2+", "H3+", "H4+")
#> 0.0403 0.0403
Enhancing meta-analysis
We consider the data in Konstantopoulos (2011), where the effect of a
modified school calendar (with more frequent but shorter breaks) on
student achievement is examined. The meta-analysis uses studies of \(n=56\) schools in 11 districts. The
experimental design is a two-group comparison (modified calendar vs
traditional calendar) which involves computing the standardized mean
difference \(X_i \sim N(\theta_i,1)\),
where \(\theta_i>0\) indicates a
positive effects, i.e. a higher level of achievement in the group
following the modified school calendar.
We analyze the data by providing the lower and upper bounds on the
number of studies with positive effect in all schools and in the
(data-driven) top \(k\) schools.
Importing data from the R package metafor:
require("metafor") || install.packages("metafor")
dat <- get(data(dat.konstantopoulos2011))
x <- as.vector(dat$yi / sqrt(dat$vi))
n <- length(x)
p_right_tail <- pnorm(x, lower.tail = FALSE)
Partitioning procedure with adaptive tests (P) and directional closed
testing (DCT) procedures: \((1-\alpha)=95\%\) confidence bounds
[LO+,UP+] on the number of studies with
positive effects along with the number of positive discoveries
(D+) and \(95\%\)
confidence bounds [LO-,UP-] on the number of
studies with nonpositive/negative effects along with the number of
nonpositive/negative discoveries (D-) for different
combining functions: Fisher, Simes, modified Simes and adaptive LRT:
alpha = 0.05
S_minus = which(p_right_tail <= 0.5)
S_plus = which(p_right_tail > 0.5)
RES = matrix(NA, ncol=8, nrow=8)
RES[,1] = rep(c("Fisher","Simes","ASimes","ALRT"),each=2)
RES[,2] = rep(c("P","DCT"),4)
colnames(RES) <- c("Local test","Procedure", "LO+","UP+","D+", "LO-","UP-","D-")
for (i in 1:4){
res_bound <- nplus_bound(p_right_tail, method=RES[2*i,1], alpha=alpha)
res_d <- nplus_fwer(p_right_tail, method=RES[2*i,1], alpha=alpha)
res_dct_lo <- min(which(nplus_pvalue(p_right_tail, method=RES[2*i,1], ix = S_minus) > alpha)) - 1
res_dct_up <- n - (length(S_plus) - max(which(nplus_pvalue(p_right_tail, method=RES[2*i,1], ix = S_plus) > alpha)) + 1)
RES[(2*i)-1,3] <- res_bound$lo
RES[2*i,3] <- res_dct_lo
RES[(2*i)-1,4] <- res_bound$up
RES[2*i,4] <- res_dct_up
RES[(2*i)-1,5] <- RES[2*i,5] <- length(res_d$positive)
RES[(2*i)-1,6] <- n-res_bound$up
RES[2*i,6] <- n - res_dct_up
RES[(2*i)-1,7] <- n-res_bound$lo
RES[2*i,7] <- n - res_dct_lo
RES[(2*i)-1,8] <- RES[2*i,8] <- length(res_d$nonpositive)
}
knitr::kable(RES)
|
Local test
|
Procedure
|
LO+
|
UP+
|
D+
|
LO-
|
UP-
|
D-
|
|
Fisher
|
P
|
10
|
53
|
4
|
3
|
46
|
0
|
|
Fisher
|
DCT
|
9
|
56
|
4
|
0
|
47
|
0
|
|
Simes
|
P
|
8
|
55
|
8
|
1
|
48
|
0
|
|
Simes
|
DCT
|
8
|
56
|
8
|
0
|
48
|
0
|
|
ASimes
|
P
|
8
|
54
|
8
|
2
|
48
|
0
|
|
ASimes
|
DCT
|
8
|
56
|
8
|
0
|
48
|
0
|
|
ALRT
|
P
|
11
|
53
|
5
|
3
|
45
|
0
|
|
ALRT
|
DCT
|
9
|
55
|
5
|
1
|
47
|
0
|
Discoveries by Guo and Romano (2015) procedures with FWER and FDR
control at \(\alpha=0.05\) level:
GuoRomano2015 <- function(p_right_tail, alpha=0.05){
p <- p_right_tail
n <- length(p)
sp <- sort( c(p,1-p), index.return=TRUE )
FWER_id <- max(which(cumsum( sp$x <= alpha / (n - 1:n + 1L + alpha) )==1:(2*n)))
FWER_REJ <- rep(FALSE,2*n)
if (is.finite(FWER_id)) FWER_REJ[1:FWER_id] <- TRUE
FWER_NONPOS <- 1:n %in% sp$ix[FWER_REJ]
FWER_POS <- (n+1):(2*n) %in% sp$ix[FWER_REJ]
FDR_id <- max(which(sp$x <= alpha * (1:n) / n))
FDR_REJ <- rep(FALSE,2*n)
if (is.finite(FDR_id)) FDR_REJ[1:FDR_id] <- TRUE
FDR_NONPOS <- 1:n %in% sp$ix[FDR_REJ]
FDR_POS <- (n+1):(2*n) %in% sp$ix[FDR_REJ]
return(list(FWER = list(n_discoveries_positive = sum(FWER_NONPOS), n_discoveries_nonpositive = sum(FWER_POS)),
FDR= list(n_discoveries_positive = sum(FDR_NONPOS), n_discoveries_nonpositive = sum(FDR_POS))))
}
GuoRomano2015(p_right_tail, alpha = alpha)
#> $FWER
#> $FWER$n_discoveries_positive
#> [1] 8
#>
#> $FWER$n_discoveries_nonpositive
#> [1] 1
#>
#>
#> $FDR
#> $FDR$n_discoveries_positive
#> [1] 12
#>
#> $FDR$n_discoveries_nonpositive
#> [1] 4
Adaptive local test with partitioning: \(95\%\) post-hoc bounds for \(n^+(I_k)\) with Fisher’s combining
function, where \(I_k\) are the indexes
of the top \(k=|I_k|\) schools with
largest (in absolute value) effects. For example, among the top 38
schools, we observe at least 9 positive effects and at least 1
non-positive effect, or among the top 51 schools, we observe at least 10
positive effects and at least 2 non-positive effects, and so on.
ix_best <- sort(pmin(p_right_tail,1-p_right_tail), index.return=TRUE)$ix
bestk <- sapply(1:n, function(i)
range(which( nplus_pvalue(p_right_tail, method = "Fisher", ix = ix_best[1:i]) > 0.05)) -1
)
plot(asp=1, type="s", ylim=c(0,n),xlim=c(1,n),
bestk[2,], xlab="k", ylab="95% CI for the top k schools")
lines(bestk[1,],type="s")
abline(h=0)
abline(a=0,b=1)
abline(v=n)
for (i in 1:n){
segments(x0=i,x1=i,y0=bestk[1,i], y1=bestk[2,i])
}
segments(x0=37,x1=37,y0=bestk[1,37], y1=bestk[2,37], col="red")
segments(x0=50,x1=50,y0=bestk[1,50], y1=bestk[2,50], col="blue")
Konstantopoulos, S. (2011). Fixed effects and variance components
estimation in three-level meta-analysis. Research Synthesis Methods,
2(1):61–76.
Guo, W. and Romano, J. P. (2015). On stepwise control of directional
errors under independence and some dependence. Journal of Statistical
Planning and Inference, 163:21–33.
Power for \(n=2\)
Reproducing Fig. 2 in Heller & Solari (2023):
The partitioning procedure with Dunnett’s local tests for
many-to-one comparisons
To illustrate this approach, we consider the recovery data from Bretz
et al. (2016). A company conducted a study to evaluate the effectiveness
of specialized heating blankets in assisting post-surgical body heat.
Four types of blankets, labeled b0, b1, b2, and b3, were tested on
surgical patients to measure their impact on recovery times. The b0
blanket served as the standard option already utilized in different
hospitals. The main focus was to assess the recovery time, measured in
minutes, of patients randomly assigned to one of the four treatments. A
shorter recovery time would indicate a more effective treatment. The key
question addressed in this study is whether blanket types b1, b2, or b3
modify recovery time compared to b0. To perform a formal analysis of the
data, we assume an one-way layout \[y_{ij} =
\gamma + \mu_i + \epsilon_{ij}\] where \(\gamma + \mu_i\) represents the expected
recovery time for blanket b\(i\), \(i=0,1,2,3\), with i.i.d. errors \(\epsilon_{ij} \sim N(0,\sigma^2)\).
rm(list=ls())
library(mvtnorm)
library(multcomp)
#> Loading required package: survival
#> Loading required package: TH.data
#> Loading required package: MASS
#>
#> Attaching package: 'TH.data'
#> The following object is masked from 'package:MASS':
#>
#> geyser
data("recovery", package = "multcomp")
plot(minutes ~ blanket, data = recovery)
fit <- aov(minutes ~ blanket, data = recovery)
sigma <- sqrt( sum(fit$residuals^2)/fit$df.residual )
Dunnett’s procedure is the standard procedure used to address the
many-to-one comparisons problem. In the absence of specific information
about the direction of the effect, we typically formulate two-sided null
hypotheses of the form \(H_i: \theta_i =
0\), where \(\theta_i = \mu_i -
\mu_0\) and \(i=1,2,3\).
set.seed(123)
# Dunnett two-sided single step
dun_twosided <- glht(fit,
linfct = mcp(blanket = "Dunnett"),
alternative = "two.sided")
summary(dun_twosided)
#>
#> Simultaneous Tests for General Linear Hypotheses
#>
#> Multiple Comparisons of Means: Dunnett Contrasts
#>
#>
#> Fit: aov(formula = minutes ~ blanket, data = recovery)
#>
#> Linear Hypotheses:
#> Estimate Std. Error t value Pr(>|t|)
#> b1 - b0 == 0 -2.1333 1.6038 -1.330 0.456
#> b2 - b0 == 0 -7.4667 1.6038 -4.656 <0.001 ***
#> b3 - b0 == 0 -1.6667 0.8848 -1.884 0.182
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> (Adjusted p values reported -- single-step method)
plot(confint(dun_twosided, level=0.9))
Step-down Dunnett procedure for two-sided hypotheses. Although the
step-down Dunnett procedure is more powerful, the single-step Dunnett’s
procedure provides simultaneous confidence intervals.
# Dunnett two-sided step-down
summary(dun_twosided, test = adjusted(type = "free"))
#>
#> Simultaneous Tests for General Linear Hypotheses
#>
#> Multiple Comparisons of Means: Dunnett Contrasts
#>
#>
#> Fit: aov(formula = minutes ~ blanket, data = recovery)
#>
#> Linear Hypotheses:
#> Estimate Std. Error t value Pr(>|t|)
#> b1 - b0 == 0 -2.1333 1.6038 -1.330 0.19160
#> b2 - b0 == 0 -7.4667 1.6038 -4.656 0.00012 ***
#> b3 - b0 == 0 -1.6667 0.8848 -1.884 0.12729
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> (Adjusted p values reported -- free method)
Dunnett’s local test for the orthant hypothesis \(J_K: \{\bigcap_{i \in K} K_i\} \cap \{\bigcap_{i
\notin K}H^-_i\}\) is given by \[\psi_{K} = 1\{ \max(\{\hat{\theta}_i, i\in
K^c \}\cup \{-\hat{\theta}_j, j\in K \}) > c_{\alpha}(K)
\}\] where \(\hat{\theta}_i\) is
the estimate of \(\theta_i\) and \(c_{\alpha}( K)\) represents the \((1-\alpha)\)-quantile of the maximum of an
\(n\)-variate \(t\) distribution with a correlation matrix
that has off-diagonal elements whose sign depends on \(K\).
For example, the Dunnett’s local test for the orthant hypothesis
\(J_K\) with \(K=\{1,2\}\):
# Contrast matrix
contr <- rbind("b1 -b0" = c(-1, 1, 0, 0),
"b2 -b0" = c(-1, 0, 1, 0),
"b3 -b0" = c(1, 0, 0, -1))
# Index set K
as.numeric(which(contr[,1]==-1))
#> [1] 1 2
RES <- summary(glht(fit, linfct = mcp(blanket = contr),alternative = "less"))
# P-value
min(RES$test$pvalues)
#> [1] 6.055887e-05
# Correlation matrix
VC <- RES$lin %*% RES$vcov %*% t(RES$lin)
D <- diag(1/sqrt(diag(VC)))
COR <- D %*% VC %*% D
COR
#> [,1] [,2] [,3]
#> [1,] 1.0000000 0.1304348 -0.2364331
#> [2,] 0.1304348 1.0000000 -0.2364331
#> [3,] -0.2364331 -0.2364331 1.0000000
# Critical value
alpha = 0.1
qmvt(1-alpha, corr = COR, df = RES$df, tail = "lower.tail")$quantile
#> [1] 1.873676
# Test statistic
max(-RES$test$tstat)
#> [1] 4.655648
The partitioning procedure with Dunnett’s local tests:
|
K
|
Statistic
|
Critical value
|
P-value
|
|
c(1, 2, 3)
|
4.656
|
1.843
|
0
|
|
c(1, 2)
|
4.656
|
1.873
|
0
|
|
c(1, 3)
|
1.884
|
1.867
|
0.097
|
|
c(1, 3)
|
1.884
|
1.867
|
0.097
|
|
1
|
1.33
|
1.867
|
0.259
|
|
2
|
4.656
|
1.866
|
0
|
|
3
|
1.884
|
1.874
|
0.098
|
|
numeric(0)
|
-1.33
|
1.843
|
0.996
|
Adjusted \(p\)-values for base
hypotheses:
|
Hypothesis
|
Adjusted p-value
|
Hypothesis
|
Adjusted p-value
|
|
H1-
|
0.996
|
K1
|
0.259
|
|
H2-
|
0.996
|
K2
|
0
|
|
H3-
|
0.996
|
K3
|
0.098
|