Categorical Data Analysis: Chapter 4
data arth;
input treat $ response $ count @@;
datalines;
active none 13 active some 7 active marked 21
placebo none 29 placebo some 7 placebo marked 7
;
proc freq data=arth order=data;
weight count;
tables treat*response / chisq nocol nopct;
run;
data arth;
input gender $ treat $ response $ count @@;
datalines;
female active none 6 female active some 5 female active marked 16
female placebo none 19 female placebo some 7 female placebo marked 6
male active none 7 male active some 2 male active marked 5
male placebo none 10 male placebo some 0 male placebo marked 1
;
proc freq data=arth order=data;
weight count;
tables gender*treat*response / cmh nocol nopct;
run;
proc freq data=arth order=data;
weight count;
tables gender*treat*response/cmh scores=modridit nocol nopct;
run;
data arth2;
input gender $ treat $ response $ count @@;
datalines;
female placebo none 19 female placebo some 7 female placebo marked 6
female active none 6 female active some 5 female active marked 16
male placebo none 10 male placebo some 0 male placebo marked 1
male active none 7 male active some 2 male active marked 5
;
proc freq data=arth2 order=data;
weight count;
tables gender*treat*response / measures scores=modridit;
output out=somerout smdcr;
ods output CrosstabFreqs=myFreqs;
run;
data ns;
set myFreqs; where _type_="110";
run;
proc iml;
use somerout var{_SMDCR_ E_SMDCR};
read all into combined;
use ns var{Frequency};
read all into total;
print total;
WH=(total[1,]#total[2,]) / (total[1,]+total[2,]) // (total[3,]#total[4,]) / (total[3,]+total[4,]);
print wh;
GH=(combined[,1]+1)/2;
SEGH=(combined[,2])/2;
print wh GH SEGH;
QH= (GH[1] - GH[2])**2 / SSQ(SEGH);
pvalue=1-probchi(QH,1);
print GH SEGH QH pvalue;
gbar1=sum(GH / SEGH##2);
gbar2= sum (1/SEGH##2);
gtilde1=sum(GH # WH);
gtilde2=sum(WH);
gbar=gbar1/gbar2;
gtilde=gtilde1/gtilde2;
gtildevar=(sum(WH##2 # SEGH##2)) /((sum(WH))**2);
segtilde=sqrt(gtildevar);
qgtilde=(gtilde-.5)**2/gtildevar;
gbarvar=inv (sum(1/SEGH##2));
qgbar=(gbar-.5)**2/gbarvar;
qgtilde=(gtilde-.5)**2/gtildevar;
print gbar gbarvar qgbar;
print gtilde gtildevar qgtilde;
quit;
data colds;
input gender $ residence $ per_cold count @@;
datalines;
female urban 0 45 female urban 1 64 female urban 2 71
female rural 0 80 female rural 1 104 female rural 2 116
male urban 0 84 male urban 1 124 male urban 2 82
male rural 0 106 male rural 1 117 male rural 2 87
;
proc freq data=colds order=data;
weight count;
tables gender*residence*per_cold / all nocol nopct;
run;
proc glm;
class gender residence;
freq count;
model per_cold = gender residence ;
estimate 'd' residence -1 1;
run;
data tobacco;
length risk $11. ;
input f_usage $ risk $ usage $ count @@;
datalines;
no minimal no 59 no minimal yes 25
no moderate no 169 no moderate yes 29
no substantial no 196 no substantial yes 9
yes minimal no 11 yes minimal yes 8
yes moderate no 33 yes moderate yes 11
yes substantial no 22 yes substantial yes 2
;
proc freq;
weight count;
tables f_usage*risk*usage /cmh chisq measures trend;
tables f_usage*risk*usage /cmh scores=modridit;
run;
proc freq;
weight count;
tables f_usage*risk*usage /cmh;
tables f_usage*risk*usage /cmh scores=modridit;
run;
data pain;
input diagnosis $ treatment $ response $ count @@;
datalines;
I placebo no 26 I placebo yes 6
I dosage1 no 26 I dosage1 yes 7
I dosage2 no 23 I dosage2 yes 9
I dosage3 no 18 I dosage3 yes 14
I dosage4 no 9 I dosage4 yes 23
II placebo no 26 II placebo yes 6
II dosage1 no 12 II dosage1 yes 20
II dosage2 no 13 II dosage2 yes 20
II dosage3 no 1 II dosage3 yes 31
II dosage4 no 1 II dosage4 yes 31
;
proc freq order=data;
weight count;
tables treatment*response / chisq;
tables diagnosis*treatment*response / chisq cmh;
tables diagnosis*treatment*response / scores=modridit cmh;
run;
proc freq order=data;
weight count;
tables diagnosis*response*treatment / cmh;
tables diagnosis*treatment*response / cmh;
run;
data mice;
input LogRank Treatment $ count @@;
datalines;
0.909 Ca 1 0.909 Ce 1
0.709 Ca 3 0.709 Ce 1
0.334 Ca 5 0.334 Ce 1
0.234 Ca 1 0.234 Ce 0
-0.099 Ca 1 -0.099 Ce 2
-0.433 Ca 0 -0.433 Ce 2
-0.933 Ca 1 -0.933 Ce 1
-1.433 Ca 0 -1.433 Ce 1
-2.433 Ca 0 -2.433 Ce 1
;
proc sort data=mice;
by LogRank;
run;
proc freq;
weight count;
tables LogRank*Treatment / norow nocol nopct scorout chisq;
tables LogRank*Treatment / noprint scores=rank scorout chisq;
exact mhchi;
run;