Multiple comparisons

JIMRC & package multcomp

For the starting point see the  multiple comparison page.  We follow the same steps until we have the list "fm1" by the command:

  fm1= aov(Rust ~ Brand)

Remark: the TukeyHSD( ) function requires a list made by the aov( ) function.

Multiple comparisons by TukeyHSD:  see previous multiple comparison page.

Multiple comparisons by three major methods as described in NKNW 96 for one-way ANOVA.

1. Bonferroni (a specified number of comparisons)
2. Scheffé (all contrasts)
3. Tukey (all pairswise comparisons, same method as TukeyHSD)

1. all.pairs(r)     which is used indirectly by next functions
2. tukeyCI(fitted, nis, df, MSE, conf.level=.95)
3. scheffeCI(fitted, nis, df, MSE, conf.level=.95)
4. bonferroniCI(fitted, nis, df, MSE, conf.level=.95)

Copy/paste the entire text below next horizontal line and before the subsequent one in R  

(remark: R-code is normally provided in a file: MultiplecompJIMRC.R ; 

save this file somewhere and read it into R by "File/Source R code" )

 all.pairs <- function(r)
  list(first = rep(1:r,rep(r,r))[lower.tri(diag(r))],
       second = rep(1:r, r)[lower.tri(diag(r))])

tukeyCI <- function(fitted, nis, df, MSE, conf.level=.95){
  ## fitted is a sequence of means
  ## nis is a corresponding sequence of sample sizes for each mean
  ## df is the residual df from the ANOVA table
  ## MSE = mean squared error from the ANOVA table
  ## conf.level is the family-wise confidence level, defaults to .95
  r <- length(fitted)
  pairs <- all.pairs(r)
  diffs <- fitted[pairs$first] - fitted[pairs$second]
  df <- sum(nis) - r
  T <- qtukey(conf.level, r, df)/sqrt(2)
  hwidths <-  T*sqrt(MSE*(1/nis[pairs$first] + 1/nis[pairs$second]))
  val <- cbind(diffs - hwidths, diffs, diffs + hwidths)
  dimnames(val) <- list(paste("mu",pairs$first," - mu", pairs$second,
  sep=""), c("Lower", "Diff","Upper"))
  val
}

scheffeCI <- function(fitted, nis, df, MSE, conf.level=.95){
  ## fitted is a sequence of means
  ## nis is a corresponding sequence of sample sizes for each mean
  ## df is the residual df from the ANOVA table
  ## MSE = mean squared error from the ANOVA table
  ## conf.level is the family-wise confidence level, defaults to .95
  r <- length(fitted)
  pairs <- all.pairs(r)
  diffs <- fitted[pairs$first] - fitted[pairs$second]
  T <- sqrt((r-1)*qf(conf.level,r-1,df))
  hwidths <-  T*sqrt(MSE*(1/nis[pairs$first] + 1/nis[pairs$second]))
  val <- cbind(diffs - hwidths, diffs, diffs + hwidths)
  dimnames(val) <- list(paste("mu",pairs$first," - mu", pairs$second,
  sep=""), c("Lower", "Diff","Upper"))
  val
}

bonferroniCI <- function(fitted, nis, df, MSE, conf.level=.95){
  ## fitted is a sequence of means
  ## nis is a corresponding sequence of sample sizes for each mean
  ## df is the residual df from the ANOVA table
  ## MSE = mean squared error from the ANOVA table
  ## conf.level is the family-wise confidence level, defaults to .95
  r <- length(fitted)
  pairs <- all.pairs(r)
  diffs <- fitted[pairs$first] - fitted[pairs$second]
  T <- qt(1-(1-conf.level)/(2*r*(r-1)),df)
  hwidths <-  T*sqrt(MSE*(1/nis[pairs$first] + 1/nis[pairs$second]))
  val <- cbind(diffs - hwidths, diffs, diffs + hwidths)
  dimnames(val) <- list(paste("mu",pairs$first," - mu", pairs$second,
  sep=""), c("Lower", "Diff","Upper"))
  val
}

We have now 4 new user defined functions. If we type in 

bonferroniCI

without argumentlist we get the R-code of the function.

We prepare the agrument list for the functions: 

• averages,

• the number of replications, 

• dfMSE and 

• MSE 

for Rust as function of each factor-level Brand.

> Rust.means <- tapply(Rust,Brand,mean);Rust.means
     1      2      3      4 
 43.14  89.44  67.95  40.47 

> Rust.len <- tapply(Rust,Brand,length);Rust.len
       1   2   3   4 
      10  10  10  10 

> dfMSE=fm1$df.residual;dfMSE
[1] 36

> MSE=sum(fm1$residuals^2)/dfMSE;MSE
[1] 6.139833

Check if these last two values are OK by executing anova(fm1)
Now we are ready for user defined functions:

> tukeyCI(Rust.means, Rust.len, dfMSE, MSE, conf=.95)
                   Lower         Diff         Upper
mu1 - mu2    -49.2844635       -46.30    -43.315536
mu1 - mu3    -27.7944635       -24.81    -21.825536
mu1 - mu4     -0.3144635         2.67      5.654464
mu2 - mu3     18.5055365        21.49     24.474464
mu2 - mu4     45.9855365        48.97     51.954464
mu3 - mu4     24.4955365        27.48     30.464464

Compare this result with the base function:

fm1Tukey=TukeyHSD(fm1,"Brand") ; fm1Tukey

and with the other methods:

bonferroniCI(Rust.means, Rust.len, dfMSE, MSE, conf=.95)

scheffeCI(Rust.means, Rust.len, dfMSE, MSE, conf=.95)

Remark: Other methods for multiple comparisons and contrasts.

 (VERY ADVANCED... niet voor Proeftechniek !!!!)

Unfortunately R provides only the Tukey-procedure in the base-package. Although Tukey is the most common and important method; it more specifically handles all pairwise comparisons simultaneously.

Other methods are via a contributed package: multcomp.  When we check the requirements we see:  "Depends R (>= 1.5.1), mvtnorm (>=0.5-4)" ; so it is meant for version 1.51 and upwards and we need a supporting package: mvtnorm. So we need to install two packages.

Packages in binary (compiled ) form  can be saved as *.zip files.

Surf starting from the R binaries to the appropiate (depending on the operating system) subdirectory: 

Index of /CRAN/bin/windows/contrib

 multcomp.zip   

 mvtnorm.zip 

Right click and save.

Then install the package by "Packages/install packages from local zip drive"

Activate the two packages by "Packages/load package" and take mvtnorm and multcomp from the list.

31 March, 2003 by Guido Wyseure