data Trajan;
   input roots shoot photoperiod bap;
   lshoot=log(shoot);
   datalines;
0 40 8 17.6
0 40 8 17.6
0 30 16 2.2
0 30 16 2.2
0 30 16 2.2
0 30 16 2.2
0 30 16 2.2
0 30 16 2.2
0 30 16 2.2
0 30 16 2.2
0 30 16 2.2
0 30 16 2.2
0 30 16 2.2
0 30 16 2.2
0 30 16 2.2
0 30 16 2.2
0 30 16 2.2
0 30 16 4.4
0 30 16 4.4
0 30 16 4.4
0 30 16 4.4
0 30 16 4.4
0 30 16 4.4
0 30 16 4.4
0 30 16 4.4
0 30 16 4.4
0 30 16 4.4
0 30 16 4.4
0 30 16 4.4
0 30 16 4.4
0 30 16 4.4
0 30 16 4.4
0 30 16 4.4
0 30 16 8.8
0 30 16 8.8
0 30 16 8.8
0 30 16 8.8
0 30 16 8.8
0 30 16 8.8
0 30 16 8.8
0 30 16 8.8
0 30 16 8.8
0 30 16 8.8
0 30 16 8.8
0 30 16 8.8
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
0 40 16 17.6
1 30 8 2.2
1 30 8 2.2
1 30 8 2.2
1 30 16 4.4
1 30 16 4.4
1 30 16 8.8
1 30 16 8.8
1 30 16 8.8
1 40 16 17.6
1 40 16 17.6
2 30 8 2.2
2 30 8 2.2
2 30 8 4.4
2 30 8 4.4
2 30 8 4.4
2 40 8 8.8
2 30 16 2.2
2 30 16 2.2
2 30 16 4.4
2 30 16 8.8
2 30 16 8.8
2 40 16 17.6
2 40 16 17.6
3 30 8 2.2
3 30 8 2.2
3 30 8 2.2
3 40 8 8.8
3 40 8 8.8
3 40 8 17.6
3 40 8 17.6
3 30 16 2.2
3 30 16 2.2
3 30 16 4.4
3 30 16 8.8
3 40 16 17.6
3 40 16 17.6
3 40 16 17.6
3 40 16 17.6
4 30 8 2.2
4 30 8 2.2
4 30 8 2.2
4 30 8 2.2
4 30 8 2.2
4 30 8 2.2
4 30 8 4.4
4 40 8 8.8
4 40 8 8.8
4 40 8 8.8
4 40 8 8.8
4 40 8 17.6
4 40 8 17.6
4 30 16 2.2
4 30 16 4.4
4 30 16 4.4
4 30 16 8.8
4 30 16 8.8
4 40 16 17.6
4 40 16 17.6
4 40 16 17.6
5 30 8 2.2
5 30 8 2.2
5 30 8 2.2
5 40 8 8.8
5 40 8 8.8
5 40 8 8.8
5 40 8 8.8
5 40 8 17.6
5 40 8 17.6
5 40 8 17.6
5 40 8 17.6
5 40 8 17.6
5 30 16 2.2
5 30 16 2.2
5 30 16 4.4
5 30 16 8.8
5 30 16 8.8
5 40 16 17.6
6 30 8 2.2
6 30 8 2.2
6 30 8 4.4
6 30 8 4.4
6 30 8 4.4
6 40 8 8.8
6 40 8 8.8
6 40 8 8.8
6 40 8 8.8
6 40 8 17.6
6 40 8 17.6
6 40 8 17.6
6 40 8 17.6
6 40 8 17.6
6 30 16 2.2
6 30 16 4.4
6 30 16 4.4
6 30 16 8.8
6 30 16 8.8
6 30 16 8.8
6 40 16 17.6
6 40 16 17.6
6 40 16 17.6
6 40 16 17.6
7 30 8 2.2
7 30 8 2.2
7 30 8 4.4
7 30 8 4.4
7 30 8 4.4
7 30 8 4.4
7 30 8 4.4
7 30 8 4.4
7 30 8 4.4
7 40 8 8.8
7 40 8 8.8
7 40 8 8.8
7 40 8 8.8
7 40 8 17.6
7 40 8 17.6
7 40 8 17.6
7 40 8 17.6
7 30 16 8.8
7 40 16 17.6
7 40 16 17.6
7 40 16 17.6
8 30 8 2.2
8 30 8 2.2
8 30 8 2.2
8 30 8 4.4
8 30 8 4.4
8 30 8 4.4
8 40 8 8.8
8 40 8 8.8
8 40 8 8.8
8 40 8 8.8
8 40 8 8.8
8 40 8 8.8
8 40 8 8.8
8 40 8 17.6
8 40 8 17.6
8 40 8 17.6
8 40 8 17.6
8 40 8 17.6
8 40 8 17.6
8 40 8 17.6
8 40 8 17.6
8 30 16 2.2
8 30 16 4.4
9 30 8 2.2
9 30 8 4.4
9 30 8 4.4
9 30 8 4.4
9 30 8 4.4
9 30 8 4.4
9 40 8 8.8
9 40 8 8.8
9 40 8 8.8
9 40 8 8.8
9 40 8 8.8
9 40 8 17.6
9 40 8 17.6
9 40 8 17.6
9 30 16 2.2
9 30 16 2.2
9 30 16 2.2
9 30 16 8.8
9 30 16 8.8
9 40 16 17.6
9 40 16 17.6
10 30 8 2.2
10 30 8 2.2
10 30 8 4.4
10 30 8 4.4
10 30 8 4.4
10 40 8 8.8
10 40 8 8.8
10 40 8 8.8
10 40 8 8.8
10 40 8 17.6
10 40 8 17.6
10 40 8 17.6
10 40 8 17.6
10 30 16 2.2
10 30 16 4.4
10 30 16 4.4
10 30 16 4.4
11 30 8 2.2
11 30 8 4.4
11 30 8 4.4
11 30 8 4.4
11 30 8 4.4
11 40 8 8.8
11 40 8 17.6
11 40 8 17.6
11 40 8 17.6
11 40 8 17.6
11 30 16 2.2
11 30 16 8.8
12 40 8 8.8
12 40 8 8.8
12 30 16 2.2
12 30 16 4.4
12 30 16 8.8
13 30 8 2.2
13 30 8 4.4
14 40 8 8.8
14 40 8 8.8
14 40 8 17.6
17 30 8 2.2
;


ods graphics on;
proc freq data=Trajan;
   table roots / plots(only)=freqplot(scale=percent);
run;


proc sort data=Trajan out=Trajan;
   by photoperiod;
run;


proc freq data=Trajan;
   table roots / plots(only)=freqplot(scale=percent);
   by photoperiod;
run;


proc sort data=Trajan out=Trajan;
   by bap;
run;


proc freq data=Trajan;
   table roots / plots(only)=freqplot(scale=percent);
   by bap;
run;


proc sort data=Trajan out=Trajan;
   by photoperiod bap;
run;


proc freq data=Trajan;
   table roots / plots(only)=freqplot(scale=percent);
   by photoperiod bap;
run;


proc genmod data=Trajan;
   class bap photoperiod;
   model roots = bap|photoperiod / dist=zip offset=lshoot;
   zeromodel photoperiod;
   output out=zip predicted=pred pzero=pzero;
   ods output Modelfit=fit;
run;


data fit;
   set fit(where=(criterion="Scaled Pearson X2"));
   format pvalue pvalue6.4;
   pvalue=1-probchi(value,df);
run;


proc print data=fit noobs;
   var criterion value df pvalue;
run;


proc means data=Trajan noprint;
   var roots;
   output out=maxcount max=max N=N;
run;


data _null_;
  set maxcount;
  call symput('N',N);
  call symput('max',max);
run;


%let max=%sysfunc(strip(&max));


data zip(drop= i);
   set zip;
   lambda=pred/(1-pzero);
   array ep{0:&max} ep0-ep&max;
   array c{0:&max} c0-c&max;
   do i = 0 to &max;
      if i=0 then ep{i}= pzero + (1-pzero)*pdf('POISSON',i,lambda);
      else        ep{i}=         (1-pzero)*pdf('POISSON',i,lambda);
      c{i}=ifn(roots=i,1,0);
   end;
run;


proc means data=zip noprint;
   var ep0 - ep&max c0-c&max;
   output out=ep(drop=_TYPE_ _FREQ_) mean(ep0-ep&max)=ep0-ep&max;
   output out=p(drop=_TYPE_ _FREQ_) mean(c0-c&max)=p0-p&max;
run;


proc transpose data=ep out=ep(rename=(col1=zip) drop=_NAME_);
run;

proc transpose data=p out=p(rename=(col1=p) drop=_NAME_);
run;


data zipprob;
   merge ep p;
   zipdiff=p-zip;
   roots=_N_ -1;
   label zip='ZIP Probabilities'
         p='Relative Frequencies'
         zipdiff='Observed minus Predicted';
run;


proc sgplot data=zipprob;
   scatter x=roots y=p /
           markerattrs=(symbol=CircleFilled size=5px color=blue);
   scatter x=roots y=zip /
           markerattrs=(symbol=TriangleFilled size=5px color=red);
   xaxis type=discrete;
run;


proc sgplot data=zipprob;
   series x=roots y=zipdiff /
           lineattrs=(pattern=ShortDash  color=blue)
           markers markerattrs=(symbol=CircleFilled size=5px color=blue);
   refline 0/ axis=y;
   xaxis type=discrete;
run;


proc genmod data=Trajan;
   class bap photoperiod;
   model roots = bap|photoperiod / dist=zinb offset=lshoot;
   zeromodel photoperiod;
   output out=zinb predicted=pred pzero=pzero;
   ods output ParameterEstimates=zinbparms;
   ods output Modelfit=fit;
run;


data fit;
   set fit(where=(criterion="Scaled Pearson X2"));
   format pvalue pvalue6.4;
   pvalue=1-probchi(value,df);
run;


proc print data=fit noobs;
   var criterion value df pvalue;
run;


data zinbparms;
   set zinbparms(where=(Parameter="Dispersion"));
   keep estimate;
   call symput('k',estimate);
run;


data zinb(drop= i);
   set zinb;
   lambda=pred/(1-pzero);
   k=&k;
   array ep{0:&max} ep0-ep&max;
   array c{0:&max} c0-c&max;
   do i = 0 to &max;
      if i=0 then ep{i}= pzero + (1-pzero)*pdf('NEGBINOMIAL',i,(1/(1+k*lambda)),(1/k));
      else        ep{i}=         (1-pzero)*pdf('NEGBINOMIAL',i,(1/(1+k*lambda)),(1/k)); 
      c{i}=ifn(roots=i,1,0);  
   end;
run;


proc means data=zinb noprint;
   var ep0 - ep&max c0-c&max;
   output out=ep(drop=_TYPE_ _FREQ_) mean(ep0-ep&max)=ep0-ep&max;
   output out=p(drop=_TYPE_ _FREQ_) mean(c0-c&max)=p0-p&max;
run;


proc transpose data=ep out=ep(rename=(col1=zinb) drop=_NAME_);
run;

proc transpose data=p out=p(rename=(col1=p) drop=_NAME_);
run;


data zinbprob;
   merge ep p;
   zinbdiff=p-zinb;
   roots=_N_ -1;
   label zinb='ZINB Probabilities'
         p='Relative Frequencies'
         zinbdiff='Observed minus Predicted';
run;


proc sgplot data=zinbprob;
   scatter x=roots y=p /
           markerattrs=(symbol=CircleFilled size=5px color=blue);
   scatter x=roots y=zinb /
           markerattrs=(symbol=TriangleFilled size=5px color=red);
   xaxis type=discrete;
run;


proc sgplot data=zinbprob;
   series x=roots y=zinbdiff /
           lineattrs=(pattern=ShortDash  color=blue)
           markers markerattrs=(symbol=CircleFilled size=5px color=blue);
   refline 0/ axis=y;
   xaxis type=discrete;
run;


data compare;
   merge zipprob zinbprob;
   by roots;
run;


proc sgplot data=compare;
   series x=roots y=zinbdiff /
           lineattrs=(pattern=Solid  color=red)
           markers markerattrs=(symbol=TriangleFilled color=red)
           legendlabel="ZINB";
   series x=roots y=zipdiff /
           lineattrs=(pattern=ShortDash  color=blue)
           markers markerattrs=(symbol=StarFilled color=blue)
           legendlabel="ZIP";
   refline 0/ axis=y;
   xaxis type=discrete;
run;


proc sort data=zinb out=zinb;
   by photoperiod bap;
run;


proc means data=zinb;
   var pred;
   by photoperiod bap;
   output out=effects  mean(pred)=pred;
run;


proc sgpanel data=effects;
   panelby photoperiod;
   series x=bap y=pred / markers
          markerattrs=(symbol=CircleFilled size=5px color=red);
   colaxis type=discrete;
run;