From: Frank E Harrell Jr <f.harrell_at_vanderbilt.edu>

Date: Wed 03 Aug 2005 - 06:16:50 EST

Date: Wed 03 Aug 2005 - 06:16:50 EST

Greg Tarpinian wrote:

*> All,
**>
*

> I have been reading Dr. Harrell's excellent

*> "Regression Modeling Strategies" book and trying out
**> the exercises. I understand that contrast( ) is used
**> to obtain contrasts between two variables for given
**> levels of other nuisance variables; is there a way
**> to use contrast( ) to obtain, for example, Scheffe
**> confidence intervals / hypothesis tests for many
**> post hoc contrasts at once?
**>
**> Any help would be much appreciated,
**>
**> Greg
*

Steven Novick (cc'd) has written some nice code that I regret I haven't had time to add into contrast.Design for doing simulation-based adjustment for multiple comparisons. Below is the code he sent me.

## Program: multiple.adjustment.R ## Version: 1 ## Author: Steven Novick ## Date: May 21, 2004 ## Purpose: Compute the exact T-critical value for multiple comparisons based on the paper ## Don Edwards and Jack Berry (1987). "The Efficiency of Simulation-Based Multiple Comparison". ## Biometrics 43, pp913-928.multiple.adjustment = function(fit, cont1, cont2, alpha=.05, m=79999) {

## fit = object of class "ols" from library Design ## cont1, cont2 = data.frame with elements to be contrasted. See function "contrast" in Design. ## alpha = test significance level ## m = number of monte-carlo runs V=fit$var/fit$stats["Sigma"]^2 ## Var( coef(fit) ) = sigma^2 * V C=contrast(fit, cont1, cont2)$X ## Contrast matrix df = fit$df ## Error degrees of freedom r = (m+1)*(1-alpha) # P( W < w[r] ) = alpha if ( floor(r) != ceiling(r) ) stop("(m+1)*(1-alpha) must be an integer.") ## Create random numbers n.contrast=nrow(C) n.coef=length(coef(fit)) G=chol(V) U=apply(C, 1, function(ctr){( G%*%ctr )/as.vector(sqrt(t(ctr)%*%V%*%ctr)) }) W = sort(sapply(1:m, function(i){ z=rnorm(n.coef); y=sqrt(rchisq(1,df)/df); w=abs(t(U)%*%z)/y; return(max(w)) }))

w.alpha = W[r] ## Cut-off value for testing; 95% CI = beta.hat + c(-1, 1)*w.alpha* SE(beta.hat)

alpha.star = 2*(1-pt(w.alpha, df)) ## What is alpha if using T-statistic ?

return(list(w.alpha=w.alpha, alpha.star=alpha.star)) }

-- Frank E Harrell Jr Professor and Chair School of Medicine Department of Biostatistics Vanderbilt University ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.htmlReceived on Wed Aug 03 06:21:06 2005

*
This archive was generated by hypermail 2.1.8
: Sun 23 Oct 2005 - 15:01:56 EST
*