[R] Hmisc / Design question

[Frank E Harrell Jr]

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