[Rd] bugs in simtest (PR#8670)

rmh@temple.edu rmh at temple.edu
Thu Mar 9 02:48:44 CET 2006


# 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



More information about the R-devel mailing list