NLMIXED: difference between convergence status and non-positive Hessian matrix

[ptadgerv]

Dear SAS Comunity

Running some nonlinear models, sometimes happen to me that I can have a successful "convergence" message from SAS (Convergestatus table: F Conv), and at the same time  I can have in the Parameter Estimates: SE with infinity value and the message in the Log: 

"The final Hessian matrix is full rank but has at least one negative eigenvalue. Second-order optimality condition violated."

I have 2 questions:

1) Could this behavior (Non-positive definite H matrix, and infinity SE) happen with a G-Convergence?

2) If the previous answer is YES (in the case of NO, the response to the following question is clear now). Is it possible to have the Hessian matrix status log message as a table output?

Thanks in advance!!!!

[ptadgerv]

I would like to rewrite my question, maybe that could help that somebody help me:

I would like to ask how to avoid negative eigenvalues in the Hessian matrix (which as I understand is related to non maximum points or unstable estimation that cannot be trusted)

I mean, how to have in the automatic way when you are running a lot of models, that one of the models had a Hessian matrix with at least one negative eigenvalue.

Usually, the message is presented in the log area, but I'm aware to have it in a table output (as for example convergence status)

[Rick_SAS]

1. Yes, this can happen with G convergence. For possible reasons and references to more information, see "Convergence in mixed models: When the estimated G matrix is not positive definite."

2. In this part of a simulation study? You can monitor the convergence during a simulation study and exclude (or otherwise handle) samples that did not converge. For details, see "Monitor convergence during simulation studies in SAS."

[ptadgerv]

Thanks, Rick for your quick answer. I did read your articles, but I don't find an answer to my question.

My study is not a simulation study but is a study that requires to fit 30 models per data (again and again!), in different scenarios, which makes it almost impossible to check by hand the non-positive definite matrix log message. The models are nonlinear ones, so I'm using NLMIXED.

For sure PROC MIXED (also GLIMMIX if I remember well) is very helpful with the ConvergeStatus

but the "same" table ConvergenceStatus in NLMIXED, only says if you get convergence and which type (nothing about non-positive definite status, like PORC MIXED):

but nothing about the positive definite matrix status :'(

[Rick_SAS]

You are correct. I don't have a simple solution, but you could try something like this to monitor each model:

%macro MonitorErrors(ModelID);
   data ErrorMsg;
      length msg $200;
      ModelID = &ModelID;
      SysErr = &syserr;
      %if &syserr = 4 %then %do;
         msg = "&syswarningtext";
      %end;
      %else %do;
         msg = "OK";
      %end;
   run;
   proc append base=Errors data=ErrorMsg force; run;
%mend;


proc datasets nolist; 
  delete Errors;   /* delete this data set if it exists */
run;

proc nlmixed;
   /* FIRST MODEL */
quit;
%MonitorErrors(1);

proc nlmixed;
   /* SECOND MODEL */
quit;
%MonitorErrors(2);

proc print data=Errors; run;

The basic idea is that if NLMIXED experiences a convergence issue, it will set the SYSERR macro to 4 and store the warning message in the SYSWARNINGTEXT macro. You can write the "model number" and the text to a data set and append each result to the one that comes before it.

I haven't tested the program, so it might need some tweaking. Actually, I've never done this before, so please let me know if it works! 

[ptadgerv]

Thank you Rick!

It's a really good idea but doesn't capture the log message of the non-positive definite H matrix. In different situations (with and without error), the results is the same:

The other idea that I have is to check for an infinitive value in SD for Parm Estimates, or empty values in the variance-covariance matrix. Which almost I saw that appear at the same time with the non-positive definite behavior. 

[ptadgerv]

Fe errata:
Which almost appears at the same time as the non-positive definite H matrix behavior.

[Rick_SAS]

Please post data and some sample SAS code so that I can produce the kinds of errors you are seeing and look for a solution.

[ptadgerv]

data sham;
input confirmed 8. time 8.;
datalines;
38 1
72 2
125 3
206 4
316 5
343 6
407 7
501 8
600 9
774 10
1024 11
1362 12
1541 13
1755 14
2142 15
2564 16
3098 17
3811 18
4473 19
4942 20
5428 21
6756 22
7951 23
9150 24
10513 25
12031 26
12875 27
13558 28
15296 29
16977 30
18493 31
19971 32
21665 33
22587 34
23252 35
25186 36
26701 37
28299 38
30538 39
32874 40
33903 41
34427 42
34964 43
36524 44
38157 45
39831 46
41225 47
41947 48
42390 49
43666 50
44936 51
45713 52
46689 53
47500 54
47888 55
48093 56
48848 57
49417 58
49939 59
50525 60
50762 61
51048 62
51188 63
51858 64
52404 65
52956 66
53398 67
53881 68
54121 69
54239 70
54715 71
55110 72
55431 73
55736 74
56082 75
56229 76
56310 77
56627 78
;
run;

proc nlmixed data=sham MAXITER=500;
parms s2=100000, maximo=500, alpha=57000 , etas=40;
mu = alpha/( (1+k*exp(-gamma*(time-eta)) ) **(1/k) );
model confirmed ~ normal(mu,s2);
run;

Once again thank you for your time and help

[ptadgerv]

Dear Rick

I would provide a more simple example

This model converges and produces a positive definite H matrix:

data pump;
input y t group;
pump = _n_;
logtstd = log(t) - 2.4564900;
datalines;
5 94.320 1
1 15.720 2
5 62.880 1
14 125.760 1
3 5.240 2
19 31.440 1
1 1.048 2
1 1.048 2
4 2.096 2
22 10.480 2
;

proc nlmixed data=pump;
parms logsig 0 1 /* multiple initial guesses for each parameter */
beta1 -1 to 1 by 0.5 
beta2 -1 to 1 by 0.5
alpha1 1 5 10
alpha2 1 5; /* use BEST=10 option to see top 10 choices */
if (group = 1) then eta = alpha1 + beta1*logtstd + e;
else eta = alpha2 + beta2*logtstd + e;
lambda = exp(eta);
model y ~ poisson(lambda);
random e ~ normal(0,exp(2*logsig)) subject=pump;
run;

The following model converges but produces a non-positive definite H matrix:

data pump;
input y t group;
pump = _n_;
logtstd = log(t) - 2.4564900;
datalines;
5 94.320 1
1 15.720 1
5 62.880 1
14 125.760 1
3 5.240 1
19 31.440 1
1 1.048 1
1 1.048 1
4 2.096 1
22 10.480 2
;

proc nlmixed data=pump;
parms logsig 0 1 /* multiple initial guesses for each parameter */
beta1 1
beta2 2
alpha1 10
alpha2 1 5; /* use BEST=10 option to see top 10 choices */
if (group = 1) then eta = alpha1 + beta1*logtstd + e;
else eta = alpha2 + beta2*logtstd + e;
lambda = exp(eta);
model y ~ poisson(lambda);
random e ~ normal(0,exp(2*logsig)) subject=pump;
run;

These two are my most common scenarios...

[ptadgerv]

Dear Rick
I think that the problem with your previous code, it's that we need to delete the memory of the variable! (I don't know why is not possible to overwrite)
It's obvious that after each procedure these variables:

%put &syserr;
%put &syswarningtext;

Capture perfectly what is happening in the log area.

But the variable needs to be deleted to capture again the log.

For example, I tried it with this (between one nlmixed and the next one), and it works perfectly (but not inside your macro):

%macro deleteALL;
   	options nonotes;
  	%local vars;
  	proc sql noprint;
      	     select name into: vars separated by ' '
         	  from dictionary.macros
            	      where scope='GLOBAL' 
			   and not name contains 'SYS_SQL_IP_';
   	quit;
   	%symdel &vars;
   	options notes;
    	%put NOTE: Macro variables deleted.;
%mend deleteALL;
%deleteALL;

[Rick_SAS · Accepted solution]

Thanks for posting data and code. Very helpful. I think the following will work. Try it out.

%macro MonitorErrors(ModelID, ErrorNum, Text);
   data ErrorMsg;
      length msg $200;
      ModelID = &ModelID;
      SysErr = &ErrorNum;
      %if &ErrorNum = 0 %then %do;
         msg = "OK";
      %end;
      %else %do;
         msg = &Text;
      %end;
   run;
   proc append base=Errors data=ErrorMsg force; run;
%mend;

proc datasets nolist; 
  delete Errors;   /* delete this data set if it exists */
run;

proc nlmixed ;
/* first model */
run;
%MonitorErrors(1, &SysErr, "&SysWarningText");

proc nlmixed;
/* second model */
run;
%MonitorErrors(2, &SysErr, "&SysWarningText");

proc print data=Errors; run;

[ptadgerv]

Thank you Rick for such an important and necessary approach in nlmixed!. 
If you ask me, this has to be the default behavior of statusconverge table in nlmixed, to inform about the non-positive definite H matrix (as MIXED and GLIMMIX does).

[ptadgerv]

Dear Rick

Is it possible to adapt your macro into a proc nlmixed with the "by statement"? I tried, but the macro didn't record all the log messages, just the last one.

[Rick_SAS]

No, SYSERR only returns the results from the last BY group, I believe. That is why I use the ConvergenceStatus table to monitor convergence in simulation studies.