topic Re: How to incorporate correlation between random error (epsilon) in proc nlmixed joint MMRM model? in Statistical Procedures

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

lcw68 — Sat, 29 Jul 2023 03:28:02 GMT

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!

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

sbxkoenk — Sat, 29 Jul 2023 08:59:17 GMT

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

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

lcw68 — Sat, 29 Jul 2023 15:26:19 GMT

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!

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

sbxkoenk — Sun, 30 Jul 2023 10:08:02 GMT

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

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

lcw68 — Sun, 30 Jul 2023 18:56:12 GMT

Thanks very much, I'll read it through

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

lcw68 — Sun, 30 Jul 2023 19:55:16 GMT

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!

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

sbxkoenk — Sun, 30 Jul 2023 20:14:38 GMT

Hello @lcw68 ,

Tomorrow (Monday) is working day again.

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

Koen

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

lcw68 — Sun, 30 Jul 2023 21:22:05 GMT

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!

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

sbxkoenk — Sun, 30 Jul 2023 21:59:18 GMT

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-sas.html

Good luck,
Koen

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

lcw68 — Mon, 31 Jul 2023 06:38:16 GMT

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?

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

Rick_SAS — Mon, 31 Jul 2023 10:03:38 GMT

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

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

sbxkoenk — Mon, 31 Jul 2023 10:25:09 GMT

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

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

lcw68 — Mon, 31 Jul 2023 13:47:43 GMT

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.

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

Rick_SAS — Mon, 31 Jul 2023 14:09:12 GMT

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?

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

jiltao — Mon, 31 Jul 2023 14:27:35 GMT

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

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

lcw68 — Mon, 31 Jul 2023 14:33:00 GMT

Sure, I write the simulation in R since I also need to make calculations to obtain the true parameter. X1 is the dataset I input to SAS.

cind = c(0,c2,c2+2,c2+4) ##threshold varb <- matrix(c(sigma.b0^2,rhob*sigma.b0*sigma.b1,rhob*sigma.b0*sigma.b1,sigma.b1^2),nrow=2) vare <- matrix(c(sigma.e0^2,rhoe*sigma.e0*sigma.e1,rhoe*sigma.e0*sigma.e1,sigma.e1^2),nrow=2) k=2 ### 2 timepoints right now n = 100 ##n samples set.seed(i) trt.pure = rbinom(n,1,2/3) trt = rep(trt.pure,each=k) time = rep(0:(k-1),times = n) id = 1:n X = data.frame(trt = trt, time = time/(k-1), id = rep(id,each=k)) mu1 = model.matrix(~time+trt,X)%*%c(alpha0,alpha1,alpha2) mu2 = model.matrix(~time+trt,X)%*%c(beta0,beta1,beta2) set.seed(i) b <- rmvnorm(n,mean = c(0,0),sigma = varb) eps <- rmvnorm(n*k,mean = c(0,0),sigma = vare) y0 = mu1+ rep(b[,1],each=k) + eps[,1] y1 = mu2+ rep(b[,2],each=k) + eps[,2] y0 = floor(y0) y0[y0>upper.limit] = upper.limit y0[y0<lower.limit] = lower.limit y1.ord <- dplyr::case_when(y1 < cind[1] ~"0", y1>=cind[1] & y1< cind[2] ~ "1", y1>=cind[2] & y1< cind[3] ~ "2", y1>=cind[3] & y1< cind[4] ~ "3", y1>=cind[4] ~"4") X$continuous =y0 X$ordinal = as.numeric(y1.ord)+1 X1 = reshape2::melt(X,id.vars = c("trt","time","id")) save(X1,file="....")

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

lcw68 — Mon, 31 Jul 2023 14:49:07 GMT

Hi, thanks for your reply, what I mean is that (b1, b2) = (z1,z2) ~ N((0,0),(sb1^2, rhob*sb1*sb2,sb2^2)), (eps1, eps2) ~ N((0,0),(se1^2, rhoe*se1*se2,se2^2)).

rhoe measures the correlation between two outcomes in each time points (within one subject), and I found I can't find a place to put rhoe in the code of proc nlmixed joint model.

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

jiltao — Mon, 31 Jul 2023 15:05:27 GMT

PROC NLMIXED cannot model correlated errors. It assumes the errors are independent.

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

lcw68 — Mon, 31 Jul 2023 15:11:14 GMT

Oh I see, that's what I want to confirm. Then in SAS I wonder if you know what kind of procedure could incoporate correlated random effect and correlated error at the same time.

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

jiltao — Mon, 31 Jul 2023 15:27:39 GMT

You can use PROC GLIMMIX with the G-side random effect and R-side random effect to model both correlations. If you have two response variables with different distributions, you can use DIST=BYOBS() in the MODEL statement for that. Here is an example -- https://go.documentation.sas.com/doc/en/statcdc/14.3/statug/statug_glimmix_examples08.htm

Thanks,

Jill