/*--------
|fitting logitstic regression to smoking data
-------------*/
options ls=90 nocenter nodate;

/**********************************************************/
/* fitting to 2x2 table */
data smoke;
input s $ y n ;
cards;
smoke 816 4019
nosmoke 188 1356
;
proc logistic data=smoke descending;
class s (ref=first) / param=ref;
model y/n = s /scale=none;
output out=predict pred=prob;
title 'Logistic regression for 2x2 table';
run;
proc print data=predict;
run;

/***********************************************************/
/* fitting to 2x3 table */
/* specifying the reference level of a string variable */
data smoke;
input s $ y n ;
cards;
both 400 1780
one 416 2239
neither 188 1356
;
proc logistic data=smoke descending;
class s (ref='neither') / order=data param=ref;
model y/n = s /scale=none lackfit;
output out=predict pred=prob;
title 'Logistic regression for 2x3 table';
run;
proc print data=predict;
run;


/***********************************************************/
/* fitting 2x3 table */
/* anther version with extra diagnostics */
data smoke;
input s $ y n ;
cards;
2 400 1780
1 416 2239
0 188 1356
;
proc logistic data=smoke descending;
class s (ref=first)/ param=ref;
model y/n = s / scale=none lackfit;
output out=predict pred=prob reschi=pearson resdev=deviance;
title 'Logistic regression fro 2x3 tables with residuals';
run;

data diagnostics;
set predict;
shat = n*prob;
fhat = n*(1-prob);
run;

proc print data=diagnostics;
var s y n prob shat fhat pearson deviance;
title 'Logistic regression diagnostics for 2x3 table';
run;

/***********************************************************/
/* fitting 2x3 table */
/* anther version with dummy variables and extra diagnostics */
data smoke;
input s1 s2 $ y n ;
cards;
0 1 400 1780
1 0 416 2239
0 0 188 1356
;
proc logistic data=smoke descending;
class s1 s2 (ref=first)/ param=ref;
model y/n = s1 s2 / scale=none lackfit;
output out=predict pred=prob reschi=pearson resdev=deviance;
run;

proc logistic data=smoke descending;
class s1 s2 (ref=first)/ param=ref;
model y/n = s1 / scale=none lackfit;
output out=predict pred=prob reschi=pearson resdev=deviance;
run;

proc logistic data=smoke descending;
class s1 s2 (ref=first)/ param=ref;
model y/n = s2 / scale=none lackfit;
output out=predict pred=prob reschi=pearson resdev=deviance;
title 'Logistic regression fro 2x3 tables with residuals';
run;