Random effect estimates (Proc nlmixed)

[Yeganeh]

Dear all

I am trying to run a proc nlmixed program. I have to estimate random effects to assess goodness of fit.

1. How can I estimate random effects for each individual in proc nlmixed?

2. Is it possible to use empirical bayes approach for this purpose in sas?

[Rick_SAS]

I assume you have repeated measurements? The Getting Started example in the doc shows how to estimate random effects for a bunch of trees, where each tree has seven measurements. The basic syntax is 

RANDOM x ~ normal(0,s2) SUBJECT=SubjectID;

where x is the random effect, which is assumed to be normally distributed with variance s2. The parameter s2 is included in the PARMS statement and is estimated from the data.

> 2. Is it possible to use empirical bayes approach for this purpose in sas?

Can you explain what you are trying to accomplish? It would be helpful to see some code, some data, or both.

[Yeganeh]

I have repeated measurements of count data and want to show violation of homogeneous random effects in sex or other known collected factors in the data. I have no code for it

[Rick_SAS]

If you post some sample data and the model you are trying to fit, we might be able to help you construct the PROC NLMIXED code. I might be misunderstanding your statement, but sex is usually treated as a fixed effect.  

[Yeganeh]

Proc nlmixed  maxiter=10000;

Parms b0=1 , b1=0 , b2=0 , b3=0 , a0= 1, g0=1 , V_u1=1;

Aialpha=a0;

Gigamma=g0;

Bibeta= b0 + b1*x1 + b2 *x2 + b3*time;

Lambda=exp(b0 + b1*x1 + b2* x2 + b3*time + u1);

P1= exp(Aialpha) / (1+ exp(Aialpha) + exp (Gigamma));

P2= exp(Gigamma) / (1+ exp(Gigamma) + exp (Aialpha));

P3=1-p1-p2;

If y=0 then ll = log (p1 + p3*exp(-lambda));

If y=1 then ll= log (p2 + p3*(1)*exp(-lambda)* lambda);

If y>1 then ll= log(p3) – lambda + y*log(lambda) – lgamma (y+1);

Model y ~ general (ll);

Random u1 ~ normal (0, [v_u1]) subject=id;

Run;

I need to estimate u1 for each individual. then to assess goodness of fit, I want to draw qq plot of the random effects in the Poisson and inflated parts of the model by gender group. In a well-fitted model, the distribution of estimated random effects (in each model part) should lie along the line y=x. I couldn't find a code for this problem, Could you help me construct a code for it, please?

Thank you

[jiltao]

I do not think it is possible to estimate separate variance for each subject.

If you want to estimate different variance for each gender (valued 0/1), you might do something like the following --

Proc nlmixed  maxiter=10000;

Parms b0=1 , b1=0 , b2=0 , b3=0 , a0= 1, g0=1 , V_u0=1, V_u1=1;

Aialpha=a0;

Gigamma=g0;

Bibeta= b0 + b1*x1 + b2 *x2 + b3*time;

Lambda=exp(b0 + b1*x1 + b2* x2 + b3*time + +u0+u1*gender);

P1= exp(Aialpha) / (1+ exp(Aialpha) + exp (Gigamma));

P2= exp(Gigamma) / (1+ exp(Gigamma) + exp (Aialpha));

P3=1-p1-p2;

If y=0 then ll = log (p1 + p3*exp(-lambda));

If y=1 then ll= log (p2 + p3*(1)*exp(-lambda)* lambda);

If y>1 then ll= log(p3) – lambda + y*log(lambda) – lgamma (y+1);

Model y ~ general (ll);

Random u0 u1 ~ normal ([0,0], [v_u0, 0, v_u1]) subject=id;

Run;

The random effects estimates are by definition the empirical Bayes estimates. You can use the OUT= option in the RANDOM statement to get these estimates in the OUT= data set.

Hope this helps,

Jill

[Yeganeh]

Thanks Jill for the detailed reply.

Thank you for the all answers.

[Rick_SAS · Accepted solution]

Thank you for posting your code. On the RANDOM statement, you can use the OUT= option to output the estimates of the random effects for each subject. The code will look like

Random u1 ~ normal (0, v_u1) subject=id OUT=RandomOut;

Your code doesn't include any information about gender, but it sounds like you want to create two QQ plots, one for males and one for females. You will need to merge the RandomOut data set (which contains the estimates for each ID) and data that says which IDs are male and which are females. Then create the QQ plot for each.

Suppose WANT is the name of the data set that merges the RandomOut estimates and the genders. Then you get the QQ plots from PROC UNIVARIATE by running:

proc sort data=Want; by Sex; run;

proc univariate data=Want(where=(Effect='u1'));
by Sex;
var Estimate;     /* the estimate for u1 for each subject */
qqplot Estimate;
run;

That will create a data set that contains the estimates. You can then