How to incorporate correlation between random error (epsilon) in proc nlmixed joint MMRM model?

[lcw68]

Hi, if I have a correlated longitudinal model with one continuous outcome and one ordinal outcome (y2it*, using probit normal to define category from latent normal variale y2it), say y1it = Xalphai + b1i + epsilon_{1it}, y2it (latent normal) = Xbetai + b2i+epsilon_{2it}, (b1i,b2i) are bivariate normal distribution, and (eps_{1it},eps_{2it})are random error that also have bivariate normal distribution (correlated). (the variance parameter is  (sigcon, rhoe*sigcon*sigord, sigord)). How should I specify the "rhoe" ? I attach my SAS code here.

PROC NLMIXED data=auto2 qpoints=10 tech = NMSIMP corr ecorr;
PARMS beta0=6 beta1=-3 beta2=0 alpha0=25
alpha1=-10 alpha2=-1 sigcon = 2 sigor = 3 rhob=0.3 i1 = 2 i2=2 rhoe = 0.4;
bounds i1 > 0, i2 > 0;
if variable="continuous" then do;
pi=constant("pi");
m1=alpha0+alpha1*trt+alpha2*time+z1;
LL1=((-0.5)*log(2*pi*sigcon*sigcon)-((value1-m1)**2)/
(2*sigcon*sigcon));
END;
if variable="ordinal" then do;
eta = beta0+beta1*trt+beta2*time+z2;
if (value1=1) then p = probnorm(-eta/sigor);
else if (value1=2) then
p = probnorm((i1-eta)/sigor) - probnorm(-eta/sigor);
else if (value1=3) then
p = probnorm((i1+i2-eta)/sigor) - probnorm((i1-eta)/sigor);
else p = 1 - probnorm((i1+i2-eta)/sigor);
if (p > 1e-8) then LL2 = log(p);
else LL2 = -1e20;
END;
LL=LL1+LL2;
MODEL value1 ~ GENERAL(LL);
RANDOM z1 z2 ~ NORMAL([0,0],[1,rhob,1]) SUBJECT=id;
estimate "sigma2" sigor*sigor;
ods output ParameterEstimates=pars;
ods output FitStatistics=aic;
RUN;

My question is where should I specify rhoe here? I thought it should be similar like type = "un"...

BTW if the error " No valid parameter points were found" happens, is it only because the bad starting values? Thanks!

[sbxkoenk]

Hello,

What is MMRM in your subject line?

• Mixed-Effects Model for Repeated Measures    or
• Marginal Models for Repeated Measurements

I cannot answer your main question (off the top of my head), but about those joint models with PROC NLMIXED you can find a lot about that on the internet.

As for the answer to your side question (proc nlmixed - No valid parameter points were found) :

The error message of
ERROR: No valid parameter points were found.
might indicate that the initial value for the parameters might be unreasonable, or the parameterization of your model might be problematic.

Koen

[lcw68]

Hi, thanks for your reply. It's Mixed-Effects Model for Repeated Measures 

I did find many joint model examples but non of them include correlated random error model. So I don't know how to include it.

Thanks for the response for the second question!

[sbxkoenk]

Hello,

For this type of models I always turn to Geert and Geert (MOLENBERGHS, Geert & VERBEKE, Geert).

They are also compatriots of mine.

See this article from 2015 (first author is IVANOVA, Anna) :

Mixed model approaches for joint modeling of different types of responses

https://documentserver.uhasselt.be/bitstream/1942/20856/1/463a.pdf

See Appendix B on SAS Implementation.

As an example, the authors present the routines for the case of a joint model for continuous-ordinal responses
with different random effects, some of which are correlated.

Good luck,

Koen

[lcw68]

Thanks very much, I'll read it through

[lcw68]

Hi Koen, I read the paper and it also gives an example that only b (random effect) are correlated, but sigma (std of random error epsilon) are still independent. But thanks for your help! 

[sbxkoenk]

Hello @lcw68 , 

Tomorrow (Monday) is working day again.

I think @jiltao or @StatsMan might be able to help you out.

Koen

[lcw68]

Hello, I found "Error: No valid parameter points were found" still happens.

I tried to run the two marginal model for starting values, the code is here: (autocon contains continuous outcome, autoord contains ordinal outcome)

proc nlmixed data=autocon corr ecorr;
parms alpha1=0.2 alpha2=0.3 sigcon = 1 sd = 1;
pi=constant("pi");
m1=alpha0+alpha1*trt+alpha2*time+z1;
LL1=((-0.5)*log(2*pi*sigcon*sigcon)-((value-m1)**2)/
(2*sigcon*sigcon));
model value ~ general(LL1);
random z1 ~ normal(0,sd*sd) subject=id;
ods output ParameterEstimates=pars;
ods output FitStatistics=aic;
run;

proc nlmixed data=autoord corr ecorr;

parms b0=0 b1=0 b2=0 sd=1 i1=1 i2=1 i3 =1 sigor = 2;
bounds i1 > 0, i2 > 0, i3>0;
eta = b0 + b1*trt + b2*time + u;
if (value=1) then p = probnorm((-eta)/sigor);
else if (value=2) then
p = probnorm((i1-eta)/sigor) - probnorm((-eta)/sigor);
else if (value=3) then
p = probnorm((i1+i2-eta)/sigor) - probnorm((i1-eta)/sigor);
else if (value=4) then
p = probnorm((i1+i2+i3-eta)/sigor) - probnorm((i1+i2-eta)/sigor);
else p = 1 - probnorm((i1+i2+i3-eta)/sigor);
if (p > 1e-8) then ll = log(p);
else ll = -1e20;
model value ~ general(ll);
random u ~ normal(0,sd*sd) subject=id;

estimate 'thresh2' i1;
estimate 'thresh3' i1 + i2;
estimate 'icc' sd*sd/(1+sd*sd);
ods output ParameterEstimates=pars;
ods output FitStatistics=aic;
run;

I can obtain the results with warnings: WARNING: The final Hessian matrix is full rank but has at least one negative eigenvalue. Second-order optimality condition
violated.

But when I applied the estimated value as starting values for my joint model, I still got the same error.  I wonder if anybody know what should I do next, how to adjust it? (Since I simulate data based on my joint model, I believe using incorrect model is not a problem) Thanks!

[sbxkoenk]

Hello @lcw68 ,

When inserting SAS-code in your post , it's better to click the running man icon in the toolbar (running man in squared box with bullet in lower-right corner) ... then paste your code in the empty box that pops up.

Lay-out, coloring and formatting of your code is optimally preserved that way
 

With regard to starting values for the parameters.
You can start with using the default starting values (No PARMS statement needed!).

When you want to "help" the optimizer, ... this should help when specifying starting values for the parameters for PARMS statement in NLMIXED procedure.

• Reduced models: A way to choose initial parameters for a mixed model
  By Rick Wicklin on The DO Loop June 27, 2018
  https://blogs.sas.com/content/iml/2018/06/27/reduced-choose-parameters-mixed-model.html
• Use a grid search to find initial parameter values for regression models in SAS
  By Rick Wicklin on The DO Loop June 25, 2018
  https://blogs.sas.com/content/iml/2018/06/25/grid-search-for-parameters-sas.html
• The method of moments: A smart way to choose initial parameters for MLE
  By Rick Wicklin on The DO Loop November 27, 2017
  https://blogs.sas.com/content/iml/2017/11/27/method-of-moments-estimates-mle.html
• The DOLIST syntax: Specify a list of numerical values in SAS
  By Rick Wicklin on The DO Loop January 18, 2021
  https://blogs.sas.com/content/iml/2021/01/18/the-dolist-syntax-specify-a-list-of-numerical-values-in...

Good luck,
Koen

[lcw68]

Thank you for your prompt help! These provides lots of different ideas. I was wondering if this warning matters when I run only the ordinal model separately, it shows "WARNING: The final Hessian matrix is full rank but has at least one negative eigenvalue. Second-order optimality condition violated." Seems for ordinal model only, the estimation result could not give the std error for each estimated parameter.

I have tried different kind of starting values, tried to reduce the number of parameters (let threshold for ordinal variable be constant, even std deviation parameter be constant), and neither of them work. 

I also found no matter I changed the initial parameter, the negative LL always keep the same. May this be a problem?

[Rick_SAS]

When you get the warnings that you report, it might be the case that the model is misspecified. Can you say, in words, what type of model you are trying to fit for the ordinal response? Have you double-checked that your implementation of the models (via the LL) are implemented correctly?

In an excellent paper, Robin High (2012) discusses how to fit various kinds of ordinal response models by using PROC LNMIXED. If you haven't seen his paper, I recommend it: Models for Ordinal Response Data

[sbxkoenk]

Hello @lcw68 ,

As @Rick_SAS wrote ... (ordinal) model might be misspecified.

Another common reason for this type of warning messages in the Log for PROC NLMIXED --

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

         Second-order optimality condition violated.

is scaling issues.
For example, when the parameter estimates are on very different scales. Rescaling (i.e. changing order of magnitude) might help. It definitely helps if the parameters are on similar scales!

BR, Koen

[lcw68]

Hi, I use probit normal model to generate ordinal data, it's like y=mu+b+epsilon, y* =0*I(y<0)+ 1*I(0<=y<c1)+2*I(c1<=y<c2)+...+4*I(y1> c3). The estimation result I got from the code below is very close to the true parameter I set, everything seems good except for a warning. After that, when I tried to fit joint model with continuous outcome together, it always shows no valid parameter. 

proc nlmixed data=autoord corr ecorr;
   parms b0=0 b1=0 b2=0 sd=1 i1=1 i2=1 i3 =1 sigor = 2;
   bounds i1 > 0, i2 > 0, i3>0;
   eta = b0 + b1*trt + b2*time + u;
   if (value=1) then p = probnorm((-eta)/sigor);
   else if (value=2) then
      p = probnorm((i1-eta)/sigor) - probnorm((-eta)/sigor);
   else if (value=3) then
      p = probnorm((i1+i2-eta)/sigor) - probnorm((i1-eta)/sigor);
  else if (value=4) then
      p = probnorm((i1+i2+i3-eta)/sigor) - probnorm((i1+i2-eta)/sigor);
   else p = 1 - probnorm((i1+i2+i3-eta)/sigor);
   if (p > 1e-8) then ll = log(p);
   else ll = -1e20;
   model value ~ general(ll);
   random u ~ normal(0,sd*sd) subject=id;
 
   estimate 'thresh2' i1;
   estimate 'thresh3' i1 + i2;
   estimate 'icc' sd*sd/(1+sd*sd);
   ods output ParameterEstimates=pars;
ods output FitStatistics=aic;
run;

I am sorry that I asked many times about it,  but I am very confused and struggled with that. 

[Rick_SAS]

It sounds like you are running this code on simulated data that you believe is generated from this model. Could you share the simulation code you are using?

[jiltao]

If you wanted to specify a UN type of covariance structure for your two random effects, you might use something like the following --

RANDOM z1 z2 ~ NORMAL([0,0],[rhoe1,rhob,rhoe2]) SUBJECT=id;

And sigma2 is already your variance for the normal distribution. I am not sure exactly what rhoe measures. There is no "error" variance for an ordinal outcome, as far as I know. 

Thanks,

Jill