Fitting a Scaled Tweedie Model with Regressors

/*--------------------------------------------------------------

                    SAS Sample Library

        Name: sevex05.sas
 Description: Example Program from SAS/ETS User's Guide,
              The SEVERITY Procedure
       Title: Fitting a Scaled Tweedie Model with Regressors
     Product: SAS/ETS Software
        Keys: Severity Distribution Modeling
        PROC: SEVERITY
       Notes:

--------------------------------------------------------------*/


/*--- Simulate a Tweedie sample that is affected by regressors ---*/
%let samplesize=300;
data test_sevtw(keep=y x1-x3
                label='A Tweedie Sample Affected by Regressors');
   array x{*} x1-x3;
   array b{4} _TEMPORARY_ (5 0.25 -1 3);
   call streaminit(45678);
   label y='Response Influenced by Regressors';

   Phi = 0.5;
   P = 1.75;
   do n = 1 to &samplesize;
      Mu = 0;
      do i = 1 to dim(x);
         x(i) = rand('UNIFORM');
         Mu = Mu + b(i+1) * x(i);
      end;
      Mu = b(1) * exp(Mu); /* b(1) is base value of mean */

      cdf = rand('UNIFORM');
      y = quantile('TWEEDIE', cdf, P, Mu, Phi);
      output;
   end;
run;


/*--- Fit the scale parameter version of the Tweedie distribution ---*/
proc severity data=test_sevtw outest=estw covout print=all plots=none;
   loss y;
   scalemodel x1-x3;

   dist stweedie;
   nloptions tech=quanew;
run;


/*--- Compute the estimate, standard error, and p-value of the mean ---*/
data Mu0(keep=Parameter Estimate Stderr T probT);
   Parameter='Mu0';
   label Estimate='Estimate' Stderr='Standard;Error'
         T='t Value' probT='Approx;Pr > |t|';

   array cov{3,3} _temporary_;
   array parms{3} _temporary_;
   array parmrow{3} theta lambda p;
   set estw(where=(_MODEL_='stweedie' and
                   (_type_='COV' or _type_='EST'))) end=eof;
   n = _n_-1;
   if (n = 0) then do;
      do i=1 to 3;
         parms[i] = parmrow[i];
      end;
   end;
   else do;
      do i=1 to n;
         cov[n,i] = parmrow[i];
      end;
   end;
   if (eof or n = 3) then do;
      mu0 = parms[1]*parms[2]*(2-parms[3])/(parms[3]-1);
      Estimate = mu0;

      a = mu0/parms[1]; /* dMu0_T0 */
      b = mu0/parms[2]; /* dMu0_L */
      c = -mu0/((parms[3]-1)*(2-parms[3])); /* dMu0_P */

      varMu0 = a * (a*cov[1,1]+b*cov[2,1]+c*cov[3,1]) +
               b * (a*cov[2,1]+b*cov[2,2]+c*cov[3,2]) +
               c * (a*cov[3,1]+b*cov[3,2]+c*cov[3,3]);

      Stderr = sqrt(varMu0);

      df = &samplesize-6;
      T = mu0/Stderr;
      probT = (1-probt(T, df))*2;
      output;
      stop;
   end;
run;

proc print data=Mu0 noobs s=';';
run;