R-alpha: Re: ANOVA posthoc tests; Hotelling T^2, Bonferroni, Scheffe (fwd)

Gregory R. Warnes (warnes@biostat.washington.edu)
Tue, 21 Jan 1997 10:16:04 -0800 ()


Date: Tue, 21 Jan 1997 10:16:04 -0800 ()
From: "Gregory R. Warnes" <warnes@biostat.washington.edu>
To: r-testers <r-testers@stat.math.ethz.ch>
Subject: R-alpha: Re: ANOVA posthoc tests; Hotelling T^2, Bonferroni, Scheffe (fwd)


Attached is some code for several multiple comparison tests.  This was
originally posted to the S-news mailing list.

-------------------------------------------------------------------------------
    Gregory R. Warnes          | It is high time that the ideal of success
warnes@biostat.washington.edu  |  be replaced by the ideal of service.
                               |                       Albert Einstein
-------------------------------------------------------------------------------

---------- Forwarded message ----------
Date: Wed, 14 Aug 1996 10:33:34 -0700 (PDT)
From: John Wallace <jrw@fish.washington.edu>
To: "Grosof, David 314-362-2384" <GROSOF@am.seer.wustl.edu>
Cc: S-news <s-news@utstat.toronto.edu>
Subject: Re: ANOVA posthoc tests; Hotelling T^2, Bonferroni, Scheffe

On Tue, 13 Aug 1996, Grosof, David 314-362-2384 wrote:

> Hi,
>  In the S-news posting by Deb Carlson back on October 3, 1995, Deb Carlson 
> thanked you for providing code for Bonferroni, Hotelling T^2 and Scheffe 
> post-hoc tests. I contacted her today and her computer's hard disk drive died 
> and may not resuscitate. 
>  I would be grateful if you could send me the same.

---------------------------

Subject : Re: Hotelling T^2, Bonferroni, Scheffe 
On Wed, 7 Aug 1996, ""Jerry C. Shaw (914)945-3111"" wrote:

> Looking for a fairly extensive package in Splus. Anything
> around? I've seen the contributions from Wallace and others
> on the WEB at http://lib.stat.cmu.edu/s-news. Anything else?
>
---------------------------------

Here is what I have. I did not write any of the following functions, only 
put them together. Just the last function was not included in 
Brani's Dec. '94 summary.

-jrw

---------------------------------

Date: ~ Dec. '94

Summary of responses to "MULTIPLE COMPARISONS in Anova" question.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
>>I am looking for an S-function doing
>>multiple comparisons in Anova (Duncan,
>>Bonferroni, Scheffe,...).
>>Thanks, Brani
>>
>>Brani Vidakovic, ISDS, Duke University.

>>++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
>Please post a summary of any replies you receive.  Thanks.
>Bill Shipley
>bshipley@courrier.USherb.ca
>+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
>
>From:(Alan) ertlea@ohsu.edu
>I couldn't find a bonferroni correction and multiple comparisons
>for aov either so I wrote this little function. I
>left the comparison sets as vectors so x is a list of vectors you
>want to compare, and y is the starting
>alpha that will be corrected by the ( k over 2)
>factorial correction. It works for the comparisons that I have
>proof for (examples in books).
>Alan
>
>bonferroni_function(x, y)
>{
># x is a list of vectors
># y is the alpha 
>	num <- 0
>	denom <- 0
>	grandmean.num <- 0
>	sb2.num <- 0
>	nn <- numeric(length(x))
>	ss2 <- numeric(length(x))
>	for(i in 1:length(x)) {
>		nn[i] <- length(x[[i]])
>		ss2[i] <- (sum((x[[i]] - mean(x[[i]]))^2))/(nn[i] - 1
>			)
>		num <- ((nn[i] - 1) * ss2[i]) + num
>		denom <- nn[i] + denom
>		grandmean.num <- (nn[i] * mean(x[[i]])) + 
>			grandmean.num
>	}
>	sw2 <- num/(denom - length(x))
>	grandmean <- grandmean.num/denom
>	for(j in 1:length(x)) {
>		sb2.num <- (nn[j] * ((mean(x[[j]]) - grandmean)^2)) + 
>			sb2.num
>	}
>	sb2 <- sb2.num/(length(x) - 1)
>	FF <- sb2/sw2
>	p.val <- 1 - pf(FF, length(x) - 1, denom - length(x))
>	effort <- numeric(length(x) * (length(x) - 1))
>	effort1 <- character(length(x) * (length(x) - 1))
>	tt <- matrix(effort, ncol = length(x))
>	pp <- matrix(effort, ncol = length(x))
>	qq <- matrix(effort1, ncol = length(x))
>	alpha.star <- y/((length(x) * (length(x) - 1))/2)
>	deg.free <- denom - length(x)
>	for(k in 1:(length(x) - 1)) {
>		for(l in 2:length(x)) {
>			if(k >= l) {
>				tt[k, l] <- 0
>			}
>			else {
>				tt[k, l] <- (mean(x[[k]]) - mean(x[[l
>				  ]]))/sqrt(sw2 * ((1/nn[k]) + (1/nn[
>				  l])))
>				pp[k, l] <- 2 * (pt(tt[k, l], 
>				  deg.free))
>				pp[k, k] <- 0
>				if(pp[k, l] == 0) {
>				  qq[k, l] <- c("Zero")
>				}
>				else if(pp[k, l] > alpha.star) {
>				  qq[k, l] <- c("Not")
>				}
>				else if(pp[k, l] < alpha.star) {
>				  qq[k, l] <- c("Signif")
>				}
>			}
>		}
>	}
>	results <- list(grandmean, sw2, sb2, FF, p.val, tt, 
>		alpha.star, pp, qq)
>	results <- list(Grand.mean = grandmean, Within.var = sw2, 
>		Between.var = sb2, F.value = FF, p.value = p.val, 
>		bonferroni.t = tt, alpha.star = alpha.star, 
>		bonferroni.p = pp, bonferroni.sig = qq)
>	results
>}
>++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
>
>From: Luciano Molinari <molinari@stat.math.ethz.ch>
>I am very much interested in functions for multiple comparisons, especially
>of the Scheffe type, which make use of the object-oriented approach
>underlying the modeling environment of S-Plus.
>Please let me know what the reactions to your query were.
>I intend to write myself some simple functions to this purpose, but, as
>usual, I could not find the time so far. Below I include two very special
>functions. The second implements a Bonferroni-type procedure of Holm
>(S. Holm, Scand. J. Statist.,6:65-70,1979); the first is for the special
>case of simultaneous testing equality of the components of two vectors of
>means, as described in Morrison.
>++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
>"twohotel"<-
>function(x1, x2, alfa = 0.01)
>{
>#Hotelling T^2 test for two groups, after Morrison, p.137
>#The FF calculated here appears to be the same at the approx F
>#calculated by S-Plus Manova
>	n1 <- dim(x1)[1]
>	n2 <- dim(x2)[1]
>	pp <- dim(x1)[2]
>	dfden <- (n1 + n2 - pp - 1)
>	c1 <- (n1 * n2)/(n1 + n2)
>	c2 <- dfden/(n1 + n2 - 2)/pp
>	Falfa <- qf(1 - alfa, pp, dfden)
>	T2alfa <- Falfa/c2
>	print(T2alfa)
>	x1bar <- t(x1) %*% rep(1, dim(x1)[1])/dim(x1)[1]
>	x2bar <- t(x2) %*% rep(1, dim(x2)[1])/dim(x2)[1]
>	maha <- twomaha(x1, x2)
>	lll <- sqrt(diag(maha$S)/c1 * T2alfa)
>	T2 <- c1 * maha$Maha
>	FF <- c2 * T2
>	list(Simul = cbind(X1Bar = x1bar, X2Bar = x2bar, Delta = lll), 
>		Hotel = c(Alfa = alfa, T2 = T2, F = FF, dfnum = pp, 
>		dfden = dfden, P = 1 - pf(FF, pp, dfden)))
>}
>"twomaha"<-
>function(x1, x2, aver = mean, calcov = cov.wt)
>{
>#Mahalanobi's distance for two groups, after Morrison, p.235
>#calcov should yield the MLE of the covariance matrix
>	n1 <- dim(x1)[1]
>	n2 <- dim(x2)[1]
>	S <- (calcov(x1)$cov * n1 + calcov(x2)$cov * n2)/(n1 + n2 - 2)
>	dx <- apply(x1, 2, aver) - apply(x2, 2, aver)
>	list(Maha = sum(dx * solve(S, dx)), S = S, Delta = dx)
>}
>"holm.srB"<-
>function(pvals, w = 1, alfa, hypnam = names(pvals))
>{
>	if(!is.null(hypnam))
>		hypnam <- hypnam[!is.na(pvals)]
>	w <- rep(w, len = length(pvals))
>	w <- w[!is.na(pvals)]
>	pvals <- pvals[!is.na(pvals)]
>	n <- length(pvals)
>	if(is.null(hypnam))
>		hypnam <- 1:n
>	ooo <- order(pvals/w)
>	w <- w[ooo]
>	hypnam <- hypnam[ooo]
>	pvals <- pvals[ooo]
>	names(pvals) <- hypnam
>	qs <- pvals/w * (sum(w) - c(0, cumsum(w)[1:(n - 1)]))
>	if(!missing(alfa)) {
>		maxr <- max((1:n)[cumsum(qs <= alfa) == (1:n)])
>		if(!is.na(maxr))
>			pvals[1:maxr]
>		else NA
>	}
>	else {
>		maxr <- (1:n)[cummax(qs) == qs]
>		list(rejs = maxr, Pmax = qs[maxr])
>	}
>}
>+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# I recently wrote the following Scheffe mc function. It is
a first draft and it can be beautified.        

> Scheffe <- function(data, contrast, alpha = 0.05, r = 3)
{
        means <- apply(data, 2, mean.na)
        cat("\n means:", means, "\n")
        fit <- matrix(means, nrow = dim(data)[1], ncol = dim(data)[2], byrow =
                T, dimnames = dimnames(data))
        fit <- fit + 0 * data
        res <- data - fit
        ngroup <- apply(!is.na(data), 2, sum)
        if((contrast %*% ngroup)[1, 1] != 0)
                stop("The contrast is", " untestable for this design.", (
                        contrast %*% ngroup)[1, 1])
        cat("\n ngroup:", ngroup, "\n")
        sse <- sum.na(res^2)
        dfe <- (sum(!is.na(data)) - dim(data)[2])
        mse <- sse/dfe
        C <- (contrast %*% means)[1, 1] #dim(C) <- NULL
        cat("\n C=", round(C, r), "\n")
        SC <- sqrt(mse * sum(contrast^2/ngroup))
        Salp <- SC * sqrt((dim(data)[2] - 1) * qf(1 - alpha, dim(data)[2] - 1,
                dfe))
        cat("\n C=", round(C, r), " SC=", round(SC, r), " Salp=", round(Salp, r
                ), "\n")
        cat("\n", 1 - alpha, "*100% confidence interval is:(", round(C - Salp,
                r), ",", round(C + Salp, r), ")\n")
        if(abs(C) < Salp) {
                   cat("\n The contrast is tested to be 0 at the level", alpha,
                        "\n")
        }
        else {
                cat("\n The contrast is tested not to be 0 at the level", alpha,
                        "\n")
        }
}
#---------------------------------------------------------------
#----------------------------------------------------------------
I tested the function on Montgomery's "Design and Analysis of
Experiments", 2nd Edition  example, page 63-64. The notation 
(C, SC, Salp) is the same. If a design is unbalanced, the empty
slots in the data matrix should be filled with NA and  care
should be taken about contrasts.
To be testable contrasts have to satisfy Sum(n_i * c_i)=0, where
n_i are the group sizes and c_i the corresponding contrast components.

    
> data
     [,1] [,2] [,3] [,4] [,5]
[1,]    7   12   14   19    7
[2,]    7   17   18   25   10
[3,]   15   12   18   22   11
[4,]   11   18   19   19   15
[5,]    9   18   19   23   11
> contrast
[1]  1  0  1 -1 -1
> Scheffe(data, contrast, alpha=0.01)

 means: 9.8 15.4 17.6 21.6 10.8

 ngroup: 5 5 5 5 5

 C= -5

 C= -5  SC= 2.539  Salp= 10.69

 0.99 *100% confidence interval is:( -15.69 , 5.69 )
#----------------------------------------------------------------
Brani

--------------------------------------------------------------------------

>From fms@chemelex.comWed Aug 14 10:04:48 1996
Date: Mon, 14 Mar 94 11:31:25 -0800
From: Fred Schenkelberg <fms@chemelex.com>
To: "John R. Wallace" <jrw@fish.washington.edu>
Subject: How to do mulitiple comparisons?


Just this one response,

Begin forwarded message:

Date: Thu, 10 Mar 1994 00:03:09 -0600 (CST)
From: Ruben Smith <rasmith@umaxc.weeg.uiowa.edu>
Subject: Re: How to do mulitiple comparisons?
To: Fred Schenkelberg <fms@chemelex.com>
In-Reply-To: <9403092159.AA00519@chemelex>


I made this function. I recommed  use length level of 1.

> plot.comparisons
function(model, f1, f2 = NULL, f3 = NULL, tplot = 1, alpha  
= 0.05, 

	main = NULL, xlab = "RESPONSE", ylab = paste(f1,  
f2, f3, 

	sep = "."))
{
# DESCRIPTION 

# Produce a Least Significant Interval Plot or Bonferroni  
Significant        

# Difference Plot for main effect, two or three factor  
interaction.    

# ARGUMENT 

#          Model: aov model.
#       f1,f2,f3: factor names. (Introduce the factor names  
in double
#                 quotes in the order that these are in the  
aov model)
#          tplot: type plot (1,2 for LSI, BSI)
#          alpha: Significancy level.
# Elaborated by Ruben Smith.

	v1 <- vector("numeric", 2)
	v2 <- vector("numeric", 2)
	mse <- sum(residuals(model)^2)/model$df.residual
	m <- model.frame(model)
	repl <- replications(model, m)
	if(is.null(f3)) {
		if(is.null(f2)) {
			inte <- m[, f1]
			aux <- f1
		}
		else {
			inte <- interaction(m[, f1], m[,  
f2])
			aux <- paste(f1, "ean)
	n <- repl[aux]
	if(tplot == 1) {
		if(is.null(main))
			main <- "LEAST SIGNIFICANT INTERVAL  
PLOT"
		d <- 1
	}
	else {
		if(is.null(main))
			main <- "BONFERRONI SIGNIFICANT  
DIFFERENCE 

PLOT"
		d <- (a * (a - 1))/2
	}
	lsd <- qt(1 - alpha/(2 * d), mod seq(1, a), lev)
	for(i in 1:a) {
		v2[1] <- factor[i]
		v2[2] <- factor[i]
		v1[1] <- av[factor[i]] + lsi
		v1[2] <- av[factor[i]] - lsi
		lines(v1, v2, type = "l")
	}
	points(av, factor, type = "p")
	box()
}


------------------------------------------------
Ruben Smith         rasmith@umaxc.weeg.uiowa.edu    

702 E Washington St. Apto 5                          

Iowa City, IA 52240        Tel: (319)338-1503       

USA                                                  

                                                      

------------------------------------------------       




On Wed, 9 Mar 1994, Fred Schenkelberg wrote:

> Hi, after doing an anova on three treaments and find a  

> signifcant difference, I'd like to do multiple  
comparisons  

> using Tukey's W and/or Scheffe's procedures.
> 

> Are these procedures built into Splus?
> 

> Thanks,  Fred
> 

> 

> Fred Schenkelberg
> NeXTmail enjoyed
> fms@chemelex.com
> Redwood City, CA





-- 
John Wallace
University of Washington                    ^    ^    ^
Fisheries Research Institute               / \  / \  / \   ^
Box 357980                                 / \  / \   |   / \
Seattle, WA 98195-7980                      |    |  o__~  / \
PHONE   (206) 543-1513                  @ @         /\/\   |
FAX     (206) 685-7471                   ~    
E-MAIL        jw@u.washington.edu              o
WWW        http://pisces.fish.washington.edu/people/jrw/Wallace.html
                                                o  _///_ //
                                                <`)=  _<<
                                                    \\\  \\





=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
r-testers mailing list -- For info or help, send "info" or "help",
To [un]subscribe, send "[un]subscribe"
(in the "body", not the subject !)  To: r-testers-request@stat.math.ethz.ch
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-