Negative binomial models are count regression models that work with overdispersed data, i.e., count data in which the variance is greater than the mean. This FAQ page will show how to use proc nlmixed to analyze negative binomial models with random effects. We will look at three models beginning with an ordinary negative binomial without random effects, a negative binomial model with random intercepts and finally a model with both random intercepts and random slopes.
All of the models use the same dataset, epilepsy, and use the same predictor variables, treat, coded 1 for progabide and 0 for placebo, lbas, log of the baseline seizure rate, and visit, doctor visits linearly coded as (-.3 -.1 .1 .3). The response variable is seizure which is the number of seizures. For each of the models we will give only partial output, primarily estimates, standard errors, Wald tests, p-values, log-likelihoods, AICs and BICs to allow comparison among the three models.
Single-level negative binomial regression
The setup for negative binomial regression in proc nlmixed is relatively straightforward with xb being the linear predictor with parameters, b0 (the intercept), b1, b2 and b3 as the regression coefficients. There is an additional parameter, alpha, which captures the overdispersion. The model is parameterized using m which is 1/alpha By using proc nlmixed we will obtain an iterated maximum likelihood solution for the model.
The nlmixed procedure is a very flexible command, and there many different ways that you could estimate any given model. We have chosen to use the general.
options nocenter;
/* single-level negative binomial regression */
proc nlmixed data='D:\data\epilepsy.sas7bdat';
xb=b0 + b1*treat + b2*lbas + b3*visit;
mu = exp(xb);
m = 1/alpha;
ll = lgamma(seizures+m)-lgamma(seizures+1)-lgamma(m)
+ seizures*log(alpha*mu)-(seizures+m)*log(1+alpha*mu);
model seizures~ general(ll);
run;
/* partial output */
Fit Statistics
-2 Log Likelihood 1304.4
AIC (smaller is better) 1314.4
AICC (smaller is better) 1314.7
BIC (smaller is better) 1331.7
Parameter Estimates
Standard
Parameter Estimate Error DF t Value Pr > |t| Alpha Lower Upper Gradient
b0 1.8947 0.07229 236 26.21 <.0001 0.05 1.7522 2.0371 4.175E-6
b1 -0.2595 0.1002 236 -2.59 0.0102 0.05 -0.4570 -0.062 0.000069
b2 1.0439 0.06426 236 16.25 <.0001 0.05 0.9173 1.1705 0.000267
b3 -0.2902 0.2250 236 -1.29 0.1984 0.05 -0.7334 0.1531 -2.13E-6
alpha 0.3915 0.05353 236 7.31 <.0001 0.05 0.2861 0.4970 0.000105
Negative binomial regression with random intercepts
Our second model is a negative binomial model predicting seizures from treat, lbas and visits.
/* random intercept negative binomial */
proc nlmixed data='D:\data\epilepsy.sas7bdat';
xb=b0 + b1*treat + b2*lbas + b3*visit + u;
mu = exp(xb);
m = 1/alpha;
ll = lgamma(seizures+m)-lgamma(seizures+1)-lgamma(m)
+ seizures*log(alpha*mu)-(seizures+m)*log(1+alpha*mu);
model seizures~ general(ll);
random u ~ normal(0,s2u) subject=subject;
run;
/* partial output */
Fit Statistics
-2 Log Likelihood 1252.6
AIC (smaller is better) 1264.6
AICC (smaller is better) 1265.0
BIC (smaller is better) 1277.1
Parameter Estimates
Standard
Parameter Estimate Error DF t Value Pr > |t| Alpha Lower Upper Gradient
b0 1.8189 0.1084 58 16.78 <.0001 0.05 1.6019 2.0358 -0.00023
b1 -0.3351 0.1517 58 -2.21 0.0312 0.05 -0.6387 -0.0314 -0.00069
b2 1.0114 0.1011 58 10.01 <.0001 0.05 0.8090 1.2137 -0.00014
b3 -0.2709 0.1666 58 -1.63 0.1095 0.05 -0.6044 0.0627 -0.00055
alpha 0.1343 0.03161 58 4.25 <.0001 0.05 0.07100 0.1975 -0.00152
s2u 0.2414 0.06334 58 3.81 0.0003 0.05 0.1147 0.3682 -0.00137
The second approach is to compute a likelihood from the probability, take the log of the likelihood and then compute a general log-likelihood function using general.
/* random intercept and slope nbreg */
proc nlmixed data='D:\data\epilepsy.sas7bdat';
parms b0 1.5 b1 -.3 b2 1 b3 -.2 alpha .1 s2u .1 s2v .5 cov 0;
xb=b0 + b1*treat + b2*lbas + b3*(visit+rv) + u;
mu = exp(xb);
m = 1/alpha;
ll = lgamma(seizures+m)-lgamma(seizures+1)-lgamma(m)
+ seizures*log(alpha*mu)-(seizures+m)*log(1+alpha*mu);
model seizures~ general(ll);
random u rv ~ normal([0,0],[s2u,cov,s2v]) subject=subject;
run;
[partial output]
Fit Statistics
-2 Log Likelihood 1252.9
AIC (smaller is better) 1268.9
AICC (smaller is better) 1269.6
BIC (smaller is better) 1285.6
Parameter Estimates
Standard
Parameter Estimate Error DF t Value Pr > |t| Alpha Lower Upper Gradient
b0 1.8188 0.1081 57 16.82 <.0001 0.05 1.6023 2.0354 0.004175
b1 -0.3348 0.1513 57 -2.21 0.0310 0.05 -0.6378 -0.03176 0.00519
b2 1.0114 0.1008 57 10.03 <.0001 0.05 0.8095 1.2134 -0.0004
b3 -0.2709 0.1665 57 -1.63 0.1093 0.05 -0.6044 0.06257 -0.00142
alpha 0.1340 0.03158 57 4.24 <.0001 0.05 0.07072 0.1972 -0.01416
s2u 0.2069 0.05850 57 3.54 0.0008 0.05 0.08973 0.3240 0.000486
s2v 0.4930 0.004293 57 114.84 <.0001 0.05 0.4844 0.5016 0.000036
cov 0.005711 0.03169 57 0.18 0.8577 0.05 -0.05776 0.06918 -0.00026