Likelihood in PROC NLMIXED with the warning of Hessian matrix not positive definite_Revised version

[bienco88]

Hi,

While running my non-linear mixed effects model in SAS for 1000 simulated data sets, I sometimes faced the warning of:

"The convergence criterion (GCONV = 1E-8) satisfied.

At least one element of the gradient is greater than 1e-3.

Moore-Penrose inverse is used in covariance matrix.

The final Hessian matrix is not positive definite, and therefore the estimated covariance matrix is not full rank and may be unreliable.  The variance of some parameter estimates is zero or some parameters are linearly related to other parameters."

Regarding the model specification: I have the data set based on the simulation. So the model specification is the model I have used in Simulx to generate the data. So I know exactly how the model should be. Still, the problem came out.

My question is: I only interest in the likelihood of the obtained model. So in this case, can I just ignore this warning and go on with the value of the -2log likelihood.

If yes, could you back up me with a good reference. I would like to read more about this.

I am using SAS 9.4

Thank you a lot.

P/s: Here is my example code:

In the attachment, I provide some portions of my data.

/*Random A0 - Null model - No censoring - Formulate to calculate the sd for A0*/
 proc nlmixed data = sim1 xtol = 1E-12 method = GAUSS qpoints = 100;

   by replicate;
   parms A0 =13 beta = 0.6 gamma = 0.6 alpha = 0.4 h = 6.27
         sigma_A0 = 2.85 betaA0 = 5.01
         sigma = 0.4;
   bounds A0 beta gamma alpha > 0;
   pi = constant('pi');
   mu1 = (A0+A0_i+betaA0*Kind2)*exp(-beta*t);
   mu2 = (A0+A0_i+betaA0*Kind2)*exp(-beta*2)*exp((t-2)*gamma);
   mu3 = (A0+A0_i+betaA0*Kind2)*exp(-beta*2)*exp((h-2)*gamma)*exp(-alpha*(t-h));
   sigma2 = sigma*sigma;
   if t<= 2  then
     ll = (1/(sqrt(2*pi*sigma2)))*exp(-(y - log10(mu1))**2/ (2*sigma2));
   if t>2 and t <= h  then
     ll = (1/(sqrt(2*pi*sigma2)))*exp(-(y - log10(mu2))**2/ (2*sigma2));
   if t > h  then
     ll = (1/(sqrt(2*pi*sigma2)))*exp(-(y - log10(mu3))**2/ (2*sigma2));
   L = log(ll); /*Log-likelihood*/
   model y ~ general(L);
   random A0_i ~ normal(0,sigma_A0*sigma_A0)
          subject = ID;

  ods output FitStatistics = Fit_Out.csv;
run; quit;  /*95% CI of all parameters contain the true values*/

[Rick_SAS]

No. It is generally a bad idea to use the last result of an optimization that did not converge. Try turning on the ITDETAILS option and using a WHERE claus to restrict to one of the BY-groups that fails. That might give you insight as to why the computation is failing. Is the iteration leading off to inifity? Maybe sigma keeps getting closer and closer to 0 (which can happen for flat data)? Are your initial estimates poor for the data that fail?

It might not make a difference, but consider rewriting the computations in the code directly in terms of the log-likelihood and the log-mu parameter.  Currently you are computing quantities and then taking the LOG of the result. It seems like the computation will be more stable if you compute the log-quantities.  (Translation: you can explore larger regions of the parameter space.)

Along those lines, use LOG(mu) instead of LOG10(mu) so that you can compute logMu directly.

[bienco88]

Hi,

Thank you for your quick response.

I see also in the log: "Convergence criterion (GCONV = 1E-8) satisfied". Does it mean that  the model already converged?

Other things related to the model specification: I find it strange since I used the simulated data. So I have 1000 simulated data sets which basicly are based on the model that I specified in the code (I use Simulx so I can generate the data based on a predefined model).

Thao.

[Rick_SAS · Accepted solution]

That doesn't surprise me. Depending on the sample size and the magnitude of the noise, the simulated data might not look much like the model.

See the article "Monitor convergence during simulation studies in SAS" for a discussion and several references. The example is logistic regression, but the general ideas are applicable to your situation.

I don't know whether or not the convergence note is applicable. You get one note for each BO group. Look in the log to see if that note is for the same BY group as the WARNING. It could be from the previous BY group.

I notice that you didn't include SIGMA or SIMA_A0 on the BOUNDS statement. I assume that these shape parameters must be positive.

[bienco88]

Hi,
Thank you.

I will look into the paper you refered.

1. Regarding to your comment of sigma and sigma_A0: Since they are standard error so I did not put them into the bound statement.

I model sigma, and the variance is sigma2 = sigma*sigma, the same thing for sigma_A0.

2. Regarding the message in SAS: I will recheck.

Thao.

[bienco88]

Yes,

I have checked again.
Here is the whole message for a replicate:

"NOTE: Convergence criterion (GCONV=1E-8) satisfied.
NOTE: Moore-Penrose inverse is used in covariance matrix.
WARNING: The final Hessian matrix is not positive definite, and therefore the estimated
         covariance matrix is not full rank and may be unreliable.  The variance of some
         parameter estimates is zero or some parameters are linearly related to other
         parameters."

Is the convergence satisifed indicating the convergence of the model already?

In this case, can I use the log-likelihood then?

Thank you for your input.
Thao.

[bienco88]

Hi Rick_SAS,

1. I have seen that for some simulated data set, reparameterization helps.

For example, instead of modeling the variance of the random effect of A0 as: sigma_A0*sigma_A0, instead I used: exp(2*logsigmaA0), so A0_i ~ normal(0, exp(2*logsigmaA0)), then the convergence is fine, there is also no warning message and the standard errors are well achieved.

2. Of course, it does not work for all data sets.

I have tried to re-write the code, where I directly estimated the variance of the random effect of A0: A0_i ~ normal(0, sigma2_A0) and then put sigma2_A0 >=0 into a bound statement, now there is another warning:

"

NOTE: FCONV convergence criterion satisfied.

NOTE: At least one element of the (projected) gradient is greater than 1e-3.

WARNING: The final Hessian matrix is full rank but has at least one negative eigenvalue.

         Second-order optimality condition violated.

"

3. I am thinking of other re-parameterization for the variance of the random effect as it seems that this causes the problem, but have not come up with any positive solution now.

Thanks.