Categorical Data Analysis: Chapter 13
options nodate nonumber ps=200 ls=80 formdlim=' ';
data colds;
input sex $ residence $ periods count @@;
datalines;
female rural 0 45 female rural 1 64 female rural 2 71
female urban 0 80 female urban 1 104 female urban 2 116
male rural 0 84 male rural 1 124 male rural 2 82
male urban 0 106 male urban 1 117 male urban 2 87
;
run;
proc catmod;
weight count;
response means;
model periods = sex residence sex*residence /freq prob;
run;
proc catmod;
weight count;
response means;
model periods = sex residence;
run;
proc catmod;
population sex residence;
weight count;
response means;
model periods = sex;
run;
data cpain;
input dstatus $ invest $ treat $ status $ count @@;
datalines;
I A active poor 3 I A active fair 2 I A active moderate 2
I A active good 1 I A active excel 0
I A placebo poor 7 I A placebo fair 0 I A placebo moderate 1
I A placebo good 1 I A placebo excel 1
I B active poor 1 I B active fair 6 I B active moderate 1
I B active good 5 I B active excel 3
I B placebo poor 5 I B placebo fair 4 I B placebo moderate 2
I B placebo good 3 I B placebo excel 3
II A active poor 1 II A active fair 0 II A active moderate 1
II A active good 2 II A active excel 2
II A placebo poor 1 II A placebo fair 1 II A placebo moderate 0
II A placebo good 1 II A placebo excel 1
II B active poor 0 II B active fair 1 II B active moderate 1
II B active good 1 II B active excel 6
II B placebo poor 3 II B placebo fair 1 II B placebo moderate 1
II B placebo good 5 II B placebo excel 0
III A active poor 2 III A active fair 0 III A active moderate 3
III A active good 3 III A active excel 2
III A placebo poor 5 III A placebo fair 0 III A placebo moderate 0
III A placebo good 8 III A placebo excel 1
III B active poor 2 III B active fair 4 III B active moderate 1
III B active good 10 III B active excel 3
III B placebo poor 2 III B placebo fair 5 III B placebo moderate 1
III B placebo good 4 III B placebo excel 2
IV A active poor 8 IV A active fair 1 IV A active moderate 3
IV A active good 4 IV A active excel 0
IV A placebo poor 5 IV A placebo fair 0 IV A placebo moderate 3
IV A placebo good 3 IV A placebo excel 0
IV B active poor 1 IV B active fair 5 IV B active moderate 2
IV B active good 3 IV B active excel 1
IV B placebo poor 3 IV B placebo fair 4 IV B placebo moderate 3
IV B placebo good 4 IV B placebo excel 2
;
proc catmod order=data;
weight count;
response 1 2 3 4 5;
model status=dstatus|invest|treat;
run;
proc catmod order=data;
weight count;
response 1 2 3 4 5;
model status=dstatus|invest|treat@2;
run;
proc catmod order=data;
weight count;
response 1 2 3 4 5;
model status=dstatus invest treat;
run;
proc catmod order=data;
weight count;
response 1 2 3 4 5;
model status=dstatus invest treat;
contrast 'I versus II' dstatus 1 -1 0;
contrast 'I versus III' dstatus 1 0 -1;
contrast 'I versus IV' dstatus 2 1 1;
contrast 'II versus III' dstatus 0 1 -1;
contrast 'II versus IV' dstatus 1 2 1;
contrast 'III versus IV' dstatus 1 1 2;
contrast 'dstatus' dstatus 1 0 0 ,
dstatus 0 1 0 ,
dstatus 0 0 1 ;
run;
data byss;
input workplace $ em_years $ smoking $ status $ count @@;
datalines ;
dusty <10 yes yes 30 dusty <10 yes no 203
dusty <10 no yes 7 dusty <10 no no 119
dusty >=10 yes yes 57 dusty >=10 yes no 161
dusty >=10 no yes 11 dusty >=10 no no 81
notdusty <10 yes yes 14 notdusty <10 yes no 1340
notdusty <10 no yes 12 notdusty <10 no no 1004
notdusty >=10 yes yes 24 notdusty >=10 yes no 1360
notdusty >=10 no yes 10 notdusty >=10 no no 986
;
run;
proc catmod order=data;
weight count;
response marginals;
model status = workplace|em_years|smoking;
run;
proc catmod order=data;
weight count;
response marginals;
model status = workplace|em_years workplace|smoking
em_years|smoking;
run;
proc catmod order=data;
weight count ;
response marginals;
model status = workplace|em_years workplace|smoking;
run;
proc catmod order=data;
weight count;
response marginals;
model status = workplace
em_years(workplace='dusty')
em_years(workplace='notdusty')
smoking(workplace='dusty')
smoking(workplace='notdusty') ;
run;
proc catmod order=data;
weight count ;
response marginals;
model status = workplace
em_years(workplace='dusty')
smoking(workplace='dusty') / pred;
run;
data pain (drop=h0-h8);
input center initial $ treat $ h0-h8;
array hours h0-h8;
do i=1 to 9;
no_hours=i-1; count=hours(i); output;
end;
datalines;
1 some placebo 1 0 3 0 2 2 4 4 2
1 some treat_a 2 1 0 2 1 2 4 5 1
1 some treat_b 0 0 0 1 0 3 7 6 2
1 some treat_ba 0 0 0 0 1 3 5 4 6
1 lot placebo 6 1 2 2 2 3 7 3 0
1 lot treat_a 6 3 1 2 4 4 7 1 0
1 lot treat_b 3 1 0 4 2 3 11 4 0
1 lot treat_ba 0 0 0 1 1 7 9 6 2
2 some placebo 2 0 2 1 3 1 2 5 4
2 some treat_a 0 0 0 1 1 1 8 1 7
2 some treat_b 0 2 0 1 0 1 4 6 6
2 some treat_ba 0 0 0 1 3 0 4 7 5
2 lot placebo 7 2 3 2 3 2 3 2 2
2 lot treat_a 3 1 0 0 3 2 9 7 1
2 lot treat_b 0 0 0 1 1 5 8 7 4
2 lot treat_ba 0 1 0 0 1 2 8 9 5
3 some placebo 5 0 0 1 3 1 4 4 5
3 some treat_a 1 0 0 1 3 5 3 3 6
3 some treat_b 3 0 1 1 0 0 3 7 11
3 some treat_ba 0 0 0 1 1 4 2 4 13
3 lot placebo 6 0 2 2 2 6 1 2 1
3 lot treat_a 4 2 1 5 1 1 3 2 3
3 lot treat_b 5 0 2 3 1 0 2 6 7
3 lot treat_ba 3 2 1 0 0 2 5 9 4
4 some placebo 1 0 1 1 4 1 1 0 10
4 some treat_a 0 0 0 1 0 2 2 1 13
4 some treat_b 0 0 0 1 1 1 1 5 11
4 some treat_ba 1 0 0 0 0 2 2 2 14
4 lot placebo 4 0 1 3 2 1 1 2 2
4 lot treat_a 0 1 3 1 1 6 1 3 6
4 lot treat_b 0 0 0 0 2 7 2 2 9
4 lot treat_ba 1 0 3 0 1 2 3 4 8
;
proc print data=pain(obs=9);
run;
proc catmod;
weight count;
response 0 .125 .25 .375 .5 .625 .75 .875 1;
model no_hours = center initial treat
treat*initial;
run;
proc catmod;
weight count;
response 0 .125 .25 .375 .5 .625 .75 .875 1;
model no_hours = center initial
treat(initial);
run;
proc catmod;
weight count;
response 0 .125 .25 .375 .5 .625 .75 .875 1;
model no_hours = center initial
treat(initial);
contrast 'lot: a-placebo' treat(initial) -1 1 0 0 0 0 ;
contrast 'lot: b-placebo' treat(initial) -1 0 1 0 0 0 ;
contrast 'lot: ba-placebo' treat(initial) -2 -1 -1 0 0 0 ;
contrast 'lot: ba-a' treat(initial) -1 -2 -1 0 0 0 ;
contrast 'lot: ba-b' treat(initial) -1 -1 -2 0 0 0 ;
contrast 'some:a-placebo' treat(initial) 0 0 0 -1 1 0 ;
contrast 'some:b-placebo' treat(initial) 0 0 0 -1 0 1 ;
contrast 'some:ba-placebo' treat(initial) 0 0 0 -2 -1 -1 ;
contrast 'some:ba-a' treat(initial) 0 0 0 -1 -2 -1 ;
contrast 'some:ba-b' treat(initial) 0 0 0 -1 -1 -2 ;
run;
proc catmod;
weight count;
response 0 .125 .25 .375 .5 .625 .75 .875 1;
model no_hours = center initial
treat(initial);
contrast 'interact:a-placebo' treat(initial) -1 1 0 1 -1 0 ;
contrast 'interact:b-placebo' treat(initial) -1 0 1 1 0 -1 ;
contrast 'interact:ba-placebo' treat(initial) -2 -1 -1 2 1 1 ;
contrast 'average:a-placebo' treat(initial) -1 1 0 -1 1 0 ;
contrast 'average:b-placebo' treat(initial) -1 0 1 -1 0 1 ;
contrast 'average:ba-placebo' treat(initial) -2 -1 -1 -2 -1 -1 ;
contrast 'average:ba-a' treat(initial) -1 -2 -1 -1 -2 -1 ;
contrast 'average:ba-b' treat(initial) -1 -1 -2 -1 -1 -2 ;
contrast 'interaction' treat(initial) -1 1 0 1 -1 0 ,
treat(initial) -1 0 1 1 0 -1 ,
treat(initial) -2 -1 -1 2 1 1 ;
contrast 'treatment effect' treat(initial) -1 1 0 -1 1 0 ,
treat(initial) -1 0 1 -1 0 1 ,
treat(initial) -2 -1 -1 -2 -1 -1 ;
run;
data wbeing;
input #1 b1-b5 _type_ $ _name_ $8. #2 b6-b10;
datalines;
7.24978 7.18991 7.35960 7.31937 7.55184 parms
7.93726 7.92509 7.82815 7.73696 8.16791
0.01110 0.00101 0.00177 -0.00018 -0.00082 cov b1
0.00189 -0.00123 0.00434 0.00158 -0.00050
0.00101 0.02342 0.00144 0.00369 0.25300 cov b2
0.00118 -0.00629 -0.00059 0.00212 -0.00098
0.00177 0.00144 0.01060 0.00157 0.00226 cov b3
0.00140 -0.00088 -0.00055 0.00211 0.00239
-0.00018 0.00369 0.00157 0.02298 0.00918 cov b4
-0.00140 -0.00232 0.00023 0.00066 -0.00010
-0.00082 0.00253 0.00226 0.00918 0.01921 cov b5
0.00039 0.00034 -0.00013 0.00240 0.00213
0.00189 0.00118 0.00140 -0.00140 0.00039 cov b6
0.00739 0.00019 0.00146 -0.00082 0.00076
-0.00123 -0.00629 -0.00088 -0.00232 0.00034 cov b7
0.00019 0.01172 0.00183 0.00029 0.00083
0.00434 -0.00059 -0.00055 0.00023 -0.00013 cov b8
0.00146 0.00183 0.01050 -0.00173 0.00011
0.00158 0.00212 0.00211 0.00066 0.00240 cov b9
-0.00082 0.00029 -0.00173 0.01335 0.00140
-0.00050 -0.00098 0.00239 -0.00010 0.00213 cov b10
0.00076 0.00083 0.00011 0.00140 0.01430
;
proc catmod data=wbeing;
response read b1-b10;
factors sex $ 2, age $ 5 /
_response_ = sex|age
profile = (female '25-34',
female '35-44',
female '45-54',
female '55-64',
female '65-74',
male '25-34',
male '35-44',
male '45-54',
male '55-64',
male '65-74');
model _f_ = _response_;
run;
proc catmod data=wbeing;
response read b1-b10;
factors sex $ 2, age $ 5 /
_response_ = sex age
profile = (male '25-34' ,
male '35-44',
male '45-54' ,
male '55-64',
male '65-74' ,
female '25-34',
female '35-44',
female '45-54',
female '55-64' ,
female '65-74');
model _f_ = _response_;
contrast '25-34 vs. 35-44' all_parms 0 0 1 -1 0 0;
contrast '25-34 vs. 45,54' all_parms 0 0 1 0 -1 0;
contrast '25-34 vs. 55,64' all_parms 0 0 1 0 0 -1;
contrast '25-34 vs. 65,74' all_parms 0 0 2 1 1 1;
run;
proc catmod data=wbeing;
response read b1-b10;
factors sex $ 2, age $ 5 /
_response_ = sex age
profile = (male '25-34' ,
male '35-44',
male '45-54' ,
male '55-64',
male '65-74' ,
female '25-34',
female '35-44',
female '45-54',
female '55-64' ,
female '65-74');
model _f_ = _response_;
contrast '25-64 the same' all_parms 0 0 1 -1 0 0,
all_parms 0 0 1 0 -1 0,
all_parms 0 0 1 0 0 -1;
run;
proc catmod data=wbeing;
response read b1-b10;
factors sex $ 2, age $ 5 /
_response_ = sex age
profile = (female '25-34',
female '35-44',
female '45-54',
female '55-64',
female '65-74',
male '25-34',
male '35-44',
male '45-54',
male '55-64',
male '65-74');
model _f_ = ( 1 0 0 ,
1 0 0 ,
1 0 0 ,
1 0 0 ,
1 0 1 ,
1 1 0 ,
1 1 0 ,
1 1 0 ,
1 1 0 ,
1 1 1 ) (1='Intercept', 2='Sex', 3='65-74')
/ pred;
proc freq data=cpain;
weight count;
tables dstatus*invest*treat*status/ measures;
run;
data MannWhitney;
input b1-b8 _type_ $ _name_ $8.;
datalines;
.6000 .6011 .6042 .8389 .5130 .5947 .5000 .4922 parms
.03091 .0000 .0000 .0000 .0000 .0000 .0000 .0000 cov b1
.0000 .00918 .0000 .0000 .0000 .0000 .0000 .0000 cov b2
.0000 .0000 .3280 .0000 .0000 .0000 .0000 .0000 cov b3
.0000 .0000 .0000 .0084 .0000 .0000 .0000 .0000 cov b4
.0000 .0000 .0000 .0000 .0129 .0000 .0000 .0000 cov b5
.0000 .0000 .0000 .0000 .0000 .0093 .0000 .0000 cov b6
.0000 .0000 .0000 .0000 .0000 .0000 .0101 .0000 cov b7
.0000 .0000 .0000 .0000 .0000 .0000 .0000 .0112 cov b8
;
proc catmod data=MannWhitney;
response read b1-b8;
factors diagnosis $ 4 , invest $ 2 /
_response_ = diagnosis invest
profile = (I A,
I B,
II A,
II B,
III A,
III B,
IV A,
IV B);
model _f_ = _response_ / cov;
run;
proc catmod data=MannWhitney;
response read b1-b8;
factors diagnosis $ 4 , invest $ 2 /
_response_ = diagnosis invest
profile = (I A,
I B,
II A,
II B,
III A,
III B,
IV A,
IV B);
model _f_ = ( 1 ,
1 ,
1 ,
1 ,
1 ,
1 ,
1 ,
1 ) ;
run;
Statistics and Operations
Research
