2
$\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$ Jul 10 '19 at 9:25
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
    Jul 20 '19 at 3:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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