/***********************************************************/
/*Data from Stokes, Davis, and Koch (1995).  Dose is a     */
/*bacteria count, and response is whether animal survives. */
/*File biostat1macros.sas has the doci macro from last     */
/*class.  Data are at                                      */
/* support.sas.com/rnd/app/da/cat/samples/chapter11.html   */
/***********************************************************/
options ls=64;
data out; run;
%include 'biostat1macros.sas';
data assay; infile 'stokes.dat'; 
   input drug $ dose status $ count; ldose=log(dose); run;
proc genmod data=assay descending; class status drug; 
   freq count;
   model status=drug ldose ldose*ldose ldose*drug 
      ldose*ldose*drug/dist=bin; run;
title1 'Small Model';
proc genmod data=assay descending; class status drug; 
   freq count; model status=drug ldose/dist=bin corrb ; run;
title1 'CI for LD50 for standard preparation';
/* Slope coefficient below is negated to make barY */
data coef; input barW barY m n tau sigma df rho; cards;
-2.4476 -1.2545 1 1 0.4532 0.5033 .  0.6837
;
run;
%doci(coef,.05,'LD50');
title1 'CI for LD75 for standard preparation';
data coef; set coef; barY=barY+log(.75)-log(.25); run;
%doci(coef,.05,'LD75');
title1 'CI for potency';
data coef; input barW barY m n tau sigma df rho; cards;
0.7234 1.2545 1 1 0.1177 0.5033 .  0.4496
;
run;
%doci(coef,.05,'Potency');
proc print data=out noobs; 
   var runn name barW barY tau sigma rho alpha bot 
   elower eupper; run;
/* Repeat for Probit regression */
data out; run;
proc genmod data=assay descending; class status drug; 
   freq count; 
   model status=drug ldose/dist=bin link=probit corrb ; run;
title1 'CI for LD50 for standard preparation';
/* Slope coefficient below is negated to make barY; */
/* rho also neg.*/
data coef; input barW barY m n tau sigma df rho; cards;
 -1.4164 -0.4096 1 1 0.2480 0.0614 .  0.6133
;
run;
%doci(coef,.05,'LD50');
title1 'CI for LD75 for standard preparation';
data coef; set coef;
barY=barY+probit(.75);
run;
%doci(coef,.05,'LD75');
title1 'CI for potency';
data coef; input barW barY m n tau sigma df rho; cards;
0.4096 0.7682 1 1 0.0614 0.2911 .  .4611
;
run;
%doci(coef,.05,'Potency');
proc print data=out noobs; 
   var runn name barW barY tau sigma rho alpha bot 
      elower eupper; run;