Categorical Data Analysis: Chapter 8
data coronary;
input sex ecg ca count @@;
datalines;
0 0 0 11 0 0 1 4
0 1 0 10 0 1 1 8
1 0 0 9 1 0 1 9
1 1 0 6 1 1 1 21
;
proc logistic;
freq count;
model ca(event='1')=sex ecg / scale=none aggregate;
run;
proc logistic descending;
freq count;
model ca=sex ecg;
output out=predict pred=prob;
run;
proc print data=predict;
run;
ods select FitStatistics ParameterEstimates;
proc logistic descending;
freq count;
model ca=sex ecg sex*ecg;
run;
data sentence;
input type $ prior $ sentence $ count @@;
datalines;
nrb some y 42 nrb some n 109
nrb none y 17 nrb none n 75
other some y 33 other some n 175
other none y 53 other none n 359
;
proc logistic;
class type prior(ref=first) / param=ref;
freq count;
model sentence(event='y') = type prior / scale=none aggregate;
run;
ods select GoodnessOfFit;
proc logistic;
class type prior (ref=first) / param=ref;
freq count;
model sentence(event='y') = type / scale=none aggregate=(type prior);
run;
ods select ClassLevelInfo GoodnessOfFit
ParameterEstimates OddsRatios;
proc logistic data=sentence;
class type prior(ref='none');
freq count;
model sentence(event='y')= type prior / scale=none aggregate;
run;
data uti;
input diagnosis : $13. treatment $ response $ count @@;
datalines;
complicated A cured 78 complicated A not 28
complicated B cured 101 complicated B not 11
complicated C cured 68 complicated C not 46
uncomplicated A cured 40 uncomplicated A not 5
uncomplicated B cured 54 uncomplicated B not 5
uncomplicated C cured 34 uncomplicated C not 6
;
proc logistic;
freq count;
class diagnosis treatment /param=ref;
model response = diagnosis|treatment;
run;
proc logistic;
freq count;
class diagnosis treatment;
model response = diagnosis treatment /
scale=none aggregate;
run;
data;
p=1-probchi(2.515,2);
put p= ;
run;
ods graphics on;
proc logistic data=uti
plots(only)=(effect(clband yrange=(.5,1) x=treatment*diagnosis)
oddsratio(logbase=2));
freq count;
class diagnosis treatment;
model response = diagnosis treatment /
scale=none aggregate;
run;
ods graphics off;
ods select ContrastTest ContrastEstimate;
proc logistic data=uti;
freq count;
class diagnosis treatment /param=ref;
model response = diagnosis treatment;
contrast 'A versus B' treatment 1 -1
/ estimate=exp;
contrast 'A' treatment 1 0;
contrast 'joint test' treatment 1 0,
treatment 0 1;
run;
ods graphics on;
proc logistic data=uti plots(only)=effect(x=treatment
sliceby=diagnosis clbar
connect yrange=(0.5));
freq count;
class diagnosis treatment /param=ref;
model response = diagnosis treatment;
oddsratio treatment / cl=pl;
oddsratio treatment / cl=pl;
run;
ods graphics off;
data byss;
input work $ years $ smoke $ 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
not <10 yes yes 14 not <10 yes no 1340
not <10 no yes 12 not <10 no no 1004
not >=10 yes yes 24 not >=10 yes no 1360
not >=10 no yes 10 not >=10 no no 986
;
proc logistic data=byss;
freq count;
class work years(ref=first) smoke(ref=first) /param=ref;
model status(event=last) = work|years|smoke@2 /
scale=none aggregate;
run;
ods graphics on;
proc logistic plots(only)=(oddsratio(logbase=2));
freq count;
class work years(ref=first) smoke(ref=first) /param=ref;
model status(event=last) = work years smoke
work*years work*smoke
/scale=none aggregate;
effectplot interaction (x=work) / at(smoke=all years=all) link noobs;
oddsratios work;
run;
ods graphics off;
data coronary;
input sex ecg age ca @@ ;
datalines;
0 0 28 0 1 0 42 1 0 1 46 0 1 1 45 0
0 0 34 0 1 0 44 1 0 1 48 1 1 1 45 1
0 0 38 0 1 0 45 0 0 1 49 0 1 1 45 1
0 0 41 1 1 0 46 0 0 1 49 0 1 1 46 1
0 0 44 0 1 0 48 0 0 1 52 0 1 1 48 1
0 0 45 1 1 0 50 0 0 1 53 1 1 1 57 1
0 0 46 0 1 0 52 1 0 1 54 1 1 1 57 1
0 0 47 0 1 0 52 1 0 1 55 0 1 1 59 1
0 0 50 0 1 0 54 0 0 1 57 1 1 1 60 1
0 0 51 0 1 0 55 0 0 2 46 1 1 1 63 1
0 0 51 0 1 0 59 1 0 2 48 0 1 2 35 0
0 0 53 0 1 0 59 1 0 2 57 1 1 2 37 1
0 0 55 1 1 1 32 0 0 2 60 1 1 2 43 1
0 0 59 0 1 1 37 0 1 0 30 0 1 2 47 1
0 0 60 1 1 1 38 1 1 0 34 0 1 2 48 1
0 1 32 1 1 1 38 1 1 0 36 1 1 2 49 0
0 1 33 0 1 1 42 1 1 0 38 1 1 2 58 1
0 1 35 0 1 1 43 0 1 0 39 0 1 2 59 1
0 1 39 0 1 1 43 1 1 0 42 0 1 2 60 1
0 1 40 0 1 1 44 1
;
proc logistic descending;
model ca=sex ecg age ecg*ecg age*age
sex*ecg sex*age ecg*age /
selection=forward include=3 details lackfit;
run;
proc logistic descending;
model ca=sex ecg age;
units age=10;
run;
data uti2;
input diagnosis : $13. treatment $ response trials;
datalines;
complicated A 78 106
complicated B 101 112
complicated C 68 114
uncomplicated A 40 45
uncomplicated B 54 59
uncomplicated C 34 40
;
ods graphics on;
proc logistic data=uti2 plots(label)=influence(unpack stdres);
class diagnosis treatment / param=ref;
model response/trials = diagnosis treatment;
run;
ods graphics off;
data preclinical;
input treatA $ treatB $ response $ count @@;
datalines;
no no yes 0 no no no 0
yes no yes 8 yes no no 0
no yes yes 0 no yes no 2
yes yes yes 21 yes yes no 6
;
proc logistic descending;
freq count;
class treatA treatB /param=ref reference=first;
model response = treatA treatB;
run;
proc logistic descending exactonly;
freq count;
class treatA treatB /param=ref reference=first;
model response = treatA treatB / alpha=.025;
exact treatA treatB /joint estimate=both;
run;
data complete;
input gender region count response @@;
datalines;
0 0 0 1 0 0 5 0
0 1 1 1 0 1 0 0
1 0 0 1 1 0 175 0
1 1 53 1 1 1 0 0
;
proc logistic data=complete descending;
freq count;
model response = gender region / firth clparm=pl;
exact gender region;
run;
data liver;
input time $ group $ status $ count @@;
datalines;
early antidote severe 6 early antidote not 12
early control severe 6 early control not 2
delayed antidote severe 3 delayed antidote not 4
delayed control severe 3 delayed control not 0
late antidote severe 5 late antidote not 1
late control severe 6 late control not 0
;
proc logistic descending;
freq count;
class time(ref='early') group(ref='control') /param=ref;
model status = time group / clparm=wald;
run;
proc logistic descending exactonly;
freq count;
class time(ref='early') group(ref='control') /param=ref;
model status = time group / clparm=wald;
exact 'Model 1' intercept time group /
estimate=both;
exact 'Joint Test' time group / jointonly;
run;
data exercise;
input location $ program $ outcome $ count @@;
datalines;
downtown office good 12 downtown office not 5
downtown home good 3 downtown home not 5
satellite office good 6 satellite office not 1
satellite home good 1 satellite home not 3
;
proc logistic;
freq count;
class location program(ref=first) /param=ref;
model outcome = location program;
exact program / estimate=both;
run;
proc genmod data=uti;
freq count;
class diagnosis treatment;
model response = diagnosis treatment /
link=logit dist=binomial type3 aggregate;
run;
proc genmod data=uti;
freq count;
class diagnosis treatment;
model response = diagnosis treatment /
link=logit dist=binomial;
contrast 'treatment' treatment 1 0 -1 ,
treatment 0 1 -1;
contrast 'A-B' treatment 1 -1 0;
contrast 'A-C' treatment 1 0 -1;
run;