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

[lcw68]

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.

[jiltao · Accepted solution]

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

[lcw68]

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.

[jiltao]

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

[lcw68]

Thanks!

[lcw68]

Hi, I found that in proc glimmix, it only specifys some specific distribution in Byobs, how could I specify "ordinal - probit normal" using byobs?

[jiltao]

do you have two response variables? Can you send me a sample data and I will be happy to provide some syntax.

[lcw68]

Sure, here is my code, trt are group effect and time are time effect. variable indicates what type of data it is.

data sample;
  input trt time id variable $ value;
 datalines;
0 0 1 continuous 21.2080566
0 1 1 continuous 11.800
0 0 1 ordinal 5
0 1 1 ordinal 4
0 0 2 continuous 23.77209112
0 1 2 continuous 6.668272081
0 0 2 ordinal 5
0 1 2 ordinal 4
0 0 3 continuous 27.54257141
0 1 3 continuous 13.36067935
0 0 3 ordinal 5
0 1 3 ordinal 5
0 0 4 continuous 26.86563107
0 1 4 continuous 9.048611888
0 0 4 ordinal 5
0 1 4 ordinal 4
0 0 5 continuous 22.64034662
0 1 5 continuous 7.429331484
0 0 5 ordinal 5
0 1 5 ordinal 3
1 0 6 continuous 24.34331053
1 1 6 continuous 8.31258158
1 0 6 ordinal 5
1 1 6 ordinal 4
0 0 7 continuous 25.70524738
0 1 7 continuous 8.856146754
0 0 7 ordinal 5
0 1 7 ordinal 4
1 0 8 continuous 24.83456925
1 1 8 continuous 10.28819607
1 0 8 ordinal 5
1 1 8 ordinal 4
1 0 9 continuous 24.23330388
1 1 9 continuous 8.515241755
1 0 9 ordinal 5
1 1 9 ordinal 4
;

[jiltao]

Upon closer examination of the documentation, I realized that multinomial is not supported in dist=BYOBS(). So I am afraid that PROC GLIMMIX does not work here, unfortunately.

[lcw68]

It's OK, thank you all the same!

[lcw68]

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="....")

[lcw68]

@Rick_SAS I am totally struggled, I tried to generate two independent variable, a continuous and an ordinal one, then still got the same error. I tried to reduce the number of parameters, make more parameter as constant in the model, but it always gives the same error. BTW, the negative log likelihood at initial parameters didn't change at all, no matter what values, or how many parameters. I wonder What  I could check next...