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!
 A: 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.
