5
$\begingroup$

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

enter image description here

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!

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ Commented Jul 10, 2019 at 9:25

2 Answers 2

2
$\begingroup$

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)
})
# adjusted p-values
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)
$\endgroup$
1
  • $\begingroup$ Thank you very much, Stéphane! @Stéphane Laurent $\endgroup$
    – Lizzy
    Commented Jul 20, 2019 at 3:27
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.