This FAQ page will show how a number of simple linear and nonlinear models can be coded using SAS proc nlmixed. What is meant by “simple” here is that all of the models are fixed effects only with no random effects. All of the models shown can be estimated using specific commands in SAS, for example the binary logistic model can be estimated using proc logistic or proc genmod. In fact it is much easier to run these commands using the specific procedures. The purpose of this page is allow you to be acquainted with these simpler nlmixed, to see how they work and how they’re parameterized. With that knowledge you will have a good foundation for building more complex nonlinear mixed models.
All of the models use the same dataset, hsbdemo.sas7bdat, and the same two predictor variables, read and female. Models will use different response variables depending upon the type of response variable that is appropriate. For each of the models we will give only partial output, primarily estimates, standard errors, wald tests and p-values so that you will be able to compare the results with specific commands in SAS.
A short note about starting values. You give starting values using the parms statement. Some models are very sensitive to starting values and will not converge unless given good values. Other models are very tolerant and will work properly with starting values for all parameters set to zero. There are even models that don’t require you to set the starting values, i.e., no parms statement. They get starting values automatically from the nlmixed procedure.
Ordinary least squares regression
We begin with an ordinary least squares regression predicting write from read and female. The setup for nlmixed is very straightforward with xb being the linear predictor with parameters b0 (the intercept), b1 and b2 as the regression coefficients. There is one additional parameter, s2e, which captures the residual variability and would be equivalent to the mean square error in a traditional OLS regression. By using nlmixed we will obtain an iterated maximum likelihood solution for the model.
proc nlmixed data='D:datahsbdemo.sas7bdat';
xb=b0+b1*read+b2*female;
model write ~ normal(xb,s2e);
run;
/* partial output */
Parameter Estimates
Standard
Parameter Estimate Error DF t Value Pr > |t| Alpha Lower Upper Gradient
b0 20.2284 2.6933 200 7.51 <.0001 0.05 14.9174 25.5393 -1.76E-6
b1 0.5659 0.04901 200 11.55 <.0001 0.05 0.4692 0.6625 -0.00011
b2 5.4869 1.0066 200 5.45 <.0001 0.05 3.5019 7.4719 -1.04E-6
s2e 50.1128 5.0113 200 10.00 <.0001 0.05 40.2310 59.9945 1.12E-8
Thus, our model looks like this, yhat = 20.2284 + .5659*read + 5.4869*female.
Binary logistic regression
Our second model is a binary logistic regression predicting honors from read and female. The variable honors means that the student has been selected to participate in the honors English program. We will estimate this model two different ways; first, computing the expected probability and using the binary distribution as shown below.
proc nlmixed data='D:datahsbdemo.sas7bdat';
parms b0=5 b1=0 b2=0;
xb=b0+b1*read+b2*female;
prob = exp(xb)/(1+exp(xb));
model honors ~ binary(prob);
run;
/* partial output */
Parameter Estimates
Standard
Parameter Estimate Error DF t Value Pr > |t| Alpha Lower Upper Gradient
b0 -9.6034 1.4264 200 -6.73 <.0001 0.05 -12.4162 -6.7906 -2.71E-9
b1 0.1444 0.02333 200 6.19 <.0001 0.05 0.09835 0.1904 -1.16E-6
b2 1.1209 0.4081 200 2.75 0.0066 0.05 0.3162 1.9257 4.077E-8
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..
proc nlmixed data='D:datahsbdemo.sas7bdat';
parms b0=0 b1=0 b2=0;
xb=b0+b1*read+b2*female;
prob = exp(xb)/(1+exp(xb));
liklhd = (prob**honors)*((1-prob)**(1-honors));
ll = log(liklhd);
model honors ~ general(ll);
run;
[partial output]
Parameter Estimates
Standard
Parameter Estimate Error DF t Value Pr > |t| Alpha Lower Upper Gradient
b0 -9.6034 1.4264 200 -6.73 <.0001 0.05 -12.4162 -6.7906 2.854E-7
b1 0.1444 0.02333 200 6.19 <.0001 0.05 0.09835 0.1904 0.000016
b2 1.1209 0.4081 200 2.75 0.0066 0.05 0.3162 1.9257 -8.69E-8
Probit regression
For the probit regression model we will use the same response variable, honors, as in the binary logit model above.
proc nlmixed data='D:datahsbdemo.sas7bdat';
parms b0=0 b1=0 b2=0;
xb=b0+b1*read+b2*female;
prob = probnorm(xb);
if honors=0 then liklhd=1-prob;
else liklhd=prob;
ll = log(liklhd);
model honors ~ general(ll);
run;
[partial output]
Parameter Estimates
Standard
Parameter Estimate Error DF t Value Pr > |t| Alpha Lower Upper Gradient
b0 -5.6720 0.7798 200 -7.27 <.0001 0.05 -7.2098 -4.1343 2.66E-6
b1 0.08560 0.01301 200 6.58 <.0001 0.05 0.05996 0.1113 0.000128
b2 0.6340 0.2301 200 2.76 0.0064 0.05 0.1803 1.0877 1.713E-6
Ordered logistic regression
We will use ses as the response variable for the ordered logistic regression. The variable ses is ordered with values 1, 2 and 3. The model given below estimates a proportional odds ordered logistic model. Some statistical software calls the parameters a1 and a2 cut points or thresholds, while other packages will parameterize the model with intercepts.
proc nlmixed data='D:datahsbdemo.sas7bdat';
parms b1=0 b2=0 a1=1 a2=1;
xb=b1*read+b2*female;
if ses=1 then liklhd=1/(1+exp(-a1+xb));
if ses=2 then liklhd=1/(1+exp(-a2+xb))-1/(1+exp(-a1+xb));
if ses=3 then liklhd=1-1/(1+exp(-a2+xb));
ll = log(liklhd);
model ses ~ general(ll);
run;
[partial output]
Parameter Estimates
Standard
Parameter Estimate Error DF t Value Pr > |t| Alpha Lower Upper Gradient
b1 0.05714 0.01386 200 4.12 <.0001 0.05 0.02980 0.08448 -0.00637
b2 -0.4195 0.2715 200 -1.55 0.1239 0.05 -0.9549 0.1159 -0.00004
a1 1.4774 0.7343 200 2.01 0.0456 0.05 0.02942 2.9254 0.000131
a2 3.7306 0.7791 200 4.79 <.0001 0.05 2.1942 5.2669 -4.58E-
Generalized ordered logistic regression
Generalized ordered logit estimates separate intercepts and coefficients for each equation and so therefore does not have a proportional odds assumption. In this example, we use the same variables as in the ordered logistic regression above, however, we have to estimate several more parameters then for that model.
proc nlmixed data='D:datahsbdemo.sas7bdat';
parms b01=-1 b02=-2 b11=0 b21=0 b12=0 b22=0;
xb1=b01+b11*read+b12*female;
xb2=b02+b21*read+b22*female;
if ses=1 then liklhd=1/(1+exp(xb1));
if ses=2 then liklhd=1/(1+exp(xb2))-1/(1+exp(xb1));
if ses=3 then liklhd=1-1/(1+exp(xb2));
ll = log(liklhd);
model ses ~ general(ll);
estimate 'b11-b12' b01-b12;
estimate 'b21-b22' b21-b22;
run;
[partial output]
Parameter Estimates
Standard
Parameter Estimate Error DF t Value Pr > |t| Alpha Lower Upper Gradient
b01 -1.1132 0.9807 200 -1.14 0.2577 0.05 -3.0470 0.8206 0.000022
b02 -3.9737 0.9294 200 -4.28 <.0001 0.05 -5.8064 -2.1410 0.000032
b11 0.05335 0.01887 200 2.83 0.0052 0.05 0.01614 0.09057 0.000953
b21 0.05963 0.01636 200 3.64 0.0003 0.05 0.02736 0.09189 0.001598
b12 -0.7005 0.3591 200 -1.95 0.0525 0.05 -1.4085 0.007610 0.000026
b22 -0.2037 0.3232 200 -0.63 0.5292 0.05 -0.8410 0.4336 0.00001
Additional Estimates
Standard
Label Estimate Error DF t Value Pr > |t| Alpha Lower Upper
b11-b12 -0.4127 1.1275 200 -0.37 0.7147 0.05 -2.6361 1.8106
b21-b22 0.2633 0.3229 200 0.82 0.4158 0.05 -0.3735 0.900
Traditional statistical procedures will often organize the output by equation. The example output below shows one way the output may be displayed.
Standard
Parameter Estimate Error DF t Value Pr > |t|
equation 1
b01 -1.1132 0.9807 200 -1.14 0.2577
b11 0.05335 0.01887 200 2.83 0.0052
b12 -0.7005 0.3591 200 -1.95 0.0525
equation 2
b02 -3.9737 0.9294 200 -4.28 <.0001
b21 0.05963 0.01636 200 3.64 0.0003
b22 -0.2037 0.3232 200 -0.63 0.5292
Multinomial logistic regression
For the multinomial logistic regression prog (program type) is the response variable. If there are k level for the response variable there will be k-1 equations in the multinomial logistic regression model. Since there are three levels of prog there will be two equations in our model. The first equation is for prog=1 and the second equation for prog=3. Thus, prog=2 is our reference or base category. Each of the equations will have an intercept and coefficients for read and female. The two intercepts are b01 and b03 while the coefficients for read are b11 and b13 and for female the coefficients are b12 and b32.
proc nlmixed data='D:datahsbdemo.sas7bdat';
parms b01=0 b03=0 b11=0 b31=0 b12=0 b32=0;
xb1=b01+b11*read+b12*female;
xb3=b03+b31*read+b32*female;
expxb1=exp(xb1);
expxb2=1;
expxb3=exp(xb3);
den = expxb1+expxb2+expxb3;
if prog=1 then liklhd=expxb1/den;
if prog=2 then liklhd=expxb2/den;
if prog=3 then liklhd=expxb3/den;
ll = log(liklhd);
model prog ~ general(ll);
run;
[partial output]
Parameter Estimates
Standard
Parameter Estimate Error DF t Value Pr > |t| Alpha Lower Upper Gradient
b01 3.0193 1.0963 200 2.75 0.0064 0.05 0.8575 5.1810 3.695E-6
b03 5.3382 1.1564 200 4.62 <.0001 0.05 3.0579 7.6185 -5.42E-6
b11 -0.07121 0.02021 200 -3.52 0.0005 0.05 -0.1111 -0.03137 0.000099
b31 -0.1173 0.02244 200 -5.23 <.0001 0.05 -0.1615 -0.07304 -0.00024
b12 -0.1835 0.3727 200 -0.49 0.6230 0.05 -0.9184 0.5514 5.62E-6
b32 -0.1938 0.3801 200 -0.51 0.6106 0.05 -0.9433 0.5556 -6.56E-7
As with the generalized ordered logistic regression above you will often see the output for multinomial logistic regression from a traditional statistical procedure organized by groups similar to what is shown below.
Standard Parameter Estimate Error DF t Value Pr > |t| prog=1 b01 3.0193 1.0963 200 2.75 0.0064 b11 -0.07121 0.02021 200 -3.52 0.0005 b12 -0.1835 0.3727 200 -0.49 0.6230 prog=3 b03 5.3382 1.1564 200 4.62 <.0001 b31 -0.1173 0.02244 200 -5.23 <.0001 b32 -0.1938 0.3801 200 -0.51 0.6106 prog=2 is the base category
Poisson regression
Our next example is a count model using a poisson distribution. The response variable for this model is awards which is a count of the number of awards received by a student during high school. As with the binary logistic model above there are two ways you can parameterize this model. The first method will use the poisson distribution option.
The poisson regression model uses awards
proc nlmixed data='D:datahsbdemo.sas7bdat';
parms b0=0 b1=0 b2=0;
xb=b0+b1*read+b2*female;
mu = exp(xb);
model awards ~ poisson(mu);
run;
[partial output]
Parameter Estimates
Standard
Parameter Estimate Error DF t Value Pr > |t| Alpha Lower Upper Gradient
b0 -3.0922 0.3325 200 -9.30 <.0001 0.05 -3.7479 -2.4364 -0.00012
b1 0.06009 0.005404 200 11.12 <.0001 0.05 0.04944 0.07075 -0.01074
b2 0.4690 0.1142 200 4.11 <.0001 0.05 0.2438 0.6941 -0.00015
The second method for parameterizing this model is to compute the log-likelihood and use the general log-likelihood function. The two approached yield the same results.
proc nlmixed data='D:datahsbdemo.sas7bdat';
parms b0=0 b1=0 b2=0;
xb=b0+b1*read+b2*female;
mu = exp(xb);
ll = awards*log(mu) - mu - lgamma(awards+1);
model awards ~ general(ll);
run;
[partial output]
Parameter Estimates
Standard
Parameter Estimate Error DF t Value Pr > |t| Alpha Lower Upper Gradient
b0 -3.0922 0.3325 200 -9.30 <.0001 0.05 -3.7479 -2.4364 -0.00012
b1 0.06009 0.005404 200 11.12 <.0001 0.05 0.04944 0.07075 -0.01074
b2 0.4690 0.1142 200 4.11 <.0001 0.05 0.2438 0.6941 -0.00015
Negative binomial regression
Our final example is a negative binomial regression. We will use the same response variable, awards, as was used in the poisson example. One way of conceptualizing a negative binomial model is to think of it as a poisson model with overdispersion, that is, excess variance. In a true poisson model the mean and the variance are equal. However, variables distributed as a negative binomial have a variance that is greater than the mean. The parameterization if this model is similar to that of the poisson with the addition of an addition parameter, alpha, which measures the degree of over dispersion.
proc nlmixed data='D:datahsbdemo.sas7bdat';
parms b0=0 b1=0 b2=0 alpha=.1;
xb=b0+b1*read+b2*female;
mu = exp(xb);
m = 1/alpha;
ll = lgamma(awards+m)-lgamma(awards+1)-lgamma(m)
+ awards*log(alpha*mu)-(awards+m)*log(1+alpha*mu);
model awards ~ general(ll);
run;
[partial output]
Parameter Estimates
Standard
Parameter Estimate Error DF t Value Pr > |t| Alpha Lower Upper Gradient
b0 -3.3266 0.4111 200 -8.09 <.0001 0.05 -4.1372 -2.5160 0.000113
b1 0.06392 0.006809 200 9.39 <.0001 0.05 0.05050 0.07735 0.010318
b2 0.5065 0.1354 200 3.74 0.0002 0.05 0.2396 0.7735 0.000062
alpha 0.1828 0.08455 200 2.16 0.0318 0.05 0.01606 0.3495 -0.00016