User:Statcompute/BlogEntry: 2007 April 25 23:56:22 EDT

From sasCommunity

/*************************************************
* MODEL ZERO-INFLATED POISSON REGRESSION USING   *
* PROC NLMIXED AND PROC COUNTREG                 *
*************************************************/

data data (keep = x1 x2 y);
  do i = 1 to 1000;
    x1 = rannor(1);
    x2 = rannor(0);

    eta_lambda = 0.1 + x1 * 0.2 + x2 * 0.3;
    lambda = exp(eta_lambda);

    eta_prob = 0.4 + 0.5 * x1 + 0.6 * x2;
    p_0 = exp(eta_prob) / (1 + exp(eta_prob));
    zero_inflate = ranbin(0, 1, p_0);

    if zero_inflate = 1 then y = 0;
    else y = ranpoi(0, lambda);
    output;
  end;
run;

proc nlmixed data = data corr cov;
  parms poisn_beta0 = 0 poisn_beta1 = 0 poisn_beta2 = 0
        logit_beta0 = 0 logit_beta1 = 0 logit_beta2 = 0;

  eta_prob = logit_beta0 + logit_beta1 * x1 + logit_beta2 * x2;
  p_0 = exp(eta_prob) / (1 + exp(eta_prob));

  eta_lambda = poisn_beta0 + poisn_beta1 * x1 + poisn_beta2 * x2;
  lambda = exp(eta_lambda);

  if y = 0 then loglike = log(p_0 + (1 - p_0) * exp(-lambda));
  else loglike = log(1 - p_0) + y * log(lambda) - lambda - lgamma(y + 1);

  model y ~ general(loglike);
run;
/*                         Parameter Estimates
                           Standard
  Parameter     Estimate      Error     DF   t Value   Pr > |t|
  poisn_beta0     0.1965    0.06748   1000      2.91     0.0037
  poisn_beta1     0.1183    0.06279   1000      1.88     0.0598
  poisn_beta2     0.3587    0.07127   1000      5.03     <.0001
  logit_beta0     0.3539     0.1253   1000      2.82     0.0048
  logit_beta1     0.5757     0.1222   1000      4.71     <.0001
  logit_beta2     0.7726     0.1538   1000      5.02     <.0001
  ... ...
*/

proc countreg data = data type = zip printall;
  model y = x1 x2 / zi(link = logistic, var = x1 x2);
run;
/*                         Parameter Estimates
                                       Standard                 Approx
  Parameter            Estimate           Error    t Value    Pr > |t|
  Intercept            0.196469        0.067476       2.91      0.0036
  x1                   0.118323        0.062792       1.88      0.0595
  x2                   0.358667        0.071274       5.03      <.0001
  Inf_Intercept        0.353878        0.125317       2.82      0.0047
  Inf_x1               0.575729        0.122236       4.71      <.0001
  Inf_x2               0.772563        0.153843       5.02      <.0001
  ... ...
*/