Categorical Data Analysis: Chapter 14
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 design;
run;
proc catmod;
weight count;
response means;
model periods = sex residence;
run;
proc catmod;
weight count;
response means;
model periods = sex;
run;
proc catmod;
population sex residence;
weight count;
response means;
model periods = sex;
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
run;
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;
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;
data wbeing(type=est);
input b1-b5 _type_ $ _name_ $ b6-b10 #2;
datalines;
7.93726 7.92509 7.82815 7.73696 8.16791 parms .
7.24978 7.18991 7.35960 7.31937 7.55184
0.00739 0.00019 0.00146 -0.00082 0.00076 cov b1
0.00189 0.00118 0.00140 -0.00140 0.00039
0.00019 0.01172 0.00183 0.00029 0.00083 cov b2
-0.00123 -0.00629 -0.00088 -0.00232 0.00034
0.00146 0.00183 0.01050 -0.00173 0.00011 cov b3
0.00434 -0.00059 -0.00055 0.00023 -0.00013
-0.00082 0.00029 -0.00173 0.01335 0.00140 cov b4
0.00158 0.00212 0.00211 0.00066 0.00240
0.00076 0.00083 0.00011 0.00140 0.01430 cov b5
-0.00050 -0.00098 0.00239 -0.00010 0.00213
0.00189 -0.00123 0.00434 0.00158 -0.00050 cov b6
0.01110 0.00101 0.00177 -0.00018 -0.00082
0.00118 -0.00629 -0.00059 0.00212 -0.00098 cov b7
0.00101 0.02342 0.00144 0.00369 0.00253
0.00140 -0.00088 -0.00055 0.00211 0.00239 cov b8
0.00177 0.00144 0.01060 0.00157 0.00226
-0.00140 -0.00232 0.00023 0.00066 -0.00010 cov b9
-0.00018 0.00369 0.00157 0.02298 0.00918
0.00039 0.00034 -0.00013 0.00240 0.00213 cov b10
-0.00082 0.00253 0.00226 0.00918 0.01921
;
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_ / design;
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;
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;
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 );
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_ = ( 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;
data cpain;
data cpain;
input center $ diagnosis $ treat $ status $ count @@;
datalines;
I A active excellent 1 I A active good 3 I A active moderate 2
I A active fair 5 I A active poor 1
I A placebo excellent 2 I A placebo good 4 I A placebo moderate 3
I A placebo fair 4 I A placebo poor 3
I B active excellent 3 I B active good 10 I B active moderate 1
I B active fair 4 I B active poor 2
I B placebo excellent 2 I B placebo good 4 I B placebo moderate 1
I B placebo fair 5 I B placebo poor 2
I C active excellent 6 I C active good 1 I C active moderate 1
I C active fair 1 I C active poor 0
I C placebo excellent 0 I C placebo good 5 I C placebo moderate 1
I C placebo fair 1 I C placebo poor 3
I D active excellent 3 I D active good 5 I D active moderate 1
I D active fair 6 I D active poor 1
I D placebo excellent 3 I D placebo good 3 I D placebo moderate 2
I D placebo fair 4 I D placebo poor 5
II A active excellent 0 II A active good 4 II A active moderate 3
II A active fair 1 II A active poor 8
II A placebo excellent 0 II A placebo good 3 II A placebo moderate 3
II A placebo fair 0 II A placebo poor 5
II B active excellent 2 II B active good 3 II B active moderate 3
II B active fair 0 II B active poor 2
II B placebo excellent 1 II B placebo good 8 II B placebo moderate 0
II B placebo fair 0 II B placebo poor 5
II C active excellent 2 II C active good 2 II C active moderate 1
II C active fair 0 II C active poor 1
II C placebo excellent 1 II C placebo good 1 II C placebo moderate 0
II C placebo fair 1 II C placebo poor 1
II D active excellent 0 II D active good 1 II D active moderate 2
II D active fair 2 II D active poor 3
II D placebo excellent 1 II D placebo good 1 II D placebo moderate 1
II D placebo fair 0 II D placebo poor 7
run;
proc freq data=cpain order=data;
weight count;
tables center*diagnosis*treat*status/ measures;
run;
proc sort data=cpain; by diagnosis;
proc freq data=cpain order=data; by diagnosis;
weight count;
where (center='II' and (diagnosis='A' or diagnosis='B'));
tables treat*status/ measures;
run;
data MannWhitney;
input b1-b8 _type_ $ _name_ $8.;
datalines;
.4922 .5946 .8389 .6011 .4688 .5286 .6042 .6000 parms
.0112 .0000 .0000 .0000 .0000 .0000 .0000 .0000 cov b1
.0000 .0092 .0000 .0000 .0000 .0000 .0000 .0000 cov b2
.0000 .0000 .0084 .0000 .0000 .0000 .0000 .0000 cov b3
.0000 .0000 .0000 .0092 .0000 .0000 .0000 .0000 cov b4
.0000 .0000 .0000 .0000 .0110 .0000 .0000 .0000 cov b5
.0000 .0000 .0000 .0000 .0000 .0133 .0000 .0000 cov b6
.0000 .0000 .0000 .0000 .0000 .0000 .0328 .0000 cov b7
.0000 .0000 .0000 .0000 .0000 .0000 .0000 .0158 cov b8
run;
proc catmod data=MannWhitney;
response read b1-b8;
factors center $ 2 , diagnosis $ 4 /
_response_ = center diagnosis
profile = (I A,
I B,
I C,
I D,
II A,
II B,
II C,
II D);
model _f_ = _response_ / cov;
run;
model _f_ =( 1,
1,
1,
1,
1,
1,
1,
1 );
data drug;
input druga $ drugb $ drugc $ count;
datalines;
F F F 6
F F U 16
F U F 2
F U U 4
U F F 2
U F U 4
U U F 6
U U U 6
;
proc catmod;
weight count;
response marginals;
model druga*drugb*drugc=_response_ / oneway cov;
repeated drug 3 / _response_=drug;
run;
model druga*drugb*drugc=_response_ / oneway design cov;
repeated drug 3 / _response_=drug;
repeated drug;
contrast 'A versus C' _response_ 2 1;
contrast 'A versus C' all_parms 0 2 1;
data vision;
input gender $ right left count;
datalines;
F 1 1 1520
F 1 2 266
F 1 3 124
F 1 4 66
F 2 1 234
F 2 2 1512
F 2 3 432
F 2 4 78
F 3 1 117
F 3 2 362
F 3 3 1772
F 3 4 205
F 4 1 36
F 4 2 82
F 4 3 179
F 4 4 492
M 1 1 821
M 1 2 112
M 1 3 85
M 1 4 35
M 2 1 116
M 2 2 494
M 2 3 145
M 2 4 27
M 3 1 72
M 3 2 151
M 3 3 583
M 3 4 87
M 4 1 43
M 4 2 34
M 4 3 106
M 4 4 331
;
proc catmod;
weight count;
response marginals;
model right*left=gender _response_(gender='F')
_response_(gender='M') / design;
repeated eye 2;
run;
contrast 'Interaction' all_parms 0 0 0 0 0 0 1 0 0 -1 0 0,
all_parms 0 0 0 0 0 0 0 1 0 0 -1 0,
all_parms 0 0 0 0 0 0 0 0 1 0 0 -1;
run;
model right*left=gender|_response_;
repeated eye 2;
run;
proc catmod;
weight count;
response means;
model right*left=gender _response_(gender='F')
_response_(gender='M') / noprofile design;
repeated eye;
run;
contrast 'Interaction' all_parms 0 0 1 -1;
run;
data wheeze;
input wheeze9 $ wheeze10 $ wheeze11 $ wheeze12 $ count;
datalines;
Present Present Present Present 94
Present Present Present Absent 30
Present Present Absent Present 15
Present Present Absent Absent 28
Present Absent Present Present 14
Present Absent Present Absent 9
Present Absent Absent Present 12
Present Absent Absent Absent 63
Absent Present Present Present 19
Absent Present Present Absent 15
Absent Present Absent Present 10
Absent Present Absent Absent 44
Absent Absent Present Present 17
Absent Absent Present Absent 42
Absent Absent Absent Present 35
Absent Absent Absent Absent 572
;
proc catmod order=data;
weight count;
response marginals;
model wheeze9*wheeze10*wheeze11*wheeze12=_response_ / oneway;
repeated age;
run;
proc catmod order=data;
weight count;
response logits;
model wheeze9*wheeze10*wheeze11*wheeze12=(1 9,
1 10,
1 11,
1 12)
(1='Intercept',
2='Linear Age') / noprofile design;
run;