Re: [R] Hmisc / Design question

From: Frank E Harrell Jr <f.harrell_at_vanderbilt.edu>
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.html
Received 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