bugs in simtest (PR#8670)

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

bugs in simtest (PR#8670)

Richard M. Heiberger
# 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
#
#       [hidden email]
#
######################################################

This report is joint from Richard Heiberger <[hidden email]>
and Burt Holland <[hidden email]>.
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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: bugs in simtest (PR#8670)

Brian Ripley
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, [hidden email] 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
> #
> #       [hidden email]
> #
> ######################################################
>
> This report is joint from Richard Heiberger <[hidden email]>
> and Burt Holland <[hidden email]>.
> 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
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
>

--
Brian D. Ripley,                  [hidden email]
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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel