# 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
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