How do I model this in PROC MIXED?

[deleted_user]

I've been told that this is a non-trivial problem, and I should post it here.

Here is a link to the LaTeXed version of the problem description and the model (in matrix form) that I would like to fit:

http://farm4.static.flickr.com/3483/3834750144_e7d9560902_o.jpg

My primary interest is estimating G, or equivalently D, as shown in the link above.

My variables in SAS are:

- effect (the responses, or treatment effects I want to model)
- endp (class variable to represent the two endpoints corresponding to alpha and beta treatment effects)
- center (class variable to represent the trial where the treatment effects were measured)

To get the specific model and covariance structures shown in the link above, I have tried the following (among many others):

proc mixed data=shihco maxiter=500 update info covtest noclprint;
class center endp;
model effect=endp /noint s;
random endp / subject = center type=un g;
repeated endp /group=center type=un r;
parms /parmsdata=parms hold =4 to &nobs;
ods output covparms=covparm;
run;

I am using the parms statement to provide initial values for the elements of G, and fixed values for the elements of R (estimated previously).

When I fit this model under different variations of the above, either there is no convergence, or the Hessian is not positive definite. I've tried a variety of covariance structures for G. The only one that yields no errors is type=UNR under the random statement; however, I do not trust these G estimates.

Should I try changing something else? Am I even using the correct options here in random/repeated to get the model form and covariances that I want? Any help is appreciated. Message was edited by: phdstudent

[SteveDenham]

I have some thoughts, but I don't view any of them as a real solution. They might give some insight into what's happening if you give them a try.

First, what happens to the estimates if you lift the constraint on R to the fixed values? Does the model converge without problems, and if so, how much change do you see in those parameters.

Second, have you tried type=FA0(2) as a covariance structure? It is constrained to be non-negative definite, so maybe something good would happen.

Third, along the same lines, change over to PROC GLIMMIX, and try type=chol, for a Cholesky root parameterization, this would give at least a positive semi-definite matrix.

Fourth, now that you are in GLIMMIX, try at least a couple of different optimization methods.

I don't know that any of these will succeed as I'm not familiar with the data, but since I've long since gone through the effort of pulling out my hair on problems like this, I thought I'd give it a shot as to what I might try.

Good luck.

Steve Denham
Associate Director, Biostatistics
MPI Research, Inc.

[deleted_user]

Thank you for your response.

If I lift the constraint on R to fixed values, we lose convergence, and the initial values are not updated. Same thing if I don't supply initial values at all.

I did try FA0(2) as a covariance structure in PROC MIXED, but the hessian is still not positive definite, and the initial values do not update at all. Same thing if I don't fix R, or if I do not supply initial values (also, no convergence).

I'm now working in PROG GLIMMIX with the "chol" option for type, and things are doing a little better. However, I am getting solutions for G for which correlations calculated from G are greater than 1, etc. Still going to try a few more things here...

Thanks again.