# 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