What's the formula for the 'Adjusted p values reported -- single-step method' in glht function?

For example, in the tutorial,

data('recovery', package = 'multcomp')

summary(recovery)

shows one control group b0 and three treatment groups b1~ b3. Then

mod <- aov(minutes ~ blanket, data = recovery)

library('multcomp')

mc <- glht(mod, linfct = mcp(blanket = 'Dunnett'), alternative = 'less')

summary(mc)

give

My question is, how the last column of Pr(< t)is calculated?

To be specific, the t value for b1-b0>=0 is -1.330, then why this Pr(< t) is 0.2411? What is the adjusted p value formula here?

Thanks!

• Do you know what p-value adjustment is? It is done when doing multiple comparisons, as you are doing here. In your particular case, you are using Dunnett correction for p-values. Jul 10, 2019 at 9:25

Consider the ANOVA statistical model $$y_{ij} = \gamma + \mu_i + \sigma\epsilon_{ij}, \quad i = 0, 1, \ldots, m, \quad j = 1, \ldots, n_i.$$ where $$\epsilon_{ij} \sim_{\textit{iid}} \mathcal{N}(0,1)$$.

The number $$\gamma$$ is the overall mean and $$\mu_i$$ is the effect in group $$i$$. Group $$0$$ is considered to be a control group, and we are interested in the $$m$$ hypotheses $$H_{0i} = \{\mu_0 = \mu_i\} \quad \textit{vs} \quad H_{1i} = \{\mu_0 > \mu_i\}.$$ Let $$H_0 = \bigcap_{i=1}^m H_{0i} = \{\mu_0=\mu_1=\cdots \mu_m\}$$.

Consider the $$m$$ mean differences $$\bar{y_i} - \bar{y_0} = (\mu_i - \mu_0) + \sigma(\bar{\epsilon_i}-\bar{\epsilon_0})$$ for $$i = 1, \ldots, m$$.

Hereafter we assume the null hypothesis $$H_0$$ is fulfilled, therefore $$\bar{y_i} - \bar{y_0} = \sigma(\bar{\epsilon_i}-\bar{\epsilon_0})$$.

We will derive the covariances between the standardized mean differences. Let's start by deriving the covariances between the $$\bar{\epsilon_i}-\bar{\epsilon_0}$$.

The $$\bar{\epsilon_i}$$ for $$i = 0, 1, \ldots, m$$ are independent and $$\bar{\epsilon_i} \sim \mathcal{N}\left(0, \frac{1}{n_i}\right)$$. Consequently, for $$i = 1, \ldots, m$$,
$$\bar{\epsilon_i}-\bar{\epsilon_0} \sim \mathcal{N}\left(0, \frac{1}{n_i} + \frac{1}{n_0}\right)$$ Now, let's derive the covariance between the $$\bar{\epsilon_i}-\bar{\epsilon_0}$$. One has \begin{align*} \textrm{Cov}(\bar{\epsilon_i}-\bar{\epsilon_0}, \bar{\epsilon_j}-\bar{\epsilon_0}) & = E\bigl((\bar{\epsilon_i}-\bar{\epsilon_0})(\bar{\epsilon_j}-\bar{\epsilon_0})\bigr) \\ & = E(\bar{\epsilon_i}\bar{\epsilon_j}) + E(\bar{\epsilon_0}\bar{\epsilon_0}) - E(\bar{\epsilon_i}\bar{\epsilon_0}) - E(\bar{\epsilon_j}\bar{\epsilon_0}). \end{align*} For $$i\neq j$$, each term is equal to $$0$$ except the second one: $$E(\bar{\epsilon_0}\bar{\epsilon_0}) = \frac{1}{n_0}$$ Thus $$\textrm{Cov}(\bar{\epsilon_i}-\bar{\epsilon_0}, \bar{\epsilon_j}-\bar{\epsilon_0}) = \frac{1}{n_0}$$ for every $$i,j \in \{1, \ldots, m\}$$, when $$i \neq j$$.

Now, let's introduce the standardized mean differences $$\delta_i = \frac{\bar{y_i} - \bar{y_0}}{\sigma\sqrt{\frac{1}{n_i} + \frac{1}{n_0}}} \sim \mathcal{N}(0,1).$$ for $$i = 1, \ldots, m$$. For $$i \neq j$$, the covariances are \begin{align*} \textrm{Cov}(\delta_i, \delta_j) & = \frac{1}{n_0}\frac{1}{\sqrt{\frac{1}{n_i} + \frac{1}{n_0}}} \frac{1}{\sqrt{\frac{1}{n_j} + \frac{1}{n_0}}} \\ & = \frac{1}{n_0}\sqrt{\frac{n_0n_i}{n_i+n_0}} \sqrt{\frac{n_0n_j}{n_j+n_0}} \\ & = \sqrt{\frac{n_i}{n_i+n_0}}\sqrt{\frac{n_j}{n_j+n_0}}. \end{align*}

Now let's introduce the statistics $$t_i = \frac{\bar{y_i} - \bar{y_0}}{\hat\sigma\sqrt{\frac{1}{n_i} + \frac{1}{n_0}}} = \frac{1}{u}\frac{\bar{y_i} - \bar{y_0}}{\sigma\sqrt{\frac{1}{n_i} + \frac{1}{n_0}}}$$ where $$\hat\sigma^2 = \frac{\sum_{i=0}^m\sum_{j=1}^{n_i}{(y_{ij}-\bar{y}_i)^2}}{\nu}$$ with $$\nu = \sum_{i=0}^m n_i - (m+1)$$, and where we have set $$u = \frac{\hat\sigma}{\sigma}$$.

It is known that $$\nu u^2 = \nu\frac{{\hat{\sigma}}^2}{\sigma^2} \sim \chi^2_\nu$$

Denote by $$\mathbf{t}$$ the vector $${(t_1, \ldots, t_m)}'$$. By the above, one has $$(\mathbf{t} \mid u) \sim \mathcal{M}\mathcal{N}(\mathbf{0}, \Sigma/u^2)$$ where $$\Sigma_{ij} = \begin{cases} 1 & \text{if } i = j \\ \sqrt{\frac{n_i}{n_i+n_0}}\sqrt{\frac{n_j}{n_j+n_0}} & \text{if } i \neq j \end{cases}$$

Therefore, $$\mathbf{t}$$ follows the multivariate Student distribution with $$\nu$$ degrees of freedom and scale matrix $$\Sigma$$ (see the paper A short review of multivariate $$t$$-distribution by Kibria & Joarder).

The adjusted Dunnett $$p$$-value for $$H_{0i}$$ vs the "less" alternative hypothesis $$H_{1i}$$ is $$p_i = \Pr\left(\min_{i \in\{1, \ldots m\}}t_i < t_i^{\text{obs}}\right)$$ where $$\Pr$$ is the probability under $$H_0$$.

One has \begin{align*} \Pr\left(\min_{i \in\{1, \ldots m\}}t_i < q\right) & = \Pr\left(-\min_{i \in\{1, \ldots m\}}t_i > -q\right) \\ & = \Pr\left(\max_{i \in\{1, \ldots m\}}-t_i > -q\right) \\ & = \Pr\left(\max_{i \in\{1, \ldots m\}}t_i > -q\right), \end{align*} the last equality stemming from the symmetry of the (centered) multivariate Student distribution.

Let's write a R function computing $$\Pr\left(\max_{i \in\{1, \ldots m\}}t_i \leq q \right)$$.

pDunnett <- function(q, n0, ni){
if(q==Inf){
return(1)
}
if(q==-Inf){
return(0)
}
m <- length(ni)
Sigma <- matrix(NA_real_, nrow=m, ncol=m)
for(i in 1:(m-1)){
for(j in (i+1):m){
Sigma[i,j] <- sqrt(ni[i]*ni[j]/(n0+ni[i])/(n0+ni[j]))
}
}
Sigma[lower.tri(Sigma)] <- Sigma[upper.tri(Sigma)]
diag(Sigma) <- 1
nu <- n0 + sum(ni) - m - 1
mnormt::pmt(rep(q,m), mean=rep(0,m), S=Sigma, df=nu, maxpts=2000*m)
}


Now let's try an example.

data("recovery", package = "multcomp")
recovery
##    blanket minutes
## 1       b0      15
## 2       b0      13
## 3       b0      12
## 4       b0      16
## 5       b0      16
## 6       b0      17
## 7       b0      13
## 8       b0      13
## 9       b0      16
## 10      b0      17
## 11      b0      17
## 12      b0      19
## 13      b0      17
## 14      b0      15
## 15      b0      13
## 16      b0      12
## 17      b0      16
## 18      b0      10
## 19      b0      17
## 20      b0      12
## 21      b1      13
## 22      b1      16
## 23      b1       9
## 24      b2       5
## 25      b2       8
## 26      b2       9
## 27      b3      14
## 28      b3      16
## 29      b3      16
## 30      b3      12
## 31      b3       7
## 32      b3      12
## 33      b3      13
## 34      b3      13
## 35      b3       9
## 36      b3      16
## 37      b3      13
## 38      b3      18
## 39      b3      13
## 40      b3      12
## 41      b3      13
# samples
ys <- lapply(split(recovery, recovery$blanket), function(df){ df$minutes
})
# sample sizes
sizes <- lengths(ys)
ni <- sizes[-1]
# degrees of freedom
nu <- sum(sizes) - length(ni) - 1
# pooled variance estimate
s2 <- sum(sapply(ys, function(y){
sum((y - mean(y))^2)
})) / nu
# Student statistics
n0 <- sizes[1]
y0bar <- mean(ys[[1]])
yi <- ys[-1]
ti <- sapply(yi, function(y){
(mean(y) - y0bar) / sqrt(1/length(y)+1/n0) / sqrt(s2)
})
sapply(ti, function(t) 1 - pDunnett(q = -t, n0, ni))
##        b1.b0        b2.b0        b3.b0
## 0.2411790515 0.0000585819 0.0924376838


Let's compare with the multcomp package:

fit <- lm(minutes ~ blanket, data = recovery)
library(multcomp)
multcomps <- glht(fit, linfct = mcp(blanket = "Dunnett"),
alternative = "less")
summary(multcomps)
##
##   Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Dunnett Contrasts
##
##
## Fit: lm(formula = minutes ~ blanket, data = recovery)
##
## Linear Hypotheses:
##              Estimate Std. Error t value Pr(<t)
## b1 - b0 >= 0  -2.1333     1.6038  -1.330 0.2412
## b2 - b0 >= 0  -7.4667     1.6038  -4.656 <0.001 ***
## b3 - b0 >= 0  -1.6667     0.8848  -1.884 0.0925 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)


The adjusted Dunnett $$p$$-value for $$H_{0i}$$ vs the "greater" alternative hypothesis $$\{\mu_0 < \mu_i\}$$ is $$p_i = \Pr\left(\max_{i \in\{1, \ldots m\}}t_i > t_i^{\text{obs}}\right)$$

sapply(ti, function(t) 1 - pDunnett(q = t, n0, ni))
##     b1.b0     b2.b0     b3.b0
## 0.9958029 1.0000000 0.9995174
multcomps <- glht(fit, linfct = mcp(blanket = "Dunnett"),
alternative = "greater")
summary(multcomps)
##
##   Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Dunnett Contrasts
##
##
## Fit: lm(formula = minutes ~ blanket, data = recovery)
##
## Linear Hypotheses:
##              Estimate Std. Error t value Pr(>t)
## b1 - b0 <= 0  -2.1333     1.6038  -1.330  0.996
## b2 - b0 <= 0  -7.4667     1.6038  -4.656  1.000
## b3 - b0 <= 0  -1.6667     0.8848  -1.884  1.000
## (Adjusted p values reported -- single-step method)

• Thank you very much, Stéphane! @Stéphane Laurent Jul 20, 2019 at 3:27

As explained in the foundation for the multcomp package, "Simultaneous Inference in General Parametric Models" (Hothorn et. al, 2008), the single step method works when the model provides a asymptotically, normally distributed estimate vector

$$\hat\theta \sim \approx N(\theta, \Sigma),$$

and a consistent estimator of $$\Sigma$$, $$\hat\Sigma$$. Assume that one wants to test $$k$$ hypotheses that can all be written in form $$H_0^j: K_j\theta=m_j$$, where $$K_j$$ is a matrix with one row. Define then $$K$$ as the stacked $$K_j$$'s (think rbind function of R), the hypothesis $$H_0$$ as the hypothesis that all $$H_0^j$$ are true, $$m=(m_1, \dots, m_k)$$, and define then the test statistic vector

$$T = \text{diag}(K\hat\Sigma K^*)^{-1/2}(K\hat\theta-m) \stackrel{\text{Under }H_0}{\sim\approx} N(0,R).$$

The single-step adjusted p-value, which controls the probability of false-positive conclusion among the $$k$$ hypotheses, is

$$1-P(\max(|T|)\leq t)=1-\int_{-t}^t\cdots\int_{-t}^{t} f_{N(0,\hat R)}(t_1, \dots,t_k) \text{ d}t_1\dots\text{d}t_k.$$

where $$t$$ is an absolute test statistic value (representing a realized value of an entry of $$T$$) and $$f_{N(0,R)}$$ is the density function of a multivariate normal distribution. For many models, the asymptotic normal distribution can be replaced by an exact $$t$$-distribution, which the glht function will do.