R help - May 2004 - pmvt problem in multcomp

Hi, all:

Two examples are shown below. 

I want to use the multiple comparison of Dunnett.
It succeeded in upper case "example 1". 

However, the lower case "example 2" went wrong. 

In "example 2", the function pmvt return NaN, so I cannot show
this simtest result. Is there any solution?

(I changed the variable "maxpts" to a large number in front of
the function pmvt ... but, the function mvt returned an error. )

-- example 1 -------------------------------
require(multcomp)
Loading required package: multcomp 
Loading required package: mvtnorm 
[1] TRUE

y <- as.vector(t$int)
f <- as.factor(t$group1)

table(f)
f
    1     2     3 
20988 20988 20988 

dat <- cbind(as.data.frame(y),f)
gc()
summary(simtest(y ~ f, data=dat, type="Dunnett"))

	 Simultaneous tests: Dunnett contrasts 

Call: 
simtest.formula(formula = y ~ f, data = dat, type = "Dunnett")

	 Dunnett contrasts for factor f

Contrast matrix:
        f1 f2 f3
f2-f1 0 -1  1  0
f3-f1 0 -1  0  1


Absolute Error Tolerance:  0.001 

Coefficients:
      Estimate t value Std.Err. p raw p Bonf p adj
f2-f1    4.015  -0.677    5.934 0.499  0.997 0.722
f3-f1    2.486  -0.419    5.934 0.675  0.997 0.722
---------------------------------

-- example 2 -------------------------------
require(multcomp)
Loading required package: multcomp 
Loading required package: mvtnorm 
[1] TRUE

y <- as.vector(t$int)
f <- as.factor(t$group2)
table(f)
f
     1      2      3      4      5 
104940 104940 104940 104940 104940 

dat <- cbind(as.data.frame(y),f)
gc()
summary(simtest(y ~ f, data=dat, type="Dunnett"))

[1] "des <- model.matrix(ff, mf)"
  (Intercept) aaa1 aaa2 aaa3 aaa4
1           1    1    0    0    0
2           1    0    1    0    0
3           1    0    0    1    0
4           1    0    0    0    1
attr(,"assign")
[1] 0 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$aaa
[1] "ct"

[1] "gls     <- rep(0,nrow(contonly))"
[1] 0 0 0 0

[1] "gls[i1]    <- 1-prob"
[1] NaN   0   0   0
[1] "gls[i1]    <- 1-prob"
[1] NaN NaN   0   0
[1] "gls[i1]    <- 1-prob"
[1]          NaN          NaN -7.01661e-14  0.00000e+00
[1] "gls[i1]    <- 1-prob"
[1]           NaN           NaN -7.016610e-14  3.362133e-11

[1] "glsbig"
  aaa1 aaa2         aaa3         aaa4
1  NaN  NaN          NaN          NaN
2  NaN  NaN          NaN          NaN
3    0    0 -7.01661e-14 0.000000e+00
4    0    0  0.00000e+00 3.362133e-11

[1] "glsp"
aaa1 aaa2 aaa3 aaa4 
 NaN  NaN  NaN  NaN 
Error in if (glsp[i] < glsp[i - 1]) { : missing value where TRUE/FALSE needed
Execution halted
---------------------------------

-- 
kuroki at oak.dti.ne.jp
GnuPG fingerprint = 90FD FE79 905F 26F9 29C4  096F 8AA2 2C42 5130 1469

Yes, `pmvt' returns NaN without indicating this error. We need to check.
Thanks for the report (and *please* cc emails reporting problems with
packages to the maintainer!),

Torsten

On Thu, 20 May 2004, Chihiro Kuroki wrote:
> Hi, all:
>
> Two examples are shown below.
>
> I want to use the multiple comparison of Dunnett.
> It succeeded in upper case "example 1".
>
> However, the lower case "example 2" went wrong.
>
> In "example 2", the function pmvt return NaN, so I cannot show
> this simtest result. Is there any solution?
>
> (I changed the variable "maxpts" to a large number in front of
> the function pmvt ... but, the function mvt returned an error. )
>
> -- example 1 -------------------------------
> require(multcomp)
> Loading required package: multcomp
> Loading required package: mvtnorm
> [1] TRUE
>
> y <- as.vector(t$int)
> f <- as.factor(t$group1)
>
> table(f)
> f
>     1     2     3
> 20988 20988 20988
>
> dat <- cbind(as.data.frame(y),f)
> gc()
> summary(simtest(y ~ f, data=dat, type="Dunnett"))
>
> 	 Simultaneous tests: Dunnett contrasts
>
> Call:
> simtest.formula(formula = y ~ f, data = dat, type = "Dunnett")
>
> 	 Dunnett contrasts for factor f
>
> Contrast matrix:
>         f1 f2 f3
> f2-f1 0 -1  1  0
> f3-f1 0 -1  0  1
>
>
> Absolute Error Tolerance:  0.001
>
> Coefficients:
>       Estimate t value Std.Err. p raw p Bonf p adj
> f2-f1    4.015  -0.677    5.934 0.499  0.997 0.722
> f3-f1    2.486  -0.419    5.934 0.675  0.997 0.722
> ---------------------------------
>
> -- example 2 -------------------------------
> require(multcomp)
> Loading required package: multcomp
> Loading required package: mvtnorm
> [1] TRUE
>
> y <- as.vector(t$int)
> f <- as.factor(t$group2)
> table(f)
> f
>      1      2      3      4      5
> 104940 104940 104940 104940 104940
>
> dat <- cbind(as.data.frame(y),f)
> gc()
> summary(simtest(y ~ f, data=dat, type="Dunnett"))
>
> [1] "des <- model.matrix(ff, mf)"
>   (Intercept) aaa1 aaa2 aaa3 aaa4
> 1           1    1    0    0    0
> 2           1    0    1    0    0
> 3           1    0    0    1    0
> 4           1    0    0    0    1
> attr(,"assign")
> [1] 0 1 1 1 1
> attr(,"contrasts")
> attr(,"contrasts")$aaa
> [1] "ct"
>
> [1] "gls     <- rep(0,nrow(contonly))"
> [1] 0 0 0 0
>
> [1] "gls[i1]    <- 1-prob"
> [1] NaN   0   0   0
> [1] "gls[i1]    <- 1-prob"
> [1] NaN NaN   0   0
> [1] "gls[i1]    <- 1-prob"
> [1]          NaN          NaN -7.01661e-14  0.00000e+00
> [1] "gls[i1]    <- 1-prob"
> [1]           NaN           NaN -7.016610e-14  3.362133e-11
>
> [1] "glsbig"
>   aaa1 aaa2         aaa3         aaa4
> 1  NaN  NaN          NaN          NaN
> 2  NaN  NaN          NaN          NaN
> 3    0    0 -7.01661e-14 0.000000e+00
> 4    0    0  0.00000e+00 3.362133e-11
>
> [1] "glsp"
> aaa1 aaa2 aaa3 aaa4
>  NaN  NaN  NaN  NaN
> Error in if (glsp[i] < glsp[i - 1]) { : missing value where TRUE/FALSE
needed
> Execution halted
> ---------------------------------
>
> --
> kuroki at oak.dti.ne.jp
> GnuPG fingerprint = 90FD FE79 905F 26F9 29C4  096F 8AA2 2C42 5130 1469
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
>
>

On Thu, 20 May 2004, Chihiro Kuroki wrote:
> Hi, all:
>
> Two examples are shown below.
>
> I want to use the multiple comparison of Dunnett.
> It succeeded in upper case "example 1".
>
> However, the lower case "example 2" went wrong.
>
it was due to an error in the code underlying the `mvtnorm' package. The
problem was fixed by Alan Genz and I uploaded a revised version (0.6-7) of
the package to CRAN a few minutes ago.

Best,

Torsten

> In "example 2", the function pmvt return NaN, so I cannot show
> this simtest result. Is there any solution?
>
> (I changed the variable "maxpts" to a large number in front of
> the function pmvt ... but, the function mvt returned an error. )
>
> -- example 1 -------------------------------
> require(multcomp)
> Loading required package: multcomp
> Loading required package: mvtnorm
> [1] TRUE
>
> y <- as.vector(t$int)
> f <- as.factor(t$group1)
>
> table(f)
> f
>     1     2     3
> 20988 20988 20988
>
> dat <- cbind(as.data.frame(y),f)
> gc()
> summary(simtest(y ~ f, data=dat, type="Dunnett"))
>
> 	 Simultaneous tests: Dunnett contrasts
>
> Call:
> simtest.formula(formula = y ~ f, data = dat, type = "Dunnett")
>
> 	 Dunnett contrasts for factor f
>
> Contrast matrix:
>         f1 f2 f3
> f2-f1 0 -1  1  0
> f3-f1 0 -1  0  1
>
>
> Absolute Error Tolerance:  0.001
>
> Coefficients:
>       Estimate t value Std.Err. p raw p Bonf p adj
> f2-f1    4.015  -0.677    5.934 0.499  0.997 0.722
> f3-f1    2.486  -0.419    5.934 0.675  0.997 0.722
> ---------------------------------
>
> -- example 2 -------------------------------
> require(multcomp)
> Loading required package: multcomp
> Loading required package: mvtnorm
> [1] TRUE
>
> y <- as.vector(t$int)
> f <- as.factor(t$group2)
> table(f)
> f
>      1      2      3      4      5
> 104940 104940 104940 104940 104940
>
> dat <- cbind(as.data.frame(y),f)
> gc()
> summary(simtest(y ~ f, data=dat, type="Dunnett"))
>
> [1] "des <- model.matrix(ff, mf)"
>   (Intercept) aaa1 aaa2 aaa3 aaa4
> 1           1    1    0    0    0
> 2           1    0    1    0    0
> 3           1    0    0    1    0
> 4           1    0    0    0    1
> attr(,"assign")
> [1] 0 1 1 1 1
> attr(,"contrasts")
> attr(,"contrasts")$aaa
> [1] "ct"
>
> [1] "gls     <- rep(0,nrow(contonly))"
> [1] 0 0 0 0
>
> [1] "gls[i1]    <- 1-prob"
> [1] NaN   0   0   0
> [1] "gls[i1]    <- 1-prob"
> [1] NaN NaN   0   0
> [1] "gls[i1]    <- 1-prob"
> [1]          NaN          NaN -7.01661e-14  0.00000e+00
> [1] "gls[i1]    <- 1-prob"
> [1]           NaN           NaN -7.016610e-14  3.362133e-11
>
> [1] "glsbig"
>   aaa1 aaa2         aaa3         aaa4
> 1  NaN  NaN          NaN          NaN
> 2  NaN  NaN          NaN          NaN
> 3    0    0 -7.01661e-14 0.000000e+00
> 4    0    0  0.00000e+00 3.362133e-11
>
> [1] "glsp"
> aaa1 aaa2 aaa3 aaa4
>  NaN  NaN  NaN  NaN
> Error in if (glsp[i] < glsp[i - 1]) { : missing value where TRUE/FALSE
needed
> Execution halted
> ---------------------------------
>
> --
> kuroki at oak.dti.ne.jp
> GnuPG fingerprint = 90FD FE79 905F 26F9 29C4  096F 8AA2 2C42 5130 1469
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
>
>

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

Dear Mr.Torsten:

At Wed, 16 Jun 2004 08:04:08 +0200 (CEST),
Torsten Hothorn wrote:
> > 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
> > f22-f21    2.125  -2.079    1.113 0.045  0.045 0.045
> > ---------------------------------
> 
> Chihiro,
> 
> Frank and I used your data to check the program and example with an
> independent
> algorithm and implementation (Westfall-Young stepdown resampling
> procedure). Theory suggests that the
> results should be similar to (but not necessarily the same as) those
> obtained with multcomp in this special case. These are the adjusted
> p-values obtained with the Westfall-Young approach for 100,000
> replications:
> 
> which fit nicely with the ones obtained from multcomp.
> 
> Hope this helps & sorry for the delay,
Thank you.

At Mon, 07 Jun 2004 12:14:26 +0900,
myself-oak wrote:
> 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).
You said "adj p" values are ...
> 0.0437
> 0.0260
> 0.0204
> 0.0001
I used SPSS ver.10 and got the following result.

Dunnett t (two sided) 

 | --------------- | ---------- | -------- | -------- | ------------------ | 
 |                 | mean diff. | s.e.     | adj p    | 95% C.I.           | 
 | ------ | ------ | (I-J)      |          |          | ---------- | ----- | 
 | (I) V3 | (J) V3 |            |          |          |            |       | 
 | ------ | ------ | ---------- | -------- | -------- | ---------- | ----- | 
 | 2      | 1      | 2.125      | 1.022    | .145     | -.507      | 4.757 | 
 | ------ | ------ | ---------- | -------- | -------- | ---------- | ----- | 
 | 3      | 1      | 2.875(*)   | 1.022    | .028     | .243       | 5.507 | 
 | ------ | ------ | ---------- | -------- | -------- | ---------- | ----- | 
 | 4      | 1      | 2.625      | 1.022    | .051     | -7.325E-03 | 5.257 | 
 | ------ | ------ | ---------- | -------- | -------- | ---------- | ----- | 
 | 5      | 1      | 5.167(*)   | 1.113    | .000     | 2.301      | 8.032 | 
 | ------ | ------ | ---------- | -------- | -------- | ---------- | ----- | 

I might make a mistake in the way to use the simtest()...How should I think? 

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