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