BookmarkSubscribeRSS Feed
tbanh
Fluorite | Level 6

Hi all,

 

I'm fitting a nonlinear mixed-effects model in NLMIXED to a set of longitudinal birth weight data. Because the data is strictly positive, I wanted to model the response with a log-normal distribution. As you may be aware, the log-normal distribution is not a default option in the MODEL statement. So I went ahead and computed the log-likelihood and used the GENERAL statement. Below is the code I used to fit the model:

 

proc nlmixed data = weight gconv = 0;
	parms b0 = 1 b1 = .6 b2 = 1 var0 = .02 cov10 = -.01 var1 = .02 cov20 = -.001 cov21 = -.04 var2 = .4 s2e = .001;
	beta0 = b0 + u0; beta1 = b1 + u1; beta2 = b2 + u2; 
	predv = beta0 + beta1*(1-exp(-beta2*Years));
	pi = arcos(-1);

	ll = (-1/2)*((log(2*pi)) + (log(s2e)) + ((log(Years - predv)**2)/s2e)) + Years;
	model Weightkg ~ general(ll);
	random u0 u1 u2 ~ normal([0,0,0], [var0,cov10,var1,cov20,cov21,var2]) subject = ID;
run;

When I run the code, this is the error code I receive:

NOTE: Execution error for observation 1

 

I've looked around and a possible reason for this error message is that my starting values may be off. I used starting values from fitting the model using Bayesian estimation with PROC MCMC but to no avail. Has anyone ever dealt with this issue and can offer any solutions? Thanks!

9 REPLIES 9
JacobSimonsen
Barite | Level 11

I think you have not specified your loglikelihood function correctly. For example, the sign for the quadratic term should be negative. An easier way would be to make a variable inside NLMIXED which is log to your observations. Then use the built-in likelihood for normal distributed data. The two ways to (using log-normal with generel loglikelood, and built-in normal likelihood) is equivalent as I can illustrate with this example

data simulation;
  do i=1 to 100;
    y=exp(rand('normal',0,1));
	output;
  end;
run;

*method 1;
proc nlmixed data=simulation;
  parm mu 0.5;
  z=log(y);
  model  z ~ normal (mu,1);
run;

*method 2;
proc nlmixed data=simulation;
  parm mu 0;	
  pi = arcos(-1);
  s=1;
  ll = -log(2*pi)/2 - log(s) - ((log(y) - mu)/(2*s))**2 - log(y);
  model  y ~ general (ll);
run;
SteveDenham
Jade | Level 19

Hi Jacob,

 

Well, they both converge to the same point estimate, but the error estimates and log likelihood estimates don't match up.  My brain is sort of fried today, so if you have any ideas why the error estimates differ so drastically...

 

See, I would pick method 2 based on information criteria, but would end up trading off for wider confidence intervals.  That makes me think there may be something missing somewhere.

 

Steve Denham

JacobSimonsen
Barite | Level 11

I did a mistake, sorry. I scale with 1/2 within the squared term. That should be done outside.


proc nlmixed data=simulation;
  parm mu 0;	
  pi = arcos(-1);
  s=1;
  ll = -log(2*pi)/2 - log(s) - ((log(y) - mu)/(s))**2/2 - log(y);
  model  y ~ general (ll);
run;
JacobSimonsen
Barite | Level 11
Such mistake just emphasize why its better to use the built-in functionality whenever possible:-)
SteveDenham
Jade | Level 19

Whcih is why I don't use NLMIXED as much as I should.

 

Anybody have a spare Gompertz likelihood lying around...?

 

Steve Denham

tbanh
Fluorite | Level 6

Hi Jacob,

 

Thanks for the reply. I tried your code but I still received this message:

 

NOTE: Execution error for observation 1

 

It seems that the likelihood isn't the problem?

JacobSimonsen
Barite | Level 11

Show the code you run.

I am confused why "years" is used in the defintion of the predictor (predv), but in the likelihood function it seems that years is your observations you try to predict.

tbanh
Fluorite | Level 6

Hi Jacob,

 

Thanks for pointing that egregious error. Yes, I mistakenly put Years as my dependent variable when I should have put Weight. The likelihood function works now.

 

However, I'm running into a new problem. The model I'm running will not converge because I'm providing poor starting values. Do you have any advice on how to fix this?

 

JacobSimonsen
Barite | Level 11

It can be quite difficult to find good starting values. I can not come up with a general good answer, but a few suggestions:

 

Estimate first the parameters in a fixed effect model. Then take these estimate as starting values for the parameters to the mean in the random effect model.

 

Also, a baysian version of the random effect model can give you some good starting values. You can use proc mcmc almost as you have made proc nlmixed. A mcmc method is often very robust to starting values. The mean values (or max values) from the posterior distribution can then afterwards be used as starting values in proc nlmixed.

 

It can also help if you make it such that the parameters is expected to have same order of magnitude. Otherwise there can easily occur convergence problems.

 

sas-innovate-2024.png

Join us for SAS Innovate April 16-19 at the Aria in Las Vegas. Bring the team and save big with our group pricing for a limited time only.

Pre-conference courses and tutorials are filling up fast and are always a sellout. Register today to reserve your seat.

 

Register now!

What is ANOVA?

ANOVA, or Analysis Of Variance, is used to compare the averages or means of two or more populations to better understand how they differ. Watch this tutorial for more.

Find more tutorials on the SAS Users YouTube channel.

Discussion stats
  • 9 replies
  • 2255 views
  • 5 likes
  • 3 in conversation