How do I run a random effect zero-inflated poisson model using nlmixed? | SAS FAQ

For this example, we are using the Hospital, Doctor, Patient dataset.
The SAS program containing all the code for this page may be downloaded
here.

Poisson model

Let’s start by just considering a poisson regression model. Unlike the familiar Gaussian distribution which has two parameters \(\mathcal{N}(\mu, \sigma^{2})\), the Poisson distribution is described by a single parameter, \(lambda\) that is both the mean and variance. That is, \(\lambda = E(x)\) and \(\lambda = Var(x) = E(x^{2}) – E(x)^{2}\). This can be a strong assumption. The standard Poisson model is (see section 4.3 in Categorical Data Analysis, 2nd ed. by Alan Agresti):
\begin{equation}
ln[E(y_{i})] = ln(\mu_{i}) = \alpha + \beta x_{i}
\end{equation}
The log link is the canonical link function for the Poisson distribution, and the expected value of the response is modeled. Because the poisson distribution only has one parameter, there is no residual variance (although there are residuals). We can put the model back in the original units as:
\begin{equation}
\lambda_{i} = \mu_{i} = e^{\alpha + \beta x_{i}}
\end{equation}
and we assume the response is distributed as Poisson
\begin{equation}
y_{i} \sim Pois(\lambda_{i})
\end{equation}
The likelihood function for the Poisson model is given by:
\begin{equation}
\mathcal{l} = e^{y_{i}ln(\mu_{i}) – \mu_{i} – ln(y_{i}!)}
\end{equation}
so the log likelihood is simply:
\begin{equation}
\mathcal{ll} = y_{i}ln(\mu_{i}) – \mu_{i} – ln(y_{i}!)
\end{equation}
finally note that the natural logarithm of a n factorial is equivalent the the natural logarithm of the gamma function n plus 1, that is:
\begin{equation}
ln(y_{i}@!) = lgamma(y_{i} + 1)
\end{equation}
Now let’s look at how we can run a Poisson model in SAS. Because proc nlmixed can be twitchy about start values, we will run the model first in proc genmod and then in proc nlmixed. We would not technically need to use proc nlmixed at all at this stage; however, because we do need to use it later, we will use it now to show how each piece is built up to make the final product clearer.

Here is the code and output for the poisson model from proc genmod.

/*read data in from internet*/
filename tmp url 'http://stats.idre.ucla.edu/stat/data/hdp.csv';

data hdp;
infile tmp dlm=',' firstobs=2;
input tumorsize co2 pain wound mobility ntumors nmorphine 
  remission lungcapacity Age Married $ FamilyHx $ SmokingHx $
  Sex $ CancerStage $ LengthofStay WBC RBC BMI IL6 CRP DID $
  Experience School $ Lawsuits HID $ Medicaid;
run;

/*First because we are using proc nlmixed*/
/*we will manually dummy code some variables*/
data hdp;
  set hdp;
  if SmokingHx = "former" then SmokingHx2 = 1;
    else SmokingHx2 = 0;
  if SmokingHx = "never" then SmokingHx3 = 1;
    else SmokingHx3 = 0;
  if CancerStage = "II" then CancerStage2 = 1;
    else CancerStage2 = 0;
  if CancerStage = "III" then CancerStage3 = 1;
    else CancerStage3 = 0;
  if CancerStage = "IV" then CancerStage4 = 1;
    else CancerStage4 = 0;
run;

/*Basic poisson regression model*/

/*using genmod to get start values for nlmixed*/
proc genmod data = hdp;
  model nmorphine = Age LengthofStay BMI 
    CancerStage2 CancerStage3 CancerStage4 / dist=poisson;
run;

Now we do the same thing, just a basic poisson model, but using proc nlmixed.

/*using nlmixed*/
proc nlmixed data = hdp;
  parms b0=.839 b1=-.001 b2=.003 b3=.013 b4=.096 b5=.209 b6=.306;
  mu = exp(b0 + b1*Age + b2*LengthofStay + b3*BMI + 
    b4*CancerStage2 + b5*CancerStage3 + b6*CancerStage4);
  ll = nmorphine*log(mu) - mu - lgamma(nmorphine + 1);
  model nmorphine ~ general(ll);
run;	

Zero-inflated poisson model

In a zero-inflated poisson model, the we observe more zero values than would be expected under a poisson model. The assumption is that these excess zeros are caused by an alternative stochastic process. One way to address this is to explicitly model the probability of being a zero or not, and then modify the estimation of the poisson model based on this. Thus zero-inflated poisson models typically have two portions. The first is a logistic model for being a zero or not, the second is a poisson model for the count. For the probability, we create a linear predictor (call it Xb):
\begin{equation}
Xb_{i} = alpha + beta x_{i}
\end{equation}
of course we can have more than one predictor and they need not be the same predictors as used for the poisson model. Then we can use the inverse logit function to get the (conditional) probability:
\begin{equation}
P(y_{i} = 1 | x_{i}) = p_{i} = \frac{1}{1 + e^{-Xb_{i}}}
\end{equation}
we can use this probability and create a new log likelihood function for the zero-inflated poisson model, that is:
\begin{equation}
y_{i} = \left\{ \begin{array}{rl}
ln(p_{i} + (1 – p_{i})e^{(-\mu_{i})}) &\mbox{if $y_{i} = 0$} \\
y_{i}ln(\mu_{i}) + ln(1 – p_{i}) – \mu_{i} – ln(y_{i}!) &\mbox{if $y_{i} > 0$}
\end{array} \right.
\end{equation}
The log likelihood function changes depending on whether the observed value is zero or not. Also, note that when \(y_{i} = 0\), \(y_{i}ln(\mu_{i})\) and \(ln(y_{i}!)\) are zero and so drop out (although we could have written them if we wanted to). Again, we will run the model twice, first using proc genmod to get start values and then using proc nlmixed to show how the likelihood function is coded.

/*Complication 1: zero-inflated poisson regression model*/	

/*using genmod to get start values for nlmixed*/
proc genmod data = hdp;
  model nmorphine = Age LengthofStay BMI
    CancerStage2 CancerStage3 CancerStage4 / dist=zip;
  zeromodel Age SmokingHx2 SmokingHx3 
    CancerStage2 CancerStage3 CancerStage4 /link = logit;
run;

Now we do the same thing, a zero-inflated poisson model in proc nlmixed.

/*using nlmixed*/
proc nlmixed data = hdp;
  parms b0=1.057 b1=-.002 b2=.009 b3=.012 b4=.024 b5=.053 b6=.132
        a0=-1.677 a1=-.012 a2=.748 a3=1.259 a4=-.808 a5=-2.141 a6=-2.753;
  logit0 = a0 + a1*Age + a2*SmokingHx2 + a3*SmokingHx3 +
    a4*CancerStage2 + a5*CancerStage3 + a6*CancerStage4;
  prob0 = 1 / (1 + exp(-logit0));
  mu = exp(b0 + b1*Age + b2*LengthofStay + b3*BMI + 
    b4*CancerStage2 + b5*CancerStage3 + b6*CancerStage4);
  if nmorphine = 0 then
    ll = log(prob0 + (1 - prob0)*exp(-mu));
  else
    ll = nmorphine*log(mu) + log(1 - prob0) - mu - lgamma(nmorphine + 1);
  model nmorphine ~ general(ll);
run;

Mixed effects poisson models

If the data are clustered, we may want to run a mixed effects poisson model. The model is similar to the original poisson model:
\begin{equation}
ln(E(y_{ij} | u_{j})) = \alpha + \beta x_{ij} + u_{j}
\end{equation}
where
\begin{equation}
u_{j} \sim \mathcal{N}(0, \sigma^{2})
\end{equation}
that is, conditional on the random effect for group j (\(mu_{j}\) we assume the regular poisson model for individual i, and the random effects are assumed to follow a normal distribution with mean 0 (the overall is captured by the intercept, \(alpha\)) and some estimated variance. This can be done in SAS using proc glimmix. We will run the model twice, first using proc glimmix to get start values and then using proc nlmixed to show how the model is setup. For proc glimmix we use the Laplace approximation, which is faster, though (in some cases) slightly less accurate the Gauss Hermite approximation. This is because we are just using proc glimmix for start values and this is not our final model anyway.

/*Complication 2: poisson regression model with random effect*/	

/*Start by sorting data*/
proc sort data = hdp;
  by DID;
run;

/*using glimmix to get start values for nlmixed*/
proc glimmix data = hdp noclprint method=laplace;
  class DID;
  model nmorphine = Age LengthofStay BMI 
    CancerStage2 CancerStage3 CancerStage4 / solution dist=poisson;
  random intercept / subject = DID;
run;

Now we do the same thing, a mixed effects poisson model in proc nlmixed.

/*using nlmixed*/
proc nlmixed data = hdp;
  parms b0=.7465 b1=-.0018 b2=.0056 b3=.013 b4=.0792 b5=.2049 b6=.3166
    sigma2=.27;
  mu = exp(b0 + u + b1*Age + b2*LengthofStay + b3*BMI + 
    b4*CancerStage2 + b5*CancerStage3 + b6*CancerStage4);
  ll = nmorphine*log(mu) - mu - lgamma(nmorphine + 1);
  model nmorphine ~ general(ll);
  random u ~ normal(0, sigma2) subject=DID;
run;

Putting it together: mixed effects zero-inflated poisson model

There is no way to run a mixed effects model in proc genmod and there is no way to include zero-inflation in proc glimmix so in SAS, this model (to our knowlege) can currently only be fit using proc nlmixed. Because this model cannot be fit with another SAS procedure, there is no direct way to get starting values for all the parameters. What we did was for the fixed effects, use the estimates from the zero-inflated poisson model run in proc genmod and for the random effects, use the estimates from proc glimmix. Combined, we at least have some (hopefully reasonable) start values for all the parameters. The start values do not need to be perfect, but usually the better they are, the faster the model will converge. Also, in some cases inadequate start values may result in the optimizer searching in a region of the parameter space where it cannot progress. In this case, there are two main options. First you can try tuning your start values. Alternately, you can try a different algorithm that may be able to progress out of the region the other algorithm got “stuck in”. Essentially, change where you are looking (different start values) or change how you are looking (different algorithm). Note that for variance components (in this case just sigma2), it is often better to overestimate than to underestimate them.

Unfortunately, using the coefficients from proc genmod and proc glmmix as start values will not guarantee that the model in proc nlmixed will converge. We see that using our start values as they are will cause SAS to return an error that the model has not converged:

/*Putting it together: zero-inflated poisson regression model with random effect*/

/*to our knowledge this is only possible in SAS using nlmixed*/
proc nlmixed data = hdp method=gauss qpoints=25;
  parms b0=1.0569 b1=-0.0020 b2=0.0094 b3=0.0122 b4=0.0239 b5=0.0530 b6=0.1322
        a0=-1.6766 a1=-0.0116 a2=0.7479 a3=1.2588 a4=-0.8077 a5=-2.1410 a6=-2.7530 	
        sigma2=0.2414;
  logit0 = a0 + a1*Age + a2*SmokingHx2 + a3*SmokingHx3 +
    a4*CancerStage2 + a5*CancerStage3 + a6*CancerStage4;
  prob0 = 1 / (1 + exp(-logit0));
  mu = exp(b0 + u + b1*Age + b2*LengthofStay + b3*BMI + 
    b4*CancerStage2 + b5*CancerStage3 + b6*CancerStage4);
  if nmorphine = 0 then
    ll = log(prob0 + (1 - prob0)*exp(-mu));
  else
	ll = nmorphine*log(mu) + log(1 - prob0) - mu - lgamma(nmorphine + 1);
  model nmorphine ~ general(ll);
  random u ~ normal(0, sigma2) subject=DID;
run;

It is difficult to tell what has gone wrong with this error. At times, using the itdetails option on the proc nlmixed statement can help diagnose the error, as you can see how the parameters change at each iteration and look for problematic parameters. The detailed iteration history provides no help in this case. We recommend varying the starting value for the random effect first, as we saw that the proc nlmixed model converges fine if the random effect is excluded (see “Zero-inflated poisson model” above). Indeed, changing the start value for sigma2 from 0.2414 to 0.3 will allow the model to converge:

proc nlmixed data = hdp method=gauss qpoints=25;
  parms b0=1.0569 b1=-0.0020 b2=0.0094 	 b3=0.0122 b4=0.0239 b5=0.0530 b6=0.1322
        a0=-1.6766 a1=-0.0116 a2=0.7479 a3=1.2588 a4=-0.8077 a5=-2.1410 a6=-2.7530 	
        sigma2=0.3;
  logit0 = a0 + a1*Age + a2*SmokingHx2 + a3*SmokingHx3 +
    a4*CancerStage2 + a5*CancerStage3 + a6*CancerStage4;
  prob0 = 1 / (1 + exp(-logit0));
  mu = exp(b0 + u + b1*Age + b2*LengthofStay + b3*BMI + 
    b4*CancerStage2 + b5*CancerStage3 + b6*CancerStage4);
  if nmorphine = 0 then
    ll = log(prob0 + (1 - prob0)*exp(-mu));
  else
	ll = nmorphine*log(mu) + log(1 - prob0) - mu - lgamma(nmorphine + 1);
  model nmorphine ~ general(ll);
  random u ~ normal(0, sigma2) subject=DID;
run;

Extending the random effects

Previously we fit a mixed effects zero-inflated poisson model. Now we are going to fit the same model, but include a random slope and covariance between the intercept and slope. Again our procedure for selecting start values was to use the prior proc nlmixed model estimates, and run a model in proc glimmix with a random intercept and random slope, and use the variance estimates from proc glimmix for the variance components. Once the model converged, we updated the start values in the code to be close to the final estimates to speed future runs of the code.

/*Extending the random effects 1*/
proc nlmixed data = hdp method=gauss qpoints=25;
  parms b0=0.8364 b1=-0.00242 b2=0.01052 b3=0.01286 b4=0.04828 b5=0.1398 b6=0.2496
        a0=-4.4730 a1=-0.03124 a2=3.5560 a3=4.5860 a4=-1.3089 a5=-7.3938  a6=-4.9867
        sigma2_u0=0.2414 sigma_cov = -0.00044 sigma2_u1=0.000016;
  bounds sigma2_u0 sigma2_u1 >= 0;
  logit0 = a0 +  a1*Age + a2*SmokingHx2 + a3*SmokingHx3 +
    a4*CancerStage2 + a5*CancerStage3 + a6*CancerStage4;
  prob0 = 1 / (1 + exp(-logit0));
  mu = exp(b0 + u0 + b1*Age + b2*LengthofStay + b3*BMI + u1*BMI + 
    b4*CancerStage2 + b5*CancerStage3 + b6*CancerStage4);
  if nmorphine = 0 then
    ll = log(prob0 + (1 - prob0)*exp(-mu));
  else
	ll = nmorphine*log(mu) + log(1 - prob0) - mu - lgamma(nmorphine + 1);
  model nmorphine ~ general(ll);
  random u0 u1 ~ normal([0, 0], [sigma2_u0, sigma_cov, sigma2_u1]) subject=DID;
run;

We can also fit a model that allows a random intercept for the zero-inflation portion. We assume that the two random intercepts (for the zero-inflation and poisson portion) are uncorrelated. The start value for sigma2_u1 of 1 was chosen arbitrarily, but results in model convergence.

/*Extending the random effects 2*/
proc nlmixed data = hdp method=gauss qpoints=25 maxfunc = 3000 maxiter = 1000;
  parms b0=0.8463 b1=-0.00245 b2=0.01118 b3=0.01246 b4=0.04719 b5=0.1427 b6=0.2496
        a0=-4.5314 a1=-0.03164 a2=3.6122 a3=4.6438 a4=-1.3091 a5=-3.4101 a6=-4.8750
        sigma2_u0=0.2300 sigma2_u1=1;
  bounds sigma2_u0 sigma2_u1 >= 0;
  logit0 = a0 + u1 + a1*Age + a2*SmokingHx2 + a3*SmokingHx3 +
    a4*CancerStage2 + a5*CancerStage3 + a6*CancerStage4;
  prob0 = 1 / (1 + exp(-logit0));
  mu = exp(b0 + u0 + b1*Age + b2*LengthofStay + b3*BMI + 
    b4*CancerStage2 + b5*CancerStage3 + b6*CancerStage4);
  if nmorphine = 0 then
    ll = log(prob0 + (1 - prob0)*exp(-mu));
  else
	ll = nmorphine*log(mu) + log(1 - prob0) - mu - lgamma(nmorphine + 1);
  model nmorphine ~ general(ll);
  random u0 u1 ~ normal([0, 0], [sigma2_u0, 0, sigma2_u1]) subject=DID;
run;

Comments

When other SAS procedures will fit a model, they are often easier to use and faster. However, for uncommon or new models that are not implemented in SAS, using proc nlmixed allows a great deal of flexibility in model specification.

References

1. Voronca, D.C. et al.(2014). Analysis of Zero Inflated Longitudinal Data using Proc Nlmixed. Conference Proceedings: SouthEast SAS Users Group.
2. Min, Y. and Agresti, A.(2005). Random effect models for repeated measures of zero-inflated count data. Statistical Modeling; 5:1-19.
3. Runze, A.B. and Zucker, R.A.(2012). Statistical models for longitudinal zero-inflated count data with application to the substance abuse fields. Stat Med; 31(29):4074-4086.