How to calculate critical values for Dunnett procedure given alpha, df1 and df2 I am looking for an R function that can calculate critical values for a Dunnett multiple comparison procedure for any family-wise error rate and degrees of freedoms. 
In short, I need similar functions to the two below found in R to calculate quantiles for the fisher and studentized distribution:
f_crit <- qf(p, df1, df2, lower.tail=F)
q_crit <- qtukey(p, df1, df2, lower.tail=F)

Is there any R package that provide similar function for the critical values of Dunnett? 
d_crit <- ???(p, df1, df2)
 A: It turns out the nCDunnett package contain functions to calculate the quantiles similar to qf() and qtukey() above. 
For one-tailed distribution:
>require(nCDunnett)
>df1 = 5
>df2 = 7
>p = 0.01
>q = qNCDun(p=1-p, nu=df2, rho=(rep(0.5,times=df1)), delta=rep(0,times=df1), two.sided=F)

>q

3.958376
Giving similar output as the table values found in text-books, such as dunnett-1side.pdf. The qNCDun() function, however, operates an order of magnitude slower than qf() and qtukey().
A: qDunnett <- function (p, df, k, rho,
                      type = c("two-sided", "one-sided"))
{
  type <- match.arg(type)
  alpha <- 1 - p
  if (type == "two-sided") {
    alpha <- alpha/2
  }
  S <- matrix(rho, nrow=k, ncol=k) + (1-rho)*diag(k)
  if (type == "two-sided") {
    f <- function(d, df, k, S, p) {
      mnormt::sadmvt(df=df, lower=rep(-d,k), upper=rep(d,k),
                      mean=rep(0,k), S=S, maxpts=2000*k) - p
    }
  }
  else {
    f <- function(d, df, k, S, p) {
      mnormt::pmt(d, S=S, df=df) - p
    }
  }
  d <- uniroot(f,
               df = df, k = k, S = S, p=p,
               lower=qt(1 - alpha, df),
               upper=qt(1 - alpha/k, df),
               tol=.Machine$double.eps, maxiter=5000)$root
  return(d)
}

Test:
> p <- 0.95; df <- 24; rho <- 0.5; k <- 3
> nCDunnett::qNCDun(p=p, nu=df, rho=rho,
+                   delta=rep(0,times=k), two.sided=T)
[1] 2.506845
> qDunnett(p, df, k, rho)
[1] 2.506654

