From ripley at stats.ox.ac.uk  Thu Mar  9 08:53:45 2006
From: ripley at stats.ox.ac.uk (ripley@stats.ox.ac.uk)
Date: Thu,  9 Mar 2006 08:53:45 +0100 (CET)
Subject: [Rd] bugs in simtest (PR#8670)
Message-ID: <20060309075345.B087A3FCA3@slim.kubism.ku.dk>

There is no simtest in R!  It appears you are using contributed package 
multcomp.

Please do re-read the FAQ.  Bug reports in contributed packages should not 
be sent to the R-bugs repository, but to the package maintainer.

This report has been closed.


On Thu, 9 Mar 2006, rmh at temple.edu wrote:

> # R for Windows will not send your bug report automatically.
> # Please copy the bug report (after finishing it) to
> # your favorite email program and send it to
> #
> #       r-bugs at r-project.org
> #
> ######################################################
>
> This report is joint from Richard Heiberger <rmh at temple.edu>
> and Burt Holland <bholland at temple.edu>.
> Burt Holland is the coauthor of the paper that the ?ptukey
> documentation references.
>
>
> R was used to run an example in our elementary Stat course.  It was
> a one-way ANOVA, the factor `strategy' having 3 levels Price, Quality
> and Convenience.
>
> We issued the command
>
> summary(simint(sales ~ strategy, type="Tukey", data=Xm15.01s))
>
> and received the output
>
> Coefficients:
>                                    Estimate  2.5 % 97.5 % t value Std.Err.
> strategyPrice-strategyConvenience      31.10 -40.67 102.87   1.043   29.824
> strategyQuality-strategyConvenience    75.45   3.68 147.22   2.530   29.824
> strategyQuality-strategyPrice          44.35 -27.42 116.12   1.487   29.824
>                                    p raw p Bonf p adj
> strategyPrice-strategyConvenience   0.301  0.904 0.553
> strategyQuality-strategyConvenience 0.014  0.043 0.037
> strategyQuality-strategyPrice       0.143  0.428 0.305
>
> This gives correct 95% confidence intervals and adjusted p-values for the
> Tukey multiple comparisons procedure.
>
>
>
> Next we issued
>
> summary(simtest(sales ~ strategy, type="Tukey", data=Xm15.01s))
>
> which produced the output
>
> Coefficients:
>                                    Estimate t value Std.Err. p raw p Bonf
> strategyQuality-strategyConvenience    75.45  -2.530   29.824 0.014  0.043
> strategyQuality-strategyPrice          44.35  -1.487   29.824 0.143  0.285
> strategyPrice-strategyConvenience      31.10  -1.043   29.824 0.301  0.301
>                                    p adj
> strategyQuality-strategyConvenience 0.037
> strategyQuality-strategyPrice       0.243
> strategyPrice-strategyConvenience   0.301
>
> Notice that the simtest output gives negative t statistics.  This is
> wrong because the point estimates are positive.
>
> The simtest Bonferroni p-values and adjusted p-values differ from the
> simint values by more than trivial amounts.
>
> We are puzzled that the two functions use different conventions on
> sequencing the contrasts.  For levels A,B,C, it looks like
>
> simint is using
> B-A
> C-A
> C-B
>
> and simtest is using
> C-A
> C-B
> B-A
>
>
> We verify all the p-values from the following code
>
>
> tt <- c(31.10,75.45,44.35)/29.824
> tt
>
>         2*(1-pt(tt, 57))            ## raw
> 3     * (2*(1-pt(tt, 57)))           ## Bonferroni
>        1-ptukey(tt*sqrt(2), 3, 57)  ## tukey
>
> ## It looks like simtest is using
> (3:1) * (2*(1-pt(tt[c(2,3,1)], 57))) ## simtest Bonferroni
> ## The subscript is there to account for the different sequencing.
> ## The (3:1) multiplier is strange.
>
> ## It looks like simtest is using approximately
> (1-ptukey(tt[c(2,3,1)]*sqrt(2), 3, 57)) / (1+c(0,1,3)*.12)^2
> ## We have no idea what that divisor is doing there other than
> ## approximating the answer that simtest is giving.
>
>
>
>
> ## Here is another example, this time using one of your datasets.
> ## Again, the p.Bonf and p.adj differ.  Again, the t.values for the
> ## simtest have the wrong sign.
> ## This is from
> ## c:/Program Files/R/R-2.2.1/library/multcomp/R-ex/Rex.zip/simtest.R
>
>> data(cholesterol)
>>
>> # adjusted p-values for all-pairwise comparisons in a one-way
>> # layout (tests for restricted combinations)
>> simtest(response ~ trt, data=cholesterol, type="Tukey", ttype="logical")
>
>        Simultaneous tests: Tukey contrasts
>
> Call:
> simtest.formula(formula = response ~ trt, data = cholesterol,
>    type = "Tukey", ttype = "logical")
>
> Contrast matrix:
>                      trt1time trt2times trt4times trtdrugD trtdrugE
> trt2times-trt1time  0       -1         1         0        0        0
> trt4times-trt1time  0       -1         0         1        0        0
> trtdrugD-trt1time   0       -1         0         0        1        0
> trtdrugE-trt1time   0       -1         0         0        0        1
> trt4times-trt2times 0        0        -1         1        0        0
> trtdrugD-trt2times  0        0        -1         0        1        0
> trtdrugE-trt2times  0        0        -1         0        0        1
> trtdrugD-trt4times  0        0         0        -1        1        0
> trtdrugE-trt4times  0        0         0        -1        0        1
> trtdrugE-trtdrugD   0        0         0         0       -1        1
>
> Adjusted P-Values
>
>                    p adj
> trtdrugE-trt1time   0.000
> trtdrugE-trt2times  0.000
> trtdrugD-trt1time   0.000
> trtdrugE-trt4times  0.000
> trt4times-trt1time  0.000
> trtdrugD-trt2times  0.000
> trtdrugE-trtdrugD   0.001
> trt2times-trt1time  0.042
> trt4times-trt2times 0.042
> trtdrugD-trt4times  0.044
>>
>>
>> summary(simtest(response ~ trt, data=cholesterol, type="Tukey", ttype="logical"))
>
>         Simultaneous tests: Tukey contrasts
>
> Call:
> simtest.formula(formula = response ~ trt, data = cholesterol,
>    type = "Tukey", ttype = "logical")
>
>         Tukey contrasts for factor trt
>
> Contrast matrix:
>                      trt1time trt2times trt4times trtdrugD trtdrugE
> trt2times-trt1time  0       -1         1         0        0        0
> trt4times-trt1time  0       -1         0         1        0        0
> trtdrugD-trt1time   0       -1         0         0        1        0
> trtdrugE-trt1time   0       -1         0         0        0        1
> trt4times-trt2times 0        0        -1         1        0        0
> trtdrugD-trt2times  0        0        -1         0        1        0
> trtdrugE-trt2times  0        0        -1         0        0        1
> trtdrugD-trt4times  0        0         0        -1        1        0
> trtdrugE-trt4times  0        0         0        -1        0        1
> trtdrugE-trtdrugD   0        0         0         0       -1        1
>
>
> Absolute Error Tolerance:  0.001
>
> Coefficients:
>                    Estimate t value Std.Err. p raw p Bonf p adj
> trtdrugE-trt1time     15.166 -10.507    1.443 0.000  0.000 0.000
> trtdrugE-trt2times    11.723  -8.122    1.443 0.000  0.000 0.000
> trtdrugD-trt1time      9.579  -6.637    1.443 0.000  0.000 0.000
> trtdrugE-trt4times     8.573  -5.939    1.443 0.000  0.000 0.000
> trt4times-trt1time     6.593  -4.568    1.443 0.000  0.000 0.000
> trtdrugD-trt2times     6.136  -4.251    1.443 0.000  0.000 0.000
> trtdrugE-trtdrugD      5.586  -3.870    1.443 0.000  0.001 0.001
> trt2times-trt1time     3.443  -2.385    1.443 0.021  0.043 0.042
> trt4times-trt2times    3.150  -2.182    1.443 0.034  0.043 0.042
> trtdrugD-trt4times     2.986  -2.069    1.443 0.044  0.044 0.044
>> summary(simint(response ~ trt, data=cholesterol, type="Tukey", ttype="logical"))
>
>        Simultaneous 95% confidence intervals: Tukey contrasts
>
> Call:
> simint.formula(formula = response ~ trt, data = cholesterol,
>    type = "Tukey", ttype = "logical")
>
>         Tukey contrasts for factor trt
>
> Contrast matrix:
>                      trt1time trt2times trt4times trtdrugD trtdrugE
> trt2times-trt1time  0       -1         1         0        0        0
> trt4times-trt1time  0       -1         0         1        0        0
> trtdrugD-trt1time   0       -1         0         0        1        0
> trtdrugE-trt1time   0       -1         0         0        0        1
> trt4times-trt2times 0        0        -1         1        0        0
> trtdrugD-trt2times  0        0        -1         0        1        0
> trtdrugE-trt2times  0        0        -1         0        0        1
> trtdrugD-trt4times  0        0         0        -1        1        0
> trtdrugE-trt4times  0        0         0        -1        0        1
> trtdrugE-trtdrugD   0        0         0         0       -1        1
>
> Absolute Error Tolerance:  0.001
>
> 95 % quantile:  2.842
>
> Coefficients:
>                    Estimate  2.5 % 97.5 % t value Std.Err. p raw p Bonf p adj
> trt2times-trt1time     3.443 -0.658  7.544   2.385    1.443 0.021  0.213 0.138
> trt4times-trt1time     6.593  2.491 10.694   4.568    1.443 0.000  0.000 0.000
> trtdrugD-trt1time      9.579  5.478 13.681   6.637    1.443 0.000  0.000 0.000
> trtdrugE-trt1time     15.166 11.064 19.267  10.507    1.443 0.000  0.000 0.000
> trt4times-trt2times    3.150 -0.952  7.251   2.182    1.443 0.034  0.344 0.205
> trtdrugD-trt2times     6.136  2.035 10.238   4.251    1.443 0.000  0.001 0.001
> trtdrugE-trt2times    11.723  7.621 15.824   8.122    1.443 0.000  0.000 0.000
> trtdrugD-trt4times     2.986 -1.115  7.088   2.069    1.443 0.044  0.443 0.251
> trtdrugE-trt4times     8.573  4.471 12.674   5.939    1.443 0.000  0.000 0.000
> trtdrugE-trtdrugD      5.586  1.485  9.688   3.870    1.443 0.000  0.003 0.003
>>
>
>
>
> --please do not edit the information below--
>
> Version:
> platform = i386-pc-mingw32
> arch = i386
> os = mingw32
> system = i386, mingw32
> status =
> major = 2
> minor = 2.1
> year = 2005
> month = 12
> day = 20
> svn rev = 36812
> language = R
>
> Windows XP Home Edition (build 2600) Service Pack 2.0
>
> Locale:
> LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
> States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
>
> Search Path:
> .GlobalEnv, package:multcomp, package:mvtnorm, package:abind, package:relimp,
> file:C:/PROGRA~1/R/R-22~1.1/library/Rcmdr/etc/HH/.RData, package:leaps, package:Rcmdr,
> package:car, package:tcltk, package:methods, package:stats, package:graphics, package:grDevices,
> package:utils, package:datasets, package:rcom, RcmdrEnv, Autoloads, package:base
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595



[Rd] bugs in simtest (PR#8670)
ripley@stats.ox.ac.uk