Miscellaneous Examples for PROC LOGISTIC
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: LOGIMISC */
/* TITLE: Miscellaneous Examples for PROC LOGISTIC */
/* PRODUCT: STAT */
/* SYSTEM: ALL */
/* KEYS: logistic regression analysis, */
/* PROCS: LOGISTIC */
/* DATA: */
/* */
/* SUPPORT: YCS UPDATE: */
/* REF: */
/* chapter */
/* MISC: */
/* */
/****************************************************************/
/*****************************************************************
Example 1. Using Proc Logistic for Poisson Regression
*****************************************************************/
/*
The data, taken from Koch et. al (1986), consist of counts of
new melanoma cases reported during 1969-1971 among white males
and the corresponding estimated populations at risk for six
age groups and two geographic regions. Possion regression
was proposed to study the vital statistics factors on new
melanoma cases. Since the event frequencies are relatively
small as compared to the corresponding population sizes, the
logistic regression estimates are essentially the same as
the Poisson regression estimates.
*/
title1 'Example 1. Poisson Regression';
data;
drop case_n case_s pop_n pop_s;
length age_grp $ 5;
input age_grp $ case_n case_s pop_n pop_s;
age_grp2 = ( age_grp = '35-44' );
age_grp3 = ( age_grp = '45-54' );
age_grp4 = ( age_grp = '55-64' );
age_grp5 = ( age_grp = '65-74' );
age_grp6 = ( age_grp = '>=75' );
case = case_n; pop = pop_n; region = 0; output;
case = case_s; pop = pop_s; region = 1; output;
label region = '0 = NORTHERN 1 = SOUTHERN';
datalines;
<35 61 64 2880262 1074246
35-44 76 75 564535 220407
45-54 98 68 592983 198119
55-64 104 63 450740 134084
65-74 63 45 270908 70708
>=75 80 27 161850 34233
;
proc logistic;
model case/pop=age_grp2-age_grp6 region;
run;
/*****************************************************************
Example 2. Overdispersion in Quantal Assay Data
*****************************************************************/
/*
In bioassay, one is often interested in the median effect dose,
ED50, which is the concentration of a chemical that is expected
to produce a response in 50% of the subjects exposed to it.
By fitting simple models such as logit, probit, or complementary
log-log, ED50 is a simple function of the regression coeffi-
cients. However, the basic underlying assumption that the
quantal assay data follow a binomial distribution may not
always be appropriate. This example presents an analysis that
incorporates extra-binomial variation to the data.
The following SAS statements create the data set ASSAY, which
contains a set of dose response data. Variable LOGCONC
represents the log concentration of a chemical, N represents
the number of subjects dosed, and R represents the number of
subjects responded.
PROC LOGISTIC is first used to fit a logit model to the quantal
assay data. The option SCALE=NONE is specified to display the
goodness-of-fit statistics. Linear predictor values and
Pearson residuals are output to data set OUT1, which is
subsequently used in PROC GPLOT to display the adequacy of
the fit.
*/
title1 'Example 2. Overdispersion in Quantal Assay Data';
data assay;
input logconc r n;
logit=log(r/(n-r));
datalines;
2.68 10 31
2.76 17 30
2.82 12 31
2.90 7 27
3.02 23 26
3.04 22 30
3.13 29 31
3.20 29 30
3.21 23 30
;
proc logistic data=assay outest=dn covout;
model r/n= logconc / scale=none;
output out=out1 xbeta=xb reschi=reschi;
run;
proc sgplot data=out1 noautolegend;
title 'Plot of Data with Fitted Line';
scatter y=logit x=logconc;
series y=xb x=logconc;
label logit = 'Logit of Proportion Responded'
logconc = 'Log Concentration';
run;
proc sgplot data=out1;
title 'Residual Plot';
scatter y=reschi x=logconc;
label logit = 'Pearson Residual'
logconc = 'Log Concentration';
run;
/*
Both goodness-of-fits statistics are statistically significant,
indicating a lack of fit of the model. However, there is not
a clear pattern of Pearson residuals and the fit of the
model appears to be adequate. This suggests that extra-binomial
variation is responsible for the poor fit.
The Williams method of modeling overdispersion is then used.
This is requested by specifying the SCALE=WILLIAMS option.
*/
title2 'Scale=Williams';
proc logistic data=assay outest=dw covout;
model r/n= logconc / scale=williams;
run;
/*
With log of concentration as the only explanatory variable,
ED50 is estimated by exp(-intercept/slope). A 95% confidence
interval of ED50 is obtained by exponentiating the 95%
confidence interval for (-intercept/slope). The latter is
conveniently calculated using Fieller's theorem.
In a single DATA step, you can compute the confidence limits
for ED50 from the OUTEST= data set that contains both
parameter estimates and covariance matrix. The following SAS
statements compute the 95% confidence interval for ED50 for
each of the SCALE= option available in PROC LOGISTIC. An
OUTEST= data set that contains both estimates and dispersion
is first output for each SCALE= option. The four OUTEST= data
sets are then concatenated in data set TWO and confidence
limits are subsequently calculated for each data set.
*/
title2 'Scale=Pearson';
proc logistic data=assay outest=dp covout noprint;
model r/n= logconc / scale=p;
run;
title2 'Scale=Deviance';
proc logistic data=assay outest=dd covout noprint;
model r/n= logconc / scale=d;
run;
data two(keep=scale ed50 lower upper);
retain b0 b1 v00 v01 v11;
set dn(in=dn) dw(in=dw) dp(in=dp) dd(in=dd);
if _name_='r' then do;
b0= INTERCEPT;
b1= logconc;
end;
if _name_='Intercept' then do;
v00=INTERCEPT;
v01=logconc;
end;
if _name_='logconc ' then do;
if dn=1 then scale='NONE ';
else if dw=1 then scale='WILLIAMS';
else if dp=1 then scale='PEARSON ';
else if dd=1 then scale='DEVIANCE';
v11=logconc;
theta=-b0/b1;
t975=tinv(.975,7);
c=t975*t975*v11/(b1*b1);
term1=theta+(c/(1-c))*(theta + v01/v11);
term2=t975/(b1*(1-c));
term3=(v00+2*theta*v01+theta*theta*
v11-c*(v00-v01*v01/v11));
term4=term2*sqrt(term3);
l1=term1-term4;
u1=term1+term4;
ed50= exp(theta);
lower=exp(l1);
upper=exp(u1);
output;
end;
run;
title2 'ED50 and 95% Confidence Interval';
proc print data=two;
id scale;
run;
/*****************************************************************
Example 3. 1:M Matching
*****************************************************************/
/*
Consider the special case that (i) there is only one case in
each matched set and the number of controls are the same for all
sets, and (ii) there is only a single binary covariable with
relative risk exp(beta). The Samsu consumption data in Breslow
(1982) contains 80 matched sets. Each matched set has one case
and four controls. The distribution according to Samsu exposure
is given in the following table:
Case Total Number (Case + Controls)
0 1 2 3 4 5
Exposed . 5 19 10 6 0
Not Exposed 10 15 8 7 0 .
Total 10 20 27 17 6 .
Let n_0m be the number of sets in which exactly m controls are
exposed and the case is not; and let n_1m be the number of sets
in which the case and m controls are exposed. The likelihood of
the conditional logistic regression model is proportional to the
product of the likelihood of the binomial distribution
B(N_m, theta_m), where
N_m = n_0m + n_1m
and
theta_m = m * exp(beta) / [m * exp(beta) + 5 - m]
Since logit(theta_m) = log( m/(5-m) + beta ), beta can be
estimated as a parameter in a logistic regression model with
no intercept and an offset value of m/(5-m).
This analysis is performed with the PROC LOGISTIC invocation
below.
*/
title1 'Example 3. 1:M Matching';
data samsu;
input m r n;
samsu = 1;
off= log(m/(5-m));
datalines;
1 5 20
2 19 27
3 10 17
4 6 6
;
proc logistic data=samsu;
model r/n = samsu / offset=off noint;
run;
/*
Note that you can perform the same analysis using the STRATA statement to
perform a conditional logistic regression.
First create the SAMSU data set as a single record for each subject.
*/
data one;
input exp1-exp5 freq;
datalines;
1 0 0 0 0 5
1 1 0 0 0 19
1 1 1 0 0 10
1 1 1 1 0 6
1 1 1 1 1 0
0 0 0 0 0 10
0 1 0 0 0 15
0 1 1 0 0 8
0 1 1 1 0 7
0 1 1 1 1 0
;
data two;
set one;
retain id 0;
do i=1 to freq;
id=id+1;
output;
end;
run;
data three;
set two;
case=1; samsu=exp1; output;
case=0; samsu=exp2; output;
samsu=exp3; output;
samsu=exp4; output;
samsu=exp5; output;
keep case samsu id;
run;
/*
Then run a conditional logistic regression.
*/
proc logistic data=three;
model case(event='1')=samsu;
strata id;
run;
