[R] pmvt problem in multcomp

[Torsten Hothorn]

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
>
>