/* Use the following DATA step to obtain Outputs 10.24 through 10.27  */

Data Counts;
  input block trt a b count;
  ctl_trt=(trt>0);
datalines;
1 1 1 1 6
1 2 1 2 2
1 5 2 2 3
1 8 3 2 3
1 7 3 1 1
1 0 0 0 16
1 3 1 3 4
1 6 2 3 1
1 9 3 3 1
1 4 2 1 5
2 1 1 1 9
2 2 1 2 6
2 5 2 2 4
2 8 3 2 2
2 7 3 1 2
2 0 0 0 25
2 3 1 3 3
2 6 2 3 5
2 9 3 3 0
2 4 2 1 3
3 1 1 1 2
3 2 1 2 14
3 5 2 2 6
3 8 3 2 3
3 7 3 1 2
3 0 0 0 5
3 3 1 3 5
3 6 2 3 17
3 9 3 3 2
3 4 2 1 3
4 1 1 1 22
4 2 1 2 4
4 5 2 2 3
4 8 3 2 4
4 7 3 1 3
4 0 0 0 9
4 3 1 3 5
4 6 2 3 1
4 9 3 3 9
4 4 2 1 2
;

proc sort; 
 by block trt;
proc print;
 var block trt ctl_trt a b count;
run;

proc genmod data=counts;
 class BLOCK CTL_TRT a b;
 model count=BLOCK CTL_TRT a b a*b/dist=negbin type1 type3 wald;
 title 'uncorrected Poisson model';
run;

proc genmod data=counts;
 class BLOCK CTL_TRT a b;
 model count=BLOCK CTL_TRT a b a*b/dist=poisson type1 ObStats;
 ODS OUTPUT ObStats=check;
 title 'compute model checking statistics';
run;

data plot;
 merge check;
 adjlamda=2*sqrt(pred);
 ystar=xbeta+(count-pred)/pred;
 absres=abs(resdev);

proc plot;
 plot resdev*(pred xbeta);
 plot (resdev reschi)*adjlamda;
 plot ystar*xbeta;
 plot absres*adjlamda;
run;

/* Use the following SAS statements in conjunction with     */
/* data set DATA=Counts (shown in Output 10.24)             */
/*    to obtain Output 10.28                                */

proc genmod data=Counts;
 class BLOCK CTL_TRT a b;
 model count=BLOCK CTL_TRT a b a*b/dist=poisson dscale type1 type3;
 title ' correction for overdispersion';
run;

/* Use the following SAS statements in conjunction with     */
/* data set DATA=Counts (shown in Output 10.24)             */
/*    to obtain Output 10.29                                */

proc genmod data=Counts;
 class block trt;
 model count=block trt/dist=poisson type1 type3 wald;
 contrast 'ctl vs trt' trt 9 -1 -1 -1 -1 -1 -1 -1 -1 -1/wald;
 contrast 'a' trt 0  1 1 1  0 0 0  -1 -1 -1,
              trt 0  0 0 0  1 1 1  -1 -1 -1/wald;
 contrast 'b' trt 0  1 0 -1  1 0 -1  1 0 -1,
              trt 0  0 1 -1  0 1 -1  0 1 -1/wald;
 contrast 'a x b' trt 0  1 0 -1  0 0  0  -1 0 1,
                  trt 0  0 0  0  1 0 -1  -1 0 1,
                                  trt 0  0 1 -1  0 0  0  0 -1 1,
                                  trt 0  0 0  0  0 1 -1  0 -1 1/wald;
 estimate 'ctl lsmean' intercept 1 trt 1 0/exp;
 estimate 'treated lsm'intercept 9 trt 0 1 1 1 1 1 1 1 1 1/divisor=9 exp;
 estimate 'A=1 lsmean' intercept 3 trt 0 1 1 1 0/divisor=3 exp;
 estimate 'A=2 lsmean' intercept 3 trt 0 0 0 0 1 1 1 0/divisor=3 exp;
 estimate 'A=3 lsmean' intercept 3 trt 0 0 0 0 0 0 0 1 1 1 0/divisor=3 exp;
run; 

/* Use the following SAS statements in conjunction with     */
/* data set DATA=Counts (shown in Output 10.24)             */
/*    to obtain Outputs 10.30 and 10.31                     */

proc genmod data=Counts;
 class BLOCK CTL_TRT a b;
 model count=BLOCK CTL_TRT a b a*b/dist=negbin type1 type3;
 title 'uncorrected Poisson model';
run;

proc genmod data=Counts;
 class BLOCK TRT;
 model count=BLOCK TRT/dist=negbin type1 type3 Wald;
 contrast 'ctl vs trt' trt 9 -1 -1 -1 -1 -1 -1 -1 -1 -1/wald;
 contrast 'a' trt 0  1 1 1  0 0 0  -1 -1 -1,
              trt 0  0 0 0  1 1 1  -1 -1 -1/wald;
 contrast 'b' trt 0  1 0 -1  1 0 -1  1 0 -1,
              trt 0  0 1 -1  0 1 -1  0 1 -1/wald;
 contrast 'a x b' trt 0  1 0 -1  0 0  0  -1 0 1,
                  trt 0  0 0  0  1 0 -1  -1 0 1,
                                  trt 0  0 1 -1  0 0  0  0 -1 1,
                                  trt 0  0 0  0  0 1 -1  0 -1 1/wald;
 estimate 'ctl lsmean' intercept 1 trt 1 0/exp;
 estimate 'treated lsm'intercept 9 trt 0 1 1 1 1 1 1 1 1 1/divisor=9 exp;
 estimate 'A=1 lsmean' intercept 3 trt 0 1 1 1 0/divisor=3 exp;
 estimate 'A=2 lsmean' intercept 3 trt 0 0 0 0 1 1 1 0/divisor=3 exp;
 estimate 'A=3 lsmean' intercept 3 trt 0 0 0 0 0 0 0 1 1 1 0/divisor=3 exp;
 title 'uncorrected Poisson model';
run;

/* Use the following SAS statements in conjunction with     */
/* data set DATA=Counts (shown in Output 10.24)             */
/*    to obtain Outputs 10.32 through 10.38                 */

PROC GENMOD data=Counts;
 k=1/0.2383;
 mu=_mean_;
 eta=_xbeta_;
 fwdlink link=log(mu/(mu+k));
 invlink ilink=k*exp(eta)/(1-exp(eta));
 CLASS BLOCK TRT;
 MODEL count=BLOCK TRT/ dist=negbin type1 type3 wald;
 contrast 'ctl vs trt' trt 9 -1 -1 -1 -1 -1 -1 -1 -1 -1;
 contrast 'a' trt 0  1 1 1  0 0 0  -1 -1 -1,
              trt 0  0 0 0  1 1 1  -1 -1 -1;
 contrast 'b' trt 0  1 0 -1  1 0 -1  1 0 -1,
              trt 0  0 1 -1  0 1 -1  0 1 -1;
 contrast 'a x b' trt 0  1 0 -1  0 0  0  -1 0 1,
                  trt 0  0 0  0  1 0 -1  -1 0 1,
                                  trt 0  0 1 -1  0 0  0  0 -1 1,
                                  trt 0  0 0  0  0 1 -1  0 -1 1;
 title 'negative binomial model';
run;

PROC GENMOD data=Counts;
 k=1/0.2383;
 mu=_mean_;
 eta=_xbeta_;
 fwdlink link=log(mu/(mu+k));
 invlink ilink=k*exp(eta)/(1-exp(eta));
 CLASS BLOCK CTL_TRT a b;
 *MODEL count=BLOCK CTL_TRT/ type3 wald;
 MODEL count=BLOCK CTL_TRT a b a*b/type3 wald ;
run;

PROC GENMOD data=Counts;
 count=_resp_;
 y=count;
 if y=0 then y=0.1;
 mu=_mean_;
 eta=_xbeta_;
 K=1;
 fwdlink link=log(mu/(mu+k));
 invlink ilink=k*exp(eta)/(1-exp(eta));
 lly=y*log(y/(y+k))-k*log((k+y)/k);
 llm=y*log(mu/(mu+k))-k*log((k+mu)/k);
 d=2*(lly-llm);
 variance var=mu+(mu*mu/k);
 deviance dev=d;
 class block ctl_trt a b;
 *model y=block ctl_trt/type3 wald;
 model y=block ctl_trt a b a*b/type3 wald;

PROC GENMOD data=Counts;
 K=2.5;
 count=_resp_;
 y=count;
 if y=0 then y=0.1;
 mu=_mean_;
 eta=_xbeta_;
 fwdlink link=log(mu/(mu+k));
 invlink ilink=k*exp(eta)/(1-exp(eta));
 lly=y*log(y/(y+k))-k*log((k+y)/k);
 llm=y*log(mu/(mu+k))-k*log((k+mu)/k);
 d=2*(lly-llm);
 variance var=mu+(mu*mu/k);
 deviance dev=d;
 CLASS BLOCK TRT;
 MODEL y=trt block;
 contrast 'ctl vs trt' trt 9 -1 -1 -1 -1 -1 -1 -1 -1 -1/wald;
 contrast 'a' trt 0  1 1 1  0 0 0  -1 -1 -1,
              trt 0  0 0 0  1 1 1  -1 -1 -1/wald;
 contrast 'b' trt 0  1 0 -1  1 0 -1  1 0 -1,
              trt 0  0 1 -1  0 1 -1  0 1 -1/wald;
 contrast 'a x b' trt 0  1 0 -1  0 0  0  -1 0 1,
                  trt 0  0 0  0  1 0 -1  -1 0 1,
                                  trt 0  0 1 -1  0 0  0  0 -1 1,
                                  trt 0  0 0  0  0 1 -1  0 -1 1/wald;
 contrast 'ctl vs trt' trt 9 -1 -1 -1 -1 -1 -1 -1 -1 -1;
 contrast 'a' trt 0  1 1 1  0 0 0  -1 -1 -1,
              trt 0  0 0 0  1 1 1  -1 -1 -1;
 contrast 'b' trt 0  1 0 -1  1 0 -1  1 0 -1,
              trt 0  0 1 -1  0 1 -1  0 1 -1;
 contrast 'a x b' trt 0  1 0 -1  0 0  0  -1 0 1,
                  trt 0  0 0  0  1 0 -1  -1 0 1,
                                  trt 0  0 1 -1  0 0  0  0 -1 1,
                                  trt 0  0 0  0  0 1 -1  0 -1 1;
 estimate 'ctl lsmean' intercept 1 trt 1 0;
 estimate 'treated lsm'intercept 9 trt 0 1 1 1 1 1 1 1 1 1/divisor=9;
 estimate 'A=1 lsmean' intercept 3 trt 0 1 1 1 0/divisor=3;
 estimate 'A=2 lsmean' intercept 3 trt 0 0 0 0 1 1 1 0/divisor=3;
 estimate 'A=3 lsmean' intercept 3 trt 0 0 0 0 0 0 0 1 1 1 0/divisor=3;
 ods output estimates=lsm;
 title 'negative binomial model';
run;

data c_lsm;
 set lsm;
 *k=1;
 k=2.5;
 counthat=k*exp(estimate)/(1-exp(estimate));
 deriv=k*exp(estimate)/((1-exp(estimate))**2);
 se_count=deriv*stderr;

proc print data=c_lsm;
run;