%let location=C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\Chapter14;
PROC IMPORT OUT= WORK.program 
            DATAFILE= "&location\chapter14.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH14TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title "Programming task example. Refer to Table 14.1, page 566";run;
proc logistic;
model success(descending)=exper/influence iplots clodds=both clparm=both corrb covb;
output out=fitted predicted=plogist l=llogist u=ulogist;
title "logistic regression maximum likelihood estimates. Refer to Table 14.1, b";
title2 "See also page 579";run;
proc print;
title "predicted probabilities. Refer to Table 14.1, page 566";run;
proc sort data=fitted out=fitted;by exper;
proc logistic;
model success(descending)=exper/influence iplots clodds=both clparm=both corrb covb link=cloglog;
output out=fitted predicted=ploglog l=lloglog u=uloglog;
title "logistic regression maximum likelihood estimates. Refer to Table 14.1, b";run;
proc logistic;
model success(descending)=exper/influence iplots clodds=both clparm=both corrb covb link=probit;
output out=fitted predicted=pprobit l=lprobit u=uprobit;
title "logistic regression maximum likelihood estimates. Refer to Table 14.1, b";run;
symbol1 i=none c=red v=circle;
symbol2 i=none c=black v=dot ;
symbol2 i=none c=red v=dot ;
symbol2 i=none c=green v=dot ;
proc gplot;
plot success*exper plogist*exper ploglog*exper pprobit*exper/overlay;
title "observed, logistic, probit and loglog curves. Refer to Figure 14.6, page 568";run;
PROC IMPORT OUT= WORK.coupon 
            DATAFILE= "&location\chapter14.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH14TA02$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title "Coupon Effectiveness data= Refer to Table 14.2, page 569";run;
data coupon;set coupon;
do i=1 to coupons;
redeem=1;
output;
end;
do i=1 to household-coupons;
redeem=0;
output;
end;
drop i predeem;
proc logistic;
model redeem(descending)=discount/influence iplots clodds=both clparm=both corrb covb;
output out=fcoupon predicted=pprob L=lprob U=uprob;
title "Coupon Effectiveness Example. Refer to Table 14.2, page 569";run;
proc sort data=fcoupon out=fcoupon;by discount;
data fcoupon;set fcoupon;by discount;
predeem=coupons/household;
if first.discount;
proc print data=fcoupon;run;
symbol1 i=join v=none;
symbol2 i=non v=dot;
proc gplot;
plot pprob*discount predeem*discount/overlay;
title "Plot of proportions and fitted logistic response. Refer to Figure 14.7, page 570";run;
/*cleaning up fcoupon, given that we do this analysis again later in the chapter*/
proc delete data=fcoupon;run;
PROC IMPORT OUT= WORK.outbreak 
            DATAFILE= "&location\chapter14.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH14TA03$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means data=outbreak;
title "Summary statistics for Disease Outbreak data. Refer to Table 14.3, page 574";run;
proc logistic data=outbreak;
model disease(descending)=age ses1 ses2 sector/influence iplots clodds=both clparm=both corrb covb expb;
output out=fbreakout predicted=pprob L=lprob U=uprob;
title "Disease Outbreak Example. Refer to Table 14.4, page 574";
title2 "note that intercept is not replicated";run;
proc print data=fbreakout;
title "predicted probabilities. Refer also to Table 14.4, page 574";run;
PROC IMPORT OUT= WORK.IPO 
            DATAFILE= "&location\chapter14.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="APPENC11$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title "IPO summary data. Refer to page 575";run;
data ipo;set ipo;
clnvalue=log(value);
proc standard mean=0 data=ipo out=ipo;var clnvalue;
data ipo;set ipo;
clnvalue2=clnvalue*clnvalue;
proc loess data=ipo;
model venture=clnvalue/direct smooth= .4 .5 .6 .7
dfmethod=exact;
ods output outputstatistics=loessout;

data loessout;set loessout;if smoothingparameter=0.5;run;
proc sort data=loessout out=loessout;by smoothingparameter clnvalue;
symbol1 c=black v=dot i=rlclm95;
symbol2 c=red v=point i=join l=1;
symbol3 c=blue v=none i=rlcli95;
proc gplot data=loessout;by smoothingparameter;
plot (depvar pred depvar)*clnvalue/overlay;run;
title 'Compare with Figure 14.9' ;run;

proc logistic data=ipo;
model venture(descending)=clnvalue/influence iplots clodds=both clparm=both corrb covb expb;
output out=fipolin predicted=pprob L=lprob U=uprob;
title "IPO data example. Refer to Table 14.5, page 576";
title2 "Wald Chi-Square is square of z in Table 14.5";run;
proc sort data=loessout out=loessout;by clnvalue;
proc sort data=fipolin out=fipo;by clnvalue;
data fipolin;merge fipolin loessout;by clnvalue;

symbol1 i=join v=none c=red;
symbol2 i=join v=none c=black;
symbol3 i=non v=circle c=blue;
proc gplot data=fipolin;
plot pprob*clnvalue pred*clnvalue venture*clnvalue/overlay;
title "plot of linear and loess regressions. Refer to Figure 14.9a, page 575";run;

proc logistic data=ipo;
model venture(descending)=clnvalue clnvalue2/influence iplots clodds=both clparm=both corrb covb expb;
output out=fipo predicted=pprob L=lprob U=uprob;
title "IPO data example. Refer to Table 14.5, page 576";
title2 "Wald Chi-Square is square of z in Table 14.5";run;
proc sort data=loessout out=loessout;by clnvalue;
proc sort data=fipo out=fipo;by clnvalue;
data fipol;merge fipo loessout;by clnvalue;

symbol1 i=none c=red;
symbol2 i=join v=none c=black;
symbol3 i=none v=circle c=blue;
proc gplot data=fipol;
plot pprob*clnvalue pred*clnvalue venture*clnvalue/overlay;
title "plot of quadratic and loess regressions. Refer to Figure 14.9, page 575";run;
/*Page 579 was done earlier in the programming task example*/
proc logistic data=program;
model success(descending)=exper/influence clodds=both clparm=both corrb covb alpha=.1;
title "Test concerning a single regression weight: Wald test Refer p. 578";
title2 "Square Root of Wald Chi is Wald z test in text.";run;
proc logistic data=program;
model success(descending)=exper/influence clodds=both clparm=both corrb covb alpha=.05;
title "Interval estimation of a single regression weight: Wald confidence intervals Refer p. 579";run;

proc logistic data=program exactoptions(method=direct);
model success(descending)=exper/influence iplots clodds=both clparm=both corrb covb;
exact 'b1' exper/estimate=both;
title "logistic regression maximum likelihood exact estimates.";
title2 "Refer to page 579";run;
proc logistic data=outbreak;
model disease(descending)=age ses1 ses2 sector/clodds=both clparm=both corrb covb expb;
agetest: test age=0;
title "Disease Outbreak Example. Tests of several regression weights=0, page 581";run;

ods listing close;
ods trace on;
ods output FitStatistics=fitfull ;
proc logistic data=outbreak;
model disease(descending)=age ses1 ses2 sector/clodds=both clparm=both corrb covb expb;
title "Disease Outbreak Example. Tests of several regression weights=0, page 581";run;
run;quit;
ods trace off;

ods listing;
options nolabel;
data fitfull;set fitfull;
if criterion="-2 Log L" then Lfull=InterceptAndCovariates;

ods listing close;
ods trace on;
ods output FitStatistics=fitnoAge ;
proc logistic data=outbreak;
model disease(descending)= ses1 ses2 sector/clodds=both clparm=both corrb covb expb;
title "Disease Outbreak Example. Tests of several regression weights=0, page 581";run;
run;quit;
ods trace off;

ods listing;
data fitnoage;set fitnoage;
if criterion="-2 Log L" then LnoAge=InterceptAndCovariates;

ods listing close;
ods trace on;
ods output FitStatistics=fitnoses ;
proc logistic data=outbreak;
model disease(descending)=age sector/clodds=both clparm=both corrb covb expb;
title "Disease Outbreak Example. Tests of several regression weights=0, page 581";run;
run;quit;
ods trace off;

ods listing;
data fitnoses;set fitnoses;
if criterion="-2 Log L" then LnoSes=InterceptAndCovariates;

ods listing close;
ods trace on;
ods output FitStatistics=fitnosector ;
proc logistic data=outbreak;
model disease(descending)=age ses1 ses2 /clodds=both clparm=both corrb covb expb;
title "Disease Outbreak Example. Tests of several regression weights=0, page 581";run;
run;quit;
ods trace off;

ods listing;
data fitnosector;set fitnosector;
if criterion="-2 Log L" then LnoSector=InterceptAndCovariates;
data fitmodels;merge fitfull fitnoage fitnoses fitnosector;
if criterion="-2 Log L";
proc print data=fitmodels;run;
data _null_;set fitmodels;
file print;
agechi=lnoage-lfull;
seschi=lnoses-lfull;
sectorchi=lnosector-lfull;
pagechi=1-cdf('CHISQUARE',agechi,1);
pseschi=1-cdf('CHISQUARE',seschi,2);
psectorchi=1-cdf('CHISQUARE',sectorchi,1);
put "-2 Log Likelihoods for" @30 "Full Model " Lfull/
@30 "Model without Age " LnoAge 11.2 @60 "LDiff=" agechi 11.2 @80 "p(0)=" pagechi 8.6/
@30 "Model without SES " LnoSES 11.2 @60 "LDiff=" seschi 11.2 @80 "p(0)=" pseschi 8.6/
@30 "Model without Sector " Lnosector 11.2 @60 "LDiff=" sectorchi 11.2 @80 "p(0)=" psectorchi 8.6;run;
data outbreak;set outbreak;
ses=3-1*ses1-2*ses2;

ods listing close;
ods trace on;
ods output FitStatistics=fitinter ;
proc logistic data=outbreak;class ses;
model disease(descending)=age SES sector age*ses age*sector ses*sector/clodds=both clparm=both corrb covb expb;
title "Disease Outbreak Example. Model with Interactions";run;
run;quit;
ods trace off;
data fitinter;set fitinter;
if criterion="-2 Log L" then Linter=InterceptAndCovariates;
proc print data=fitinter;run;
data inttest;merge fitfull fitinter;
if criterion="-2 Log L";
file print;
interchi=lfull-linter;
pinterchi=1-cdf('CHISQUARE',interchi,5);
put "-2 Log Likelihoods for" @30 "Interaction Model " Linter/
@30 "Model without Interactions " Lfull 11.2 @60 "LDiff=" interchi 11.2 @80 "p(0)=" pinterchi 8.6;run;
/*The Wald test mentioned in the comment on page 582 is another way to accomplish this 
comparison. To do this, though, you'll need to make the interaction vectors and then use the test statement. 
  This is probably a little easier for many folks, so I include it here:*/
data outbreak;set outbreak;
ageses1=age*ses1;
ageses2=age*ses2;
agesect=age*sector;
ses1sect=ses1*sector;
ses2sect=ses2*sector;
proc logistic data=outbreak;
model disease(descending)=age SES1 ses2 sector ageses1 ageses2 agesect
      ses1sect ses2sect/clodds=both clparm=both corrb covb expb;
testint: test ageses1=ageses2=agesect=ses1sect=ses2sect=0;
title "Disease Outbreak Example. test statement to test interactions to calculate Wald test";
title2 "refer to page 582";run;
proc logistic data=outbreak;
model disease(descending)=age SES1 ses2 sector /clodds=both clparm=both corrb covb expb;
testage: test age=0;
title "Disease Outbreak Example. test statement to test age main effect to calculate Wald test";
title2 "refer to page 582";run;

proc logistic data=coupon;
model redeem(descending)=discount/influence iplots clodds=both clparm=both corrb covb;
output out=fcoupon predicted=predeem;
title "Coupon Effectiveness Example. Refer to Table 14.2, page 569";run;

data fcoupon;set fcoupon;
group=predeem;
proc rank data=fcoupon out=fcoupon groups=5;var group;
data fcoupon;set fcoupon;
unredeem=1-redeem;
punredeem=1-predeem;
proc means data=fcoupon sum mean;class group;
output out=cells sum=spredeem spunred sredeem sunredeem;
var predeem punredeem redeem unredeem;
title "Goodness of fit test. Refer Table14.7, page 587";run;
data cells;set cells;
chigrp=((sredeem-spredeem)**2./spredeem)+((sunredeem-spunred)**2./spunred);
proc print;
proc means data=cells sum;var chigrp;
output out=chi sum=chisq;
data chi;set chi;
file print;
pchisq=1-probchi(chisq,3);
put "Pearson Chi Square test is " chisq 12.4 " with probability " pchisq 8.4;run;
/*Although the above code is of interest if you want to see "where the numbers come from"
  you can simply create the group variable as done above using proc rank to create
  the group variable and rerun the logistic, specifying the aggregate option*/
proc logistic data=Fcoupon;
model redeem(descending)=discount/
 influence iplots clodds=both clparm=both corrb covb aggregate=(group) scale=deviance lackfit;
title "Coupon Effectiveness Example. Goodness of Fit tests, Refer to page 588";run;
data fbreakout;set fbreakout;
group=pprob;
proc rank groups=5 data=fbreakout out=fbreakout;var group;
data fbreakout;set fbreakout;
nodisease=1-disease;
pnodisease=1-pprob;
proc means min max n sum mean;class group;
var disease pprob nodisease pnodisease;
title "observed and expected values for Hosmer-Lemoshow Goodness of Fit";
title2 "refer to Table 14.8, page 590";run;
/*Alternatively, you can specify the Hosmer-Lemeshow test in logistic.
  Note that the SAS procedure splits the data into ten groups if possible, so that 
  is why you don't get exactly the values reported in the book on page 590"*/
proc logistic data=fbreakout;class group;
model disease(descending)=age ses1 ses2 sector/
aggregate=(group) scale=none lackfit rsq clodds=both clparm=both corrb covb expb;
output out=fbreakout predicted=pprob L=lprob U=uprob;
title "Disease Outbreak Example. Hosmer-Lemoshow test. Refer to page 590";
title2 "(Printout also includes r-square fit not discussed in text)";
/*If you want the solution for five groups, you'll need to calculate the Chi square
  yourself as follows:*/
proc means data=fbreakout sum mean;class group;
output out=cells sum=spdisease spnodisease sdisease snodisease;
var pprob pnodisease disease nodisease;
title "Goodness of fit test. Refer Table14.7, page 587";run;
data cells;set cells;
if group=. then delete;
chigrp=((sdisease-spdisease)**2./spdisease)+((snodisease-spnodisease)**2./spnodisease);
proc print;
proc means data=cells sum;var chigrp;
output out=chi sum=chisq;
data chi;set chi;
file print;
pchisq=1-probchi(chisq,3);
put "Hosmer-Lemeshow GoF test is " chisq 12.4 " with probability " pchisq 8.4;
title "Note that value is not replicated. Book reports Chi of 1.98";
title2 "Even though expected and observed frequencies from data set match";run;
proc logistic data=outbreak;
model disease(descending)=age ses1 ses2 sector/
influence clodds=both clparm=both corrb covb expb;
output out=fbreakout predicted=pprob reschi=pearson resdev=dev h=h;
run;
data fbreakout;set fbreakout;
linpred=log(pprob/(1-pprob));
resid=disease-pprob;
stresid=pearson/sqrt(1-h);
proc print data=fbreakout;var disease pprob resid pearson stresid dev h;
title "Logistic residuals and hat matrix. Rfer to Table 14.0, page 593";run;
symbol1 v=circle;
proc gplot;
plot resid*pprob;
title "Plot of error by predicted probability. Refer to Figure 14.12a.";run;
proc gplot;
plot pearson*pprob;
title "Plot of pearson residual by predicted probability. Refer to Figure 14.12b.";run;
proc gplot;
plot stresid*pprob;
title "Plot of studentized pearson residual by predicted probability. Refer to Figure 14.12c.";run;
proc gplot;
plot dev*pprob;
title "Plot of deviance by predicted probability. Refer to Figure 14.12d.";run;
proc sort data=fbreakout out=fbreakout;by pprob;*Data should be sorted by predictor first;
proc loess data=fbreakout;
model stresid=pprob/direct smooth=.7
dfmethod=exact;
ods output outputstatistics=loessout;
proc print data=loessout;run;
proc sort data=loessout out=loessout;by smoothingparameter pprob;
symbol1 c=black v=dot i=rlclm95;
symbol2 c=red v=point i=join l=1;
symbol3 c=blue v=none i=rlcli95;
proc gplot data=loessout;by smoothingparameter;
plot (depvar pred depvar)*pprob/overlay;run;
title 'Lowess smoothing of studentized residuals by predicted probability. Refer to Figure 14.13a page 595' ;run;

proc loess data=fbreakout;
model stresid=linpred/direct smooth=.7
dfmethod=exact;
ods output outputstatistics=loessout;
proc print data=loessout;run;
proc sort data=loessout out=loessout;by smoothingparameter linpred;
symbol1 c=black v=dot i=rlclm95;
symbol2 c=red v=point i=join l=1;
symbol3 c=blue v=none i=rlcli95;
proc gplot data=loessout;by smoothingparameter;
plot (depvar pred depvar)*linpred/overlay;run;
title 'Lowess smoothing of studentized residuals by linear predictor. Refer to Figure 14.13b page 595' ;run;

proc loess data=fbreakout;
model dev=pprob/direct smooth=.7
dfmethod=exact;
ods output outputstatistics=loessout;
proc print data=loessout;run;
proc sort data=loessout out=loessout;by smoothingparameter pprob;
symbol1 c=black v=dot i=rlclm95;
symbol2 c=red v=point i=join l=1;
symbol3 c=blue v=none i=rlcli95;
proc gplot data=loessout;by smoothingparameter;
plot (depvar pred depvar)*pprob/overlay;run;
title 'Lowess smoothing of deviance by predicted probability. Refer to Figure 14.13c page 595' ;run;

proc loess data=fbreakout;
model dev=linpred/direct smooth=.7
dfmethod=exact;
ods output outputstatistics=loessout;
proc print data=loessout;run;
proc sort data=loessout out=loessout;by smoothingparameter linpred;
symbol1 c=black v=dot i=rlclm95;
symbol2 c=red v=point i=join l=1;
symbol3 c=blue v=none i=rlcli95;
proc gplot data=loessout;by smoothingparameter;
plot (depvar pred depvar)*linpred/overlay;run;
title 'Lowess smoothing of deviance by linear predictor. Refer to Figure 14.13d page 595' ;run;
%include "&location\halfnorm.sas";
%halfnorm(fbreakout,age ses1 ses2 sector,pprob,dev,98);
proc delete data=fbreakout;
proc logistic data=outbreak;
model disease(descending)=age ses1 ses2 sector/
influence clodds=both clparm=both corrb covb expb;
output out=fbreakout predicted=pprob reschi=pearson resdev=dev h=h difdev=deltadev difchisq=deltchisq ;
data fbreakout;set fbreakout;
rsp=pearson/(sqrt(1-h));
cookd=(h*pearson**2.)/(5*(1-h)**2.);*here in denominator number of parameters is 4+1=5;
proc print;
title "Pearson residuals, studentized residuals, hat, deviance, delta chi, delta deviance and cooks's D";
title2 "refer to Table 14.11, page 599";
run;
proc gplot;
plot (deltadev deltchisq)*case;
title "Delta Chi Square and Deviance plots, refer to figure 14.15 a and b, page 600";run;
proc sort data=fbreakout out=fbreakout;by pprob;
symbol1 i=none v=circle c=red;
proc gplot;
plot (deltadev deltchisq)*pprob;
title "Delta Chi Square and Deviance Plots, refer to figure 14.15 c and d, page 600";run;
symbol1 i=join v=circle c=red;
proc sort data=fbreakout out=fbreakout;by case;
proc gplot;
plot (h cookd)*case;
title "Lever and Cook's d plots, refer to Figure 14.16 a and b, page 601";run;
symbol1 i=none v=circle c=red;

proc sort data=fbreakout out=fbreakout;by pprob;
symbol1 i=none v=circle c=red;
proc gplot;
plot deltadev*pprob;
title "delta deviances by predicted probability, Refer to Figrue 14.16c, page 601";run;
proc delete data=fbreakout;
data newobs;
age=10;ses1=0;ses2=1;sector=0;
data outbreak2;set newobs outbreak;
proc logistic data=outbreak2;
model disease(descending)=age ses1 ses2 sector/outroc=roccurve
influence clodds=both clparm=both corrb covb expb;
output out=fbreakout predicted=pprob l=lower u=upper;
data fbreakout;set fbreakout;
if disease=.;
proc print data=fbreakout;
title "Estimated confidence intervals on predicted probabilities. Refer to pages 602-603";run;
proc print data=roccurve;
title "false positive and negative weights by different cutoffs. Refer to Table 14.12, page 605";
run;
symbol1 i=join c=red v=circle;
proc sort data=roccurve out=roccurve;by _1mspec_;
proc gplot;
plot _sensit_*_1mspec_;
title "Receiver Operating Curve. Refer to Figure 14.17, page 606";run;
PROC IMPORT OUT= WORK.valid 
            DATAFILE= "&location\chapter14.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="APPENC10$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc delete data=fbreakout;
data new;set valid;if case>=99;
if ses=2 then ses1=1;else ses1=0;
if ses=3 then ses2=1;else ses2=0;
proc delete data=outval validroc fbreakout;
proc logistic data=outbreak2;
model disease=age ses1 ses2 sector/;
output out=fbreakout predicted=pprob l=lower u=upper;
score data=new out=outval clm fitstat clm outroc=validroc;
title "validation of prediction error rate. Refer to p. 607-608";run;
data validroc;set validroc;
mindist=abs(_prob_-.3);
proc sort data=validroc out=validroc;by mindist;
data validroc;set validroc;if _n_=1;
diserr=100*_falpos_/(_neg_+_falpos_);
nodiserr=100*_falneg_/(_pos_+_falneg_);
toterr=100*(_falpos_+_falneg_)/(_neg_+_falpos_+_pos_+_falneg_);
file print;
put "Percent prediction error rates "/
 "with Disease " diserr 8.2/
 "without Disease " nodiserr 8.2/
 "Total " toterr 8.2;

proc print data=outval;run;
proc print data=validroc;run;
proc freq data=outval;tables disease*(F_disease i_disease);run;
PROC IMPORT OUT= WORK.PregDur 
            DATAFILE= "&location\chapter14.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH14TA13$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc print;
title "pregnancy duration data. Refer to Table 14.13, page 609";run;
proc logistic data=pregdur;class term(ref=last);
model term(order=data)=nutrit age1 age2 alcohol smoking/
influence clodds=both clparm=both corrb covb expb link=glogit;
title "Nominal logistic regression. Refer to Figure 14.18, page 613";run;
proc logistic data=pregdur;class term;
model term=nutrit age1 age2 alcohol smoking/
influence clodds=both clparm=both corrb covb expb;
title "ordinal logistic regression. Refer to Figure 14.19, page 618";run;
PROC IMPORT OUT= WORK.Lumber 
            DATAFILE= "&location\chapter14.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH14TA14$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc print data=lumber;
title "Miller Lumber Company example. Refer to Table 14.14, page 622";run;
proc genmod data=lumber;
model customer=housing income age competit distance/dist=poisson link=log;
output out=plumber p=predcnt resdev=resdev;
title "Poisson regression. Refer Table 14.15, page 622";
data plumber;set plumber;
i+1;
proc print;var customer predcnt resdev;
title "Responses, fitted values, and deviance residuals. Refer to Table 14.16, page 623";
proc gplot;
plot resdev*i;
title "Index plot of deviance residuals. Refer to Figure 14.20, page 623";run;