[R] pmvt problem in multcomp

[Chihiro Kuroki]

I found the following function in http://aoki2.si.gunma-u.ac.jp/R/dunnett.html.

require(mvtnorm)
dunnett <- function(dat, observe, group)
{
    printf <- function(fmt, ...) cat(sprintf(fmt, ...))

    get.rho <- function(ni)
    {
        k <- length(ni)
        rho <- outer(ni, ni, function(x, y) { sqrt(x/(x+ni[1])*y/(y+ni[1])) })
        diag(rho) <- 0
        sum(rho[-1,-1])/(k-2)/(k-1)
    }

    pdunnett <- function(a,df,x,r)
    {
        corr <- diag(a-1)
        corr[lower.tri(corr)] <- r
        1-pmvt(lower=-x, upper=x, delta=rep(0, a-1), df=df, corr=corr, abseps=0.0001)
    }

    x <- dat[,observe]
    ni <- table(g <- dat[,group])
    gv <- as.numeric(names(ni))
    a <- length(ni)
    n <- sum(ni)
    mi <- vi <- rep(0, a)
    for (i in 1:a) {
        mi[i] <- mean(xi <- x[g==gv[i]])
        vi[i] <- var(xi)
    }
    Vw <- sum(vi*(ni-1))/(n-a)
    rho <- get.rho(ni)
    printf("rho=%5.3f\n", rho)
    for (i in 2:a) {
        ti <- abs(mi[i]-mi[1])/sqrt(Vw*(1/ni[i]+1/ni[1]))
        pi <- pdunnett(a, n-a, ti, rho)
        printf("group:%i t=%.3f p=%.3f\n", i, ti, pi[1])
    }
}

> > > summary(simtest(y4 ~ f2, data=dat2, type="Dunnett"))
> >
> > 	 Simultaneous tests: Dunnett contrasts
> >
> > Call:
> > simtest.formula(formula = y4 ~ f2, data = dat2, type = "Dunnett")
> >
> > 	 Dunnett contrasts for factor f2
> >
> > Contrast matrix:
> >           f21 f22 f23 f24 f25
> > f22-f21 0  -1   1   0   0   0
> > f23-f21 0  -1   0   1   0   0
> > f24-f21 0  -1   0   0   1   0
> > f25-f21 0  -1   0   0   0   1
> >
> >
> > Absolute Error Tolerance:  0.001
> >
> > Coefficients:
> >         Estimate t value Std.Err. p raw p Bonf p adj
> > f25-f21    5.167  -4.644    1.022 0.000  0.000 0.000
> > f23-f21    2.875  -2.813    1.022 0.008  0.024 0.022
> > f24-f21    2.625  -2.569    1.022 0.015  0.029 0.028 --- (A)
> > f22-f21    2.125  -2.079    1.113 0.045  0.045 0.045
> > ---------------------------------

dunnett(dat2,1,2)
rho=0.426
group:5 t=4.644 p=0.000
group:3 t=2.813 p=0.028
group:4 t=2.569 p=0.051 --- (B)
group:2 t=2.079 p=0.145
(sorted in order of p values.)

p values are different although t values are equal.

> > I got the following inequality from the appended chart of a
> > book.
> >
> 
> hm, without knowing what
> 
> > 2.558 < d(5, 35, 0.4263464, 0.05) < 2.598
> 
> means it is hard to tell what the problem is. Could you please explain it
> further?

The alternative hypothesis is "two sided".

When significant level is equal to 0.05 , number of groups=5,
df of error=35 and rho=0.426, I think that absolute t-value
should be between 2.558 and 2.598.

So, (B) is easy to understand for me than (A).
In (A), I cannot believe easily that the "p adj" value is
smaller than "p Bonf".
And I want to know why the above results(A,B) are different.

Best regards,
-- 
kuroki
GnuPG fingerprint = 90FD FE79 905F 26F9 29C4  096F 8AA2 2C42 5130 1469