BookmarkSubscribeRSS Feed
BjoernHolzhauer
Obsidian | Level 7

Dear all,

I am trying to help convergence in a repated measures model (where in a simulation study I quite often got an infinite likelihood - this is a bit odd, because with 50 and 100 observations this hardly happens, while with 400 it frequently does ~25% of simulations) by providing starting values for a completely unstructured covariance from a simpler covariance structure that can always be fitted. Since my data structure is two different variables (with different variability) at 6 different times, I tried a UN @ CS as a simpler covariance structure.

Running the following code seems to work fine.

ods output covparms=covparms;

proc mixed data=both NOCLPRINT;

  by sim;

  class trt(ref='0') week wk param;

  model value = trt week trt*week bhba1c bhba1c*week bfpg bfpg*week bfpg*bhba1c bfpg*bhba1c*week / ddfm=kr NOTEST;

  repeated param wk / subject=pat type=UN@CS group=trt;

run;

 

I then tried to get the covariance parameters estimated from the simulations

Covariance Parameter Estimates
Cov Parm      Subject  Group  Estimate
param UN(1,1) pat      trt 1  0.5046

UN(2,1)       pat      trt 1  0.9147

UN(2,2)       pat      trt 1  5.0867
wk Corr       pat      trt 1  0.5229
param UN(1,1) pat      trt 0  0.4395

UN(2,1)       pat      trt 0  0.6968
UN(2,2)       pat      trt 0  4.0510
wk Corr       pat      trt 0  0.4104

into the same format that I think is required as input for the covariance parameters for a completely unstuctured covariance matrix and get something like this:

sim trt CovParm Group Est

1    1  UN(1,1) trt 1 0.50457

1    1  UN(2,1) trt 1 0.26385

1    1  UN(2,2) trt 1 0.50457

1    1  UN(3,1) trt 1 0.26385

1    1  UN(3,2) trt 1 0.26385

1    1  UN(3,3) trt 1 0.50457

...

However, when I then run

proc mixed data=both NOCLPRINT;

  by sim;

  class trt(ref='0') week;

  parms / PARMSDATA=covparmsUN2(keep=sim Est trt);

  model value = trt week trt*week bhba1c bhba1c*week bfpg bfpg*week bfpg*bhba1c bfpg*bhba1c*week / ddfm=kr NOTEST;

  repeated week / subject=pat type=UN group=trt;

run;

I get

WARNING: Unable to make hessian positive definite.

NOTE: The above message was for the following BY group:

      sim=1

NOTE: PROCEDURE MIXED used (Total process time):

      real time           4.21 seconds

      cpu time            4.18 seconds

What is going on here? Does estimate covariance matrix for the UN@CS structure not ensure this does not happen? Am I missing something or do I not understand how SAS expects the input*?

Many thanks for any comments,

Björn

*: When I run with type=UN, take the covparms and use them directly as input, this seems to work fine, so I assume this is fine, but perhaps there is something I overlooked.

13 REPLIES 13
SteveDenham
Jade | Level 19

Not real sure here (calling ), but I think the problem may be that the 'constructed' covariance matrix has eigenvalue issues, leading to the Hessian warning.  This is a major concern when constructing simulated MVN data--specifying a particular covariance matrix is an art, and I suspect the same considerations need to be applied when passing starting values to MIXED (or any of the mixed model procedures).

Steve Denham

Rick_SAS
SAS Super FREQ

I'm not an expert on MIXED, butI believe that the parameters in the PARMSDATA= data set need to match the entries in the "Covariance Parameter Estimates" table.  Have you specified all of the parameters in the order in which they appear in the table?

The starting parameters you provide appear to be a Toeplitz structure.  I see no reason why that would be a good initial guess for an unstructured covariance matrix, so maybe these initial guesses do not converge. On the other hand, I don't see why they would be particularly awful guesses.

My recommendation is to include

WHERE sim=1;

and try to understand what is happening with that one sample of simulated data.

BjoernHolzhauer
Obsidian | Level 7

I wrote out the Covparms dataset I get when I do TYPE=UN and ordered my input dataset in the exact same order (same number of records and populated as I think it should be). The order is definitely the right one, because I also tested using TYPE=UN and reading the CovParms dataset back in as starting values and that worked perfectly, so in principle SAS does what the documentation says in this respect.

Actually, I ensure that there was only one simulated dataset (so as it happens the code turns out to do the same as if I had done where sim=1), so the "by sim" is only testing that nothing goes wrong when doing by-group processing.

If SAS does not ensure the UN@CS has a positive definite Hessian, then I guess I could start adding things to the diagnoal until it's positive definite (although to know how much to do to always ensure that may be tricky).

Rick_SAS
SAS Super FREQ

I think that the UN@CS estimate will be PD. (Your example is.) The question is whether the optimization algorithm can get from that initial guess to the TYPE=UN estimate while preserving PD. I don't know that answer. Is there an option to turn on a history of the optimization or to get more details from the procedure? 

lvm
Rhodochrosite | Level 12 lvm
Rhodochrosite | Level 12

Anytime one uses UN, one is not guaranteed of a PD matrix. A warning is given in the Log. For your example, I don't see how you are going from your outputted G matrix of UN@CS to the G matrix your are inputting to the next run.

lvm
Rhodochrosite | Level 12 lvm
Rhodochrosite | Level 12

Also, you said that you have 6 times, but your last G matrix is for 3 times.

BjoernHolzhauer
Obsidian | Level 7

I only gave an excerpt of the dataset in the previous post to keep it short. The full version would be (columns: sim, trt, CovParm, group, estimate)

The full dataset that I created is below, which is the same ordering I get when I write out the covparms dataset when fitting a model with type=UN group=trt. The implied correlation matrix is certainly positive definite. The fact that I do not even get an iteration history makes me wonder whether the initial matrix is somehow problematic or ends me up in a situation where PROC MIXED cannot iterate (but then I wonder how that worked with the UN@CS matrix...).

11UN(1,1)trt 10.50457
11UN(2,1)trt 10.26385
11UN(2,2)trt 10.50457
11UN(3,1)trt 10.26385
11UN(3,2)trt 10.26385
11UN(3,3)trt 10.50457
11UN(4,1)trt 10.26385
11UN(4,2)trt 10.26385
11UN(4,3)trt 10.26385
11UN(4,4)trt 10.50457
11UN(5,1)trt 10.26385
11UN(5,2)trt 10.26385
11UN(5,3)trt 10.26385
11UN(5,4)trt 10.26385
11UN(5,5)trt 10.50457
11UN(6,1)trt 10.26385
11UN(6,2)trt 10.26385
11UN(6,3)trt 10.26385
11UN(6,4)trt 10.26385
11UN(6,5)trt 10.26385
11UN(6,6)trt 10.50457
11UN(7,1)trt 10.91469
11UN(7,2)trt 10.47830
11UN(7,3)trt 10.47830
11UN(7,4)trt 10.47830
11UN(7,5)trt 10.47830
11UN(7,6)trt 10.47830
11UN(7,7)trt 15.08672
11UN(8,1)trt 10.47830
11UN(8,2)trt 10.91469
11UN(8,3)trt 10.47830
11UN(8,4)trt 10.47830
11UN(8,5)trt 10.47830
11UN(8,6)trt 10.47830
11UN(8,7)trt 12.65992
11UN(8,8)trt 15.08672
11UN(9,1)trt 10.47830
11UN(9,2)trt 10.47830
11UN(9,3)trt 10.91469
11UN(9,4)trt 10.47830
11UN(9,5)trt 10.47830
11UN(9,6)trt 10.47830
11UN(9,7)trt 12.65992
11UN(9,8)trt 12.65992
11UN(9,9)trt 15.08672
11UN(10,1)trt 10.47830
11UN(10,2)trt 10.47830
11UN(10,3)trt 10.47830
11UN(10,4)trt 10.91469
11UN(10,5)trt 10.47830
11UN(10,6)trt 10.47830
11UN(10,7)trt 12.65992
11UN(10,8)trt 12.65992
11UN(10,9)trt 12.65992
11UN(10,10)trt 15.08672
11UN(11,1)trt 10.47830
11UN(11,2)trt 10.47830
11UN(11,3)trt 10.47830
11UN(11,4)trt 10.47830
11UN(11,5)trt 10.91469
11UN(11,6)trt 10.47830
11UN(11,7)trt 12.65992
11UN(11,8)trt 12.65992
11UN(11,9)trt 12.65992
11UN(11,10)trt 12.65992
11UN(11,11)trt 15.08672
11UN(12,1)trt 10.47830
11UN(12,2)trt 10.47830
11UN(12,3)trt 10.47830
11UN(12,4)trt 10.47830
11UN(12,5)trt 10.47830
11UN(12,6)trt 10.91469
11UN(12,7)trt 12.65992
11UN(12,8)trt 12.65992
11UN(12,9)trt 12.65992
11UN(12,10)trt 12.65992
11UN(12,11)trt 12.65992
11UN(12,12)trt 15.08672
10UN(1,1)trt 00.43945
10UN(2,1)trt 00.18037
10UN(2,2)trt 00.43945
10UN(3,1)trt 00.18037
10UN(3,2)trt 00.18037
10UN(3,3)trt 00.43945
10UN(4,1)trt 00.18037
10UN(4,2)trt 00.18037
10UN(4,3)trt 00.18037
10UN(4,4)trt 00.43945
10UN(5,1)trt 00.18037
10UN(5,2)trt 00.18037
10UN(5,3)trt 00.18037
10UN(5,4)trt 00.18037
10UN(5,5)trt 00.43945
10UN(6,1)trt 00.18037
10UN(6,2)trt 00.18037
10UN(6,3)trt 00.18037
10UN(6,4)trt 00.18037
10UN(6,5)trt 00.18037
10UN(6,6)trt 00.43945
10UN(7,1)trt 00.69679
10UN(7,2)trt 00.28599
10UN(7,3)trt 00.28599
10UN(7,4)trt 00.28599
10UN(7,5)trt 00.28599
10UN(7,6)trt 00.28599
10UN(7,7)trt 04.05102
10UN(8,1)trt 00.28599
10UN(8,2)trt 00.69679
10UN(8,3)trt 00.28599
10UN(8,4)trt 00.28599
10UN(8,5)trt 00.28599
10UN(8,6)trt 00.28599
10UN(8,7)trt 01.66271
10UN(8,8)trt 04.05102
10UN(9,1)trt 00.28599
10UN(9,2)trt 00.28599
10UN(9,3)trt 00.69679
10UN(9,4)trt 00.28599
10UN(9,5)trt 00.28599
10UN(9,6)trt 00.28599
10UN(9,7)trt 01.66271
10UN(9,8)trt 01.66271
10UN(9,9)trt 04.05102
10UN(10,1)trt 00.28599
10UN(10,2)trt 00.28599
10UN(10,3)trt 00.28599
10UN(10,4)trt 00.69679
10UN(10,5)trt 00.28599
10UN(10,6)trt 00.28599
10UN(10,7)trt 01.66271
10UN(10,8)trt 01.66271
10UN(10,9)trt 01.66271
10UN(10,10)trt 04.05102
10UN(11,1)trt 00.28599
10UN(11,2)trt 00.28599
10UN(11,3)trt 00.28599
10UN(11,4)trt 00.28599
10UN(11,5)trt 00.69679
10UN(11,6)trt 00.28599
10UN(11,7)trt 01.66271
10UN(11,8)trt 01.66271
10UN(11,9)trt 01.66271
10UN(11,10)trt 01.66271
10UN(11,11)trt 04.05102
10UN(12,1)trt 00.28599
10UN(12,2)trt 00.28599
10UN(12,3)trt 00.28599
10UN(12,4)trt 00.28599
10UN(12,5)trt 00.28599
10UN(12,6)trt 00.69679
10UN(12,7)trt 01.66271
10UN(12,8)trt 01.66271
10UN(12,9)trt 01.66271
10UN(12,10)trt 01.66271
10UN(12,11)trt 01.66271
10UN(12,12)trt 04.05102
SteveDenham
Jade | Level 19

Wow.  That is a lot of parameters (132 if I am not mistaken).  Now I am really curious as to how you were able to get anything at all to happen with 50 or 100 observations.  With only 400 simulated values, I am not surprised that you get into problems.  Have you considered a semi-parametric approach to the repeated measures, fitting a spline structure, or even some sort of autoregressive structure?  No guarantees, but it is the first thing that popped into my mind.

Steve Denham

BjoernHolzhauer
Obsidian | Level 7

Actually it's 50 patients x 12 time points x 2 treatment groups = 1200 observations (of course not independent).

In fact we had the most problems with 400/arm (with non-positive definite Hessians, but possibly more a numeric problem?!!), and hardly any problems with 50 or 100/arm. We were trying to get out of that by giving MIXED better starting values...

SteveDenham
Jade | Level 19

<sheepish grin>.  Should have done the multiplication myself--sorry about the first statement.

So about a quarter of the simulations at 400 give NPD problems, and you don't see a problem, even occasionally (say 1 or 2) at n=50 or n=100.  Then I do worry about the numeric situation, and how the data were simulated.  Were they resampled from known data, or was a MVN simulation done?  If the latter, was any structure applied to mimic the repeated nature?  I would start looking for "pathology" in the simulated datasets that blow up.  I realize this is now way off-topic from the subject of how to get good starting values in, but I guess I am looking for a different kind of solution.

Steve Denham

lvm
Rhodochrosite | Level 12 lvm
Rhodochrosite | Level 12

The thought of using an unstructured 12x12 matrix for a repeated measure is intimidating. You have a awful lot of parameters, even for 1200 observations. Did you try method=mivque0, just to see what the G matrix looks like with a noniterative method (if this works)? By the way, the warning you get about the Hession NPD is not the same as G being NPD.

http://support.sas.com/resources/papers/proceedings12/332-2012.pdf

BjoernHolzhauer
Obsidian | Level 7

Yes, the 12x12 matrix is big, so we were doing a simulation study to compare that approach with ones with simpler covariance matrices and other alternative approaches. We had not tried method=mivque0. I realise that the Hessian being npd is not the same as the covariance matrix being npd, but that could have been one source (and the one that is the most easily checked).

The whole idea with starting values came up as an idea to get out of the problem that in about 25% of our simulations with 400 patients/arm we had an "infinite likelihood" error (but not with 50 or 100 patients per arm), which does seem rather illogical. Our guess as to what went wrong was that the likelihood ended up in bad regions due to poor covariance parameters (possibly due to poor starting values)... The error message we got is shown below.

NOTE: Convergence criteria met.
NOTE: The above message was for the following BY group:
      sim=1
WARNING: Stopped because of infinite likelihood.
NOTE: The above message was for the following BY group:
      sim=2
NOTE: Convergence criteria met.
NOTE: The above message was for the following BY group:
      sim=3
NOTE: The data set WORK.JMMRM has 2 observations and 14 variables.
43                                                         The SAS System                        04:42 Wednesday, September 10, 2014

NOTE: PROCEDURE MIXED used (Total process time):
      real time           13:40.78
      cpu time            13:41.20

SteveDenham
Jade | Level 19

So for SIM=2, I am going to assume that it never even got started, so there is no iteration history to look at.  Is there any possibility at all of duplicate subjects, or even more weird, exactly duplicate records except for subject_id?  I know I'm spending way more time on the pathology than on how to get good starting values, but...

Use the mivque0 values as starting values.  That should work unless you have a negative on the diagonal--and you well might have, given the nature of simulated data.

Consider changing to HPMIXED or GLIMMIX, where you can specify type=CHOL for the unstructured data.  Should get away from the NPD Hessian problem, if the data allow.

Steve Denham

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
  • 13 replies
  • 4012 views
  • 1 like
  • 4 in conversation