%let nObs = 5000;
%let nClass = 5;
%let nLevs = 4;
%let seed = 1234;

data tmp1;
   array c{&nClass};

   keep c1-c&nClass yTweedie d1 d2;

   /* Tweedie parms */
   phi=0.5;
   p=1.5;

   do i=1 to &nObs;

      do j=1 to &nClass;
         c{j} = int(ranuni(1)*&nLevs);
      end;

      d1 = ranuni(&seed);
      d2 = ranuni(&seed);

      xBeta = 0.5*((c2<2) - 2*(c1=1) + 0.5*c&nClass + 0.05*d1);
      mu = exp(xBeta);

      /* Poisson distributions parms */
      lambda = mu**(2-p)/(phi*(2-p));
      /* Gamma distribution parms */
      alpha = (2-p)/(p-1);
      gamma = phi*(p-1)*(mu**(p-1));

      rpoi = ranpoi(&seed,lambda);
      if rpoi=0 then yTweedie=0;
      else do;
         yTweedie=0;
         do j=1 to rpoi;
         yTweedie = yTweedie + rangam(&seed,alpha);
         end;
         yTweedie = yTweedie * gamma;
      end;
      output;
   end;
run;

11   ods listing close;ods html5 file=stdout options(bitmap_mode='inline') device=png; ods graphics on / outputfmt=png;
NOTE: Writing HTML5 Body file: STDOUT
12   
13   %let nObs = 5000;
14   %let nClass = 5;
15   %let nLevs = 4;
16   %let seed = 1234;
17   
18   data tmp1;
19      array c{&nClass};
20   
21      keep c1-c&nClass yTweedie d1 d2;
22   
23      /* Tweedie parms */
24      phi=0.5;
25      p=1.5;
26   
27      do i=1 to &nObs;
28   
29         do j=1 to &nClass;
30            c{j} = int(ranuni(1)*&nLevs);
31         end;
32   
33         d1 = ranuni(&seed);
34         d2 = ranuni(&seed);
35   
36         xBeta = 0.5*((c2<2) - 2*(c1=1) +
36 ! 0.5*c&nClass + 0.05*d1);
37         mu = exp(xBeta);
38   
39         /* Poisson distributions parms */
40         lambda = mu**(2-p)/(phi*(2-p));
41         /* Gamma distribution parms */
42         alpha = (2-p)/(p-1);
43         gamma = phi*(p-1)*(mu**(p-1));
44   
45         rpoi = ranpoi(&seed,lambda);
46         if rpoi=0 then yTweedie=0;
47         else do;
48            yTweedie=0;
49            do j=1 to rpoi;
50            yTweedie = yTweedie + rangam(&seed,alpha);
51            end;
52            yTweedie = yTweedie * gamma;
53         end;
54         output;
55      end;
56   run;
NOTE: The data set WORK.TMP1 has 5000 observations and 8 variables.
NOTE: DATA statement used (Total process time):
      real time           0.01 seconds
      cpu time            0.01 seconds
      
57   ods html5 close;ods listing;

58   

/* EDA */

%let var_char = yTweedie c1 c2 c3 c4 c5 d1 d2;

%put &var_char;


data var_char;
    set tmp1
    (keep= &var_char);     
run;

proc contents data = var_char varnum nodetails noprint 
out=var_char_names (keep=name);
run;

data var_char_names;
    set var_char_names;
    j = _n_;
run;

* Determine the number of observations;
data _NULL_;
    if 0 then set var_char_names nobs=n;
    call symputx('nrows',n);
    stop;
run;

%put &nrows;

%macro do_eda_uni;
%do obs = 1 %to &nrows;

data _null_;
    set var_char_names;
    if j = &obs  then call symputx("var", put(name, 10.));
run;
%if (%upcase(&var)=YTWEEDIE) or (%upcase(&var)=D1) or (%upcase(&var)=D2)   %then %do; 

    ods graphics on;
        proc means data=tmp1 fw=12 printalltypes chartype
            qmethod=os maxdec=2

            mean 
            min 
            max 
            mode 
            range 
            n 
            nmiss   
            p1 
            p5 
            median 
            p95 
            p99 ;
            var &var;
        run;

        title "histograms";
        proc univariate data=tmp1   noprint;
            var &var;
            histogram ;
        run; 
    ods graphics off;
    %end;

    %else %do;
    ods graphics on;
        proc freq data=tmp1
        order=internal;
        tables &var /  scores=table plots(only)=freq;
        run;
    ods graphics off;
    %end;

%end;
%mend do_eda_uni;

%do_eda_uni;

/* Variable Selection Macro: Backwards elimination */
%let p=1.5;
options mlogic;
%macro MdStmt(
        resvar = /*response variable */
       ,expvar = /*list of explanatory variables, separated by ' ' */
       ,clsvar = /*classification variables in the CLASS statement separated by ' ' */
       ,p = 
       );

        ods output Type3=pval(rename=source=parm);
        proc genmod data=tmp1 NAMELEN=50; 
            if _resp_ > 0 then 
            d = 2*(_resp_*(_resp_**(1-&p)-_mean_**(1-&p))/
            (1-&p)-(_resp_**(2-&p)-_mean_**(2-&p))/(2-&p)); 
            else d = 2* _mean_**(2-&p)/(2-&p);
            variance var = _mean_**&p;
            deviance dev = d;
            class &clsvar;  
            model &resvar =  &expvar /link=log type3 scale=pearson;                 
            *scwgt expos;
            title "&resvar = &expvar";  
        run;
        ods output close;
 %mend MdStmt;

144  ods listing close;ods html5 file=stdout options(bitmap_mode='inline') device=png; ods graphics on / outputfmt=png;
NOTE: Writing HTML5 Body file: STDOUT
145  
146  /* Variable Selection Macro: Backwards elimination */
147  %let p=1.5;
148  options mlogic;
149  %macro MdStmt(
150        resvar = /*response variable */
151      ,expvar = /*list of explanatory variables, separated by ' ' */
152        ,clsvar = /*classification variables in the CLASS statement separated by ' ' */
153       ,p =
154      );
155  
156          ods output Type3=pval(rename=source=parm);
157       proc genmod data=tmp1 NAMELEN=50;
158             if _resp_ > 0 then
159             d = 2*(_resp_*(_resp_**(1-&p)-_mean_**(1-&p))/
160            (1-&p)-(_resp_**(2-&p)-_mean_**(2-&p))/(2-&p));
161           else d = 2* _mean_**(2-&p)/(2-&p);
162            variance var = _mean_**&p;
163            deviance dev = d;
164             class &clsvar;  
165              model &resvar =  &expvar /link=log type3 scale=pearson;
166           *scwgt expos;
167              title "&resvar = &expvar";   
168          run;
169          ods output close;
170   %mend MdStmt;
171  ods html5 close;ods listing;

172  

%macro MdSelect(
       var= /*response variable */
       ,intvar= /*initial explanatory variables for full model */
       ,catvar= /*categorical explanatory variables */
       ,slstay= /*criterion for removing variable */
       ,power=
       );
    %let var=%upcase(&var);
    %let intvar=%upcase(&intvar);
    %let catvar=%upcase(&catvar);
    %let power =&power; 
%*-------------------------------------------------------------------------*;
%* Create empty dataset "step" with only one column "parm". It will be *;
%* merged with "pval" from PROC GENMOD by "parm" *;
%*-------------------------------------------------------------------------*;
 proc sql;
    create table step_&var (parm char(9));
 quit;
%*------------------------------------------------------------------------------*;
%* %do %until performs multivariate backward model selection: *;
%* In each iteration: *;
%* 1. Run the logistic regression model *;
%* 2. Update the dataset "step_&var" *;
%* 3. Create &pmax as the maximum p-value, and &varlist as the list of *;
%* variables without the one with the max p-value *;
%* 4. Check whether the max p-value <= &SLSTAY *;
%* 5. If NO, then eliminate the variable with max p-value, repeat step 1 to 4.*;
%* If YES, the loop stops *;
%*------------------------------------------------------------------------------*;
 %let i=1;
 %do %until (&pmax<=&slstay);

    %if &i = 1 %then
        %MdStmt(resvar=&var ,expvar=&intvar, clsvar=&catvar, p=&power); %*initial model;
    %else %do;
        %MdStmt(resvar=&var ,expvar=&varlist, clsvar=&catvar, p=&power); %*reduced model;
    %end;
    proc sort data=step_&var; by parm;
    proc sort data=pval; by parm;
    data step_&var;
        merge step_&var pval;
        by parm;
        p&i=put(ProbChiSq, pvalue6.3);
        drop ProbChiSq ChiSq DF;
    run;
    proc sql noprint;
        select max(ProbChiSq) into :pmax
        from pval;  
        select distinct parm into :varlist separated by ' '
        from pval
        having ProbChiSq^=max(ProbChiSq);
    quit;

    %let i=%eval(&i+1);

 %end;
 proc print data=step_&var;
    title "&var: model selection process";
 run;
%mend MdSelect; 

%MdSelect(var=yTweedie, intvar=c1 c2 c3 c4 c5 d1 d2, catvar=c1 c2 c3 c4 c5, slstay=0.05, power=1.5);