[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