CQLIM Procedure

This example uses the CQLIM procedure to process a small data table.

The following DATA step generates a data set that contains 5,000 observations from a logistic regression model. The model contains seven variables.

data nuts_ex;
   call streaminit(784923);
   array vars x1-x7;
   array slopes{8}  (0.03 0.04 0.02 0.04 -0.03 -0.05 -0.03 0.01);
   intercept = 2;

   do i = 1 to 5000;
      sum_xb = intercept;
      do j=1 to 7;
         vars[j] = rand('NORMAL',0,1);
         sum_xb = sum_xb + slopes[j] * vars[j];
         end;
      sum_xb = sum_xb + vars[1] * vars[1] * slopes[8];
      u = rand('UNIFORM');
      z = sum_xb + log(u) - log(1-u);
      if z > 0 then y = 1;
      if z <= 0 then y = 0;
      output;
      end;
   keep y x1-x7;
run;

You can load the data set nuts_ex into a data table in your session that is associated with the mylib libref. The DATA step assumes that your libref is named mylib, but you can substitute any appropriately defined libref.

data mylib.nuts_ex;
   set nuts_ex;
run;

The following statements estimate a logistic regression model in the Bayesian framework by using the No-U-Turn Sampler (NUTS). The model includes a second-order term for the covariate x1. The NCHAIN= option controls the number of Markov chains, and the NSAMPLE= option controls the size of the posterior sample per chain. NUTS is specified in the SAMPLER= option, and several values controlling the sampler can be set as arguments to SAMPLER=NUTS. The OUTPOST= option creates a CAS table that contains the posterior sample, which in this case is named nuts_post.

proc cqlim data = mylib.nuts_ex;
   model y = x1-x7 x1*x1 / discrete(distribution = logistic);
   bayes seed = 72342 nsample = 2000 nchain = 4 sampler = nuts(ntune = 2000)
         outpost = mylib.nuts_post;
   prior intercept ~ normal(mean = 0, sd = 1000);
   prior x1-x7 ~ normal(mean = 0, sd = 100);
   prior x1___x1 ~ normal(mean = 0, sd = 10);
run;

Output 12.4.1 shows the Bayesian estimation results for the logistic regression model. The "Posterior Summaries" table contains two point estimates of the parameters from the posterior samples: the posterior mean and the posterior median (50th percentile). The "MCMC Diagnostics Summaries" table contains basic convergence diagnostics to check whether the Markov chain has converged and whether the sample size is sufficient.

The following sequence of statements computes the posterior sample of the sample mean marginal effects for each of the covariates in the model. In binary choice models, the marginal effect of the jth covariate on the ith observation is given by the derivative of i’s choice probability with respect to the jth covariate, evaluated at i’s value of that covariate. For logistic regression, this is given by

The model includes a squared term for the covariate x1. Rewrite as to make the squared term explicit. Then the marginal effect of the covariate x1 is given by

The sample mean marginal effect is the mean of these individual marginal effects over all observations in the data set. These mean marginal effects must be computed for each iteration in the posterior sample to obtain the posterior distribution of the mean marginal effects.

The first set of statements computes the sample mean marginal effects for each iteration in the posterior distribution:

proc cas;
   fedSql.execDirect / query =
      "create table nuts_marginaleffects {options replace=true} as
       select avg(gradprob * (b_x1 + 2 * x1 * b_x1___x1)) as me_x1,
              avg(gradprob * b_x2) as me_x2,
              avg(gradprob * b_x3) as me_x3,
              avg(gradprob * b_x4) as me_x4,
              avg(gradprob * b_x5) as me_x5,
              avg(gradprob * b_x6) as me_x6,
              avg(gradprob * b_x7) as me_x7 from
       (
          select exp(xbeta) / (1 + exp(xbeta))**2 as gradprob, * from
          (
             select b.Intercept + a.x1 * b.x1 + a.x2 * b.x2 + a.x3 * b.x3 +
                    a.x4 * b.x4 + a.x5 * b.x5 + a.x6 * b.x6 + a.x7 * b.x7 +
                    a.x1 * a.x1 * b.x1___x1 as xbeta,
                    a.*, b.Chain, b.Iteration, b.Intercept, b.x1___x1 as b_x1___x1,
                    b.x1 as b_x1, b.x2 as b_x2, b.x3 as b_x3, b.x4 as b_x4,
                    b.x5 as b_x5, b.x6 as b_x6, b.x7 as b_x7
             from nuts_ex as a, nuts_post as b
          ) as c
       ) as d
       group by Chain, Iteration;";
run; quit;

The last set of statements uses several actions to compute summary statistics of the posterior distribution of the mean marginal effects, including posterior means, standard deviations, and percentiles:

proc cas;
   simple.summary /  table = {name = "nuts_marginaleffects"},
      inputs = {"me_x1", "me_x2", "me_x3", "me_x4", "me_x5", "me_x6", "me_x7"},
      subSet = {"N", "mean", "std"}
      casout = {name = "nuts_me_summary", replace = True};
   percentile.percentile /  table = {name = "nuts_marginaleffects"},
      inputs = {"me_x1", "me_x2", "me_x3", "me_x4", "me_x5", "me_x6", "me_x7"},
      values = {2.5, 25, 50, 75, 97.5}
      casout = {name = "nuts_me_percentiles", replace = True};
   transpose.transpose /  table = {name = "nuts_me_percentiles",
      groupby = {"_Column_"}}, transpose = {"_Value_"}, prefix = "Pct_",
      id = {"_Pctl_"},  casout = {name = "nuts_me_percentiles", replace = True};
   table.alterTable /  name = "nuts_me_percentiles",
      drop = "_NAME_",
      columns = {{name = "_Column_", rename = "Parameter"},
                 {name = "Pct_2.5", rename = "Pct_2_5",
                  label = "2.5%", format = "D8."},
                 {name = "Pct_25", label = "25%", format = "D8."},
                 {name = "Pct_50", label = "50%", format = "D8."},
                 {name = "Pct_75", label = "75%", format = "D8."},
                 {name = "Pct_97.5", rename = "Pct_97_5",
                  label = "97.5%", format = "D8."}};
   table.alterTable /  name = "nuts_me_summary",
      columns = {{name = "_Column_", rename = "Parameter"},
                 {name = "_NObs_", rename = "N"},
                 {name = "_Mean_", rename = "Mean", format = "D8."},
                 {name = "_Std_", rename = "Std",
                  label = "Standard Deviation", format = "D8."}};
   fedSql.execDirect / query =
      "create table nuts_me_post_summary{options replace = true} as
       select a.*, b.Pct_2_5, b.Pct_25, b.Pct_50, b.Pct_75, b.Pct_97_5 from
       nuts_me_summary as a, nuts_me_percentiles as b
       where a.Parameter = b.Parameter";
    table.fetch / table = {name = "nuts_me_post_summary"},
       index = False, sortby = "Parameter";
run;
quit;

Output 12.4.2 shows the posterior summary of the mean marginal effects.

An advantage of the Bayesian framework is its ability to compute the posterior probability that a particular hypothesis is true. The following statements compute the posterior probability that each mean marginal effect is positive.

data mylib.nuts_me_probs;
   set mylib.nuts_marginaleffects;
   if me_x1 > 0 then x1 = 1; else x1 = 0;
   if me_x2 > 0 then x2 = 1; else x2 = 0;
   if me_x3 > 0 then x3 = 1; else x3 = 0;
   if me_x4 > 0 then x4 = 1; else x4 = 0;
   if me_x5 > 0 then x5 = 1; else x5 = 0;
   if me_x6 > 0 then x6 = 1; else x6 = 0;
   if me_x7 > 0 then x7 = 1; else x7 = 0;
   keep x1 x2 x3 x4 x5 x6 x7;
run;
quit;

proc cas;
   simple.summary /
      table = {name = "nuts_me_probs"},
      inputs = {"x1", "x2", "x3", "x4", "x5", "x6", "x7"},
      subSet = {"mean"},
      casout = {name = "nuts_me_post_probs", replace = True};
   table.alterTable /  name = "nuts_me_post_probs",
      columns = {{name = "_Mean_", rename = "Prob",
                  label = "P(ME > 0)", format = "7.4"},
                 {name = "_Column_", rename = "Parameter"}};
   table.fetch / table = {name = "nuts_me_post_probs"},
      index = False, sortby = "Parameter";
run;
quit;

Output 12.4.3 shows the posterior probability that each mean marginal effect is greater than zero.