[Rd] bugs in simtest (PR#8670)
ripley@stats.ox.ac.uk
ripley at stats.ox.ac.uk
Thu Mar 9 08:53:45 CET 2006
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
More information about the R-devel
mailing list