Categorical Data Analysis: Chapter 11
options nodate nonumber ps=200 ls=80 formdlim=' ';
data bacteria;
input dose status $ count @@;
ldose=log(dose);
sq_ldose=ldose*ldose;
datalines;
1200 dead 0 1200 alive 5
12000 dead 0 12000 alive 5
120000 dead 2 120000 alive 3
1200000 dead 4 1200000 alive 2
12000000 dead 5 12000000 alive 1
120000000 dead 5 120000000 alive 0
;
proc print;
run;
proc logistic data=bacteria descending;
freq count;
model status = ldose sq_ldose / scale=none aggregate
selection=forward start=1 details covb;
run;
data assay;
input drug $ dose status $ count;
int_n=(drug='n');
int_s=(drug='s');
ldose=log(dose);
ldose_n=int_n*ldose;
ldose_s=int_s*ldose;
sqldose_n=int_n*ldose*ldose;
sqldose_s=int_s*ldose*ldose;
datalines;
n 0.01 dead 0
n 0.01 alive 30
n .03 dead 1
n .03 alive 29
n .10 dead 1
n .10 alive 9
n .30 dead 1
n .30 alive 9
n 1.00 dead 4
n 1.00 alive 6
n 3.00 dead 4
n 3.00 alive 6
n 10.00 dead 5
n 10.00 alive 5
n 30.00 dead 7
n 30.00 alive 3
s .30 dead 0
s .30 alive 10
s 1.00 dead 0
s 1.00 alive 10
s 3.00 dead 1
s 3.00 alive 9
s 10.00 dead 4
s 10.00 alive 6
s 30.00 dead 5
s 30.00 alive 5
s 100.00 dead 8
s 100.00 alive 2
;
proc logistic data=assay descending;
freq count;
model status = int_n int_s ldose_n ldose_s
sqldose_n sqldose_s / noint
scale=none aggregate
start=4 selection=forward details;
eq_slope: test ldose_n=ldose_s;
run;
proc logistic data=assay descending outest=estimate
(drop= intercept _link_ _lnlike_) covout;
freq count;
model status = int_n int_s ldose /
noint scale=none aggregate covb;
run;
proc iml;
use estimate;
start fieller;
title 'Confidence Intervals';
use estimate;
read all into beta where (_type_='PARMS');
beta=beta`;
read all into cov where (_type_='COV');
ratio=(k`*beta) / (h`*beta);
a=(h`*beta)**2-(3.84)*(h`*cov*h);
b=2*(3.84*(k`*cov*h)-(k`*beta)*(h`*beta));
c=(k`*beta)**2 -(3.84)*(k`*cov*k);
disc=((b**2)-4*a*c);
if (disc<=0 | a<=0) then do;
print "confidence interval can't be computed", ratio;
stop; end;
sroot=sqrt(disc);
l_b=((-b)-sroot)/(2*a);
u_b=((-b)+sroot)/(2*a);
interval=l_b||u_b;
lname={"l_bound", "u_bound"};
print "95 % ci for ratio based on fieller", ratio interval[colname=lname];
finish fieller;
k={ 1 -1 0 }`;
h={ 0 0 1 }`;
run fieller;
k={-1 0 0 }`;
h={ 0 0 1 }`;
run fieller;
k={ 0 -1 0 }`;
h={ 0 0 1 }`;
run fieller;
data adverse;
input diagnos $ dose status $ count @@;
i_diagII=(diagnos='II');
i_diagI= (diagnos='I');
doseI=i_diagI*dose;
doseII=i_diagII*dose;
diagdose=i_diagII*dose;
if doseI > 0 then ldoseI=log(doseI); else ldoseI=0;
if doseII > 0 then ldoseII=log(doseII); else ldoseII=0;
datalines;
I 1 adverse 3 I 1 no 26
I 5 adverse 7 I 5 no 26
I 10 adverse 10 I 10 no 22
I 12 adverse 14 I 12 no 18
I 15 adverse 18 I 15 no 14
II 1 adverse 6 II 1 no 26
II 5 adverse 20 II 5 no 12
II 10 adverse 26 II 10 no 6
II 12 adverse 28 II 12 no 4
II 15 adverse 31 II 15 no 1
;
proc freq data=adverse;
weight count;
tables dose*status diagnos*status diagnos*dose*status /
nopct nocol cmh;
run;
proc logistic data=adverse outest=estimate
(drop= intercept _link_ _lnlike_) covout;
freq count;
model status = i_diagI i_diagII doseI doseII /
noint scale=none aggregate;
eq_slope: test doseI=doseII;
run;
proc logistic data=adverse;
freq count;
model status = i_diagI i_diagII ldoseI ldoseII /
noint scale=none aggregate;
eq_slope: test ldoseI=ldoseII;
run;
proc iml;
start fieller;
title 'Confidence Intervals';
use estimate;
read all into beta where (_type_='PARMS');
beta=beta`;
read all into cov where (_type_='COV');
ratio=(k`*beta) / (h`*beta);
a=(h`*beta)**2-(3.84)*(h`*cov*h);
b=2*(3.84*(k`*cov*h)-(k`*beta)*(h`*beta));
c=(k`*beta)**2 -(3.84)*(k`*cov*k);
disc=((b**2)-4*a*c);
if (disc<=0 | a<=0) then do;
print "confidence interval can't be computed", ratio;
stop; end;
sroot=sqrt(disc);
l_b=((-b)-sroot)/(2*a);
u_b=((-b)+sroot)/(2*a);
interval=l_b||u_b;
lname={"l_bound", "u_bound"};
print "95 % ci for ratio based on fieller", ratio interval[colname=lname];
finish fieller;
k={ -1 0 0 0}`;
h={ 0 0 1 0}`;
run fieller;
k={ 0 -1 0 0}`;
h={ 0 0 0 1}`;
run fieller;
Statistics and Operations Research
