Re: [R] pmvt problem in multcomp

About this list Date view Thread view Subject view Author view Attachment view

From: Torsten Hothorn (Torsten.Hothorn@rzmail.uni-erlangen.de)
Date: Thu 27 May 2004 - 20:00:44 EST


Message-id: <Pine.LNX.4.51.0405271158500.1100@artemis.imbe.med.uni-erlangen.de>

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@oak.dti.ne.jp
> GnuPG fingerprint = 90FD FE79 905F 26F9 29C4 096F 8AA2 2C42 5130 1469
>
> ______________________________________________
> R-help@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
>
>

______________________________________________
R-help@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


About this list Date view Thread view Subject view Author view Attachment view

This archive was generated by hypermail 2.1.3 : Mon 31 May 2004 - 23:05:13 EST