Date:         Tue, 21 Dec 1999 15:58:41 GMT
Reply-To:     olmsted_a@MY-DEJA.COM
Sender:       "SAS(r) Discussion" <SAS-L@LISTSERV.UGA.EDU>
From:         olmsted_a@MY-DEJA.COM
Organization: Deja.com - Before you buy.
Subject:      Re: exact CMH (Cochran-Mantel-Haenszel) test for 2x2xk tables

Proc multtest to the rescue! The code does exactly what I need.
Below, for those who have Agresti's 1990 Categorical data analysis
book, is some code reproducing the Table 7.9 example on p. 233
(both the two-tailed and the one-tailed test):

data WORK (keep = dose delay cured died) ;
   infile cards stopover ;
   attrib dose     label='Penicillin level (MUs)' ;
   attrib delay    label='Injection delay' length=$4 ;
   attrib cured    label='No. of rabbits cured' length=3 ;
   attrib died     label='No. of rabbits who died' length=3 ;
   input dose delay & cured died ;
 /*
@ Agresti, Alan (1990). Categorical data analysis. New York: John Wiley
& Sons, Inc.
@ Table 7.9, p. 233
@ Example of Cochran-Mantel-Haenszel test
@ dose     delay  cured   died
 */
cards ;
   0.125  None         0       6
   0.125  1½ h         0       5
    0.25  None         3       3
    0.25  1½ h         0       6
     0.5  None         6       0
     0.5  1½ h         2       4
       1  None         5       1
       1  1½ h         6       0
       4  None         2       0
       4  1½ h         5       0
;

proc print data=WORK ;
   id dose delay ;
   title3 'WORK for checking' ;
run ;

proc transpose data=WORK out=WORK (rename=(col1=count)) name=response
label=label ;
   by dose delay notsorted ;
   var cured died ;
   attrib response label='Response (cured|died)' ;
run ;

proc print data=WORK noobs ;
   by dose delay notsorted ;
   id dose delay ;
   var response label count ;
   title3 'WORK for checking' ;
run ;

* make None sort before 1 1/2 h ;
proc sql feedback ;
   update work set delay = '0 h' where delay='None' ;

proc freq data=WORK formchar(1,2,7)='   ' ;
   tables dose*delay*response / cmh nopercent norow nocol ;
   weight count ;
   title3 'Agresti gives the large-sample CMH p-value as .047' ;
run ;

* response (column) variable must be numeric ;
proc sql feedback ;
   create table work as
   select dose, delay, response, count, response = 'CURED' as cured
label='Cured?' format=1.
   from work ;

proc freq data=WORK formchar(1,2,7)='   ' ;
   tables response * cured / missing list nopercent ;
   title3 'check recoding of response variable' ;
run ;

proc multtest data=work ;
   class delay ;
   strata dose ;
   freq count ;
   test ca(cured/permutation=100) ;
   contrast "1 vs 2" 1 -1 ; * not needed ;
   title3 'Agresti gives the exact 2-sided p-value as .040' ;
run ;
quit ;

proc multtest data=work ;
   class delay ;
   strata dose ;
   freq count ;
   test ca(cured/permutation=100 uppertailed) ; * alternative hyp is
immediate injection better than delayed ;
   contrast "1 vs 2" 1 -1 ; * not needed ;
   title3 'Agresti gives the exact 1-sided p-value as .020' ;
run ;
quit ;




Sent via Deja.com http://www.deja.com/
Before you buy.