Issues with weight statement in a multifactorial repeated measures model in glimmix

[palolix]

Dear SAS Community,

I'm having issues when trying to run a multifactorial RCBD model that includes depth as a split block factor, date as repeated measures, and a weight statement since I'm combining data of two experiments. I get WARNING: Obtaining minimum variance quadratic unbiased estimates as starting values for the covariance parameters failed.  Although I noticed that if I remove the weight statement, I stop getting that warning and SAS is able to provide results. Also, for another dependent variable I don't have any problem when using the same model.  

I went back to the model used to get the residuals for that variable and even though I got results, I noticed this warning:  MIVQUE0 estimate of profiled variance is linearly related to other covariance parameters.

This is the model I used to obtain residuals to calculate weights:

Proc glimmix data=one;
by Exp;
class blk Trt Rate Depth Date;
model RL=Trt|Rate|Depth|Date/dist=lognormal ddfm=bw/*to get residual variances*/;
random intercept Trt*Rate Depth Trt*Rate*Depth Date Trt*Date Depth*Date Trt*Rate*Depth*Date/subject=blk;
random Date/residual subject=blk(Trt*Rate Depth Trt*Rate*Depth Date Trt*Date Depth*Date Trt*Rate*Depth*Date) type=cs;
lsmeans Trt*Rate*Depth*Date;/*to get the error DF*/
run;

I wonder if the random statements are correct (including date as repeated measures, and depth as split-block factor within the model).

Thank you so much for your support!

Caroline

[jiltao]

Your R-side random effect specification seems to be redundant to the G-side random effects. I would remove the second RANDOM statement.

Hope this helps,

Jill

[jiltao]

I forgot to point out that your last effect in the RANDOM statement is not necessary -- it confounds with the residual variance. So try this --

Proc glimmix data=one;
by Exp;
class blk Trt Rate Depth Date;
model RL=Trt|Rate|Depth|Date/dist=lognormal ddfm=bw/*to get residual variances*/;
random intercept Trt*Rate Depth Trt*Rate*Depth Date Trt*Date Depth*Date /subject=blk;
lsmeans Trt*Rate*Depth*Date;/*to get the error DF*/
run;

Hope this helps,

Jill

[palolix]

Thank you Jiltao!

Ok, I changed that in the model to get new residuals (no warning this time). Now, for the final weighted analysis, should I also eliminate the last interaction (Trt*Rate*Depth*Date) in the random statement, but keeping the repeated measures random  statement? 

Then it should look like this:

Proc glimmix data=one;
class blk Trt Rate Depth Date;
model RL=Trt|Rate|Depth|Date/dist=lognormal ddfm=kr;
random intercept Trt*Rate Depth Trt*Rate*Depth Date Trt*Date Depth*Date /subject=blk;
random Date/residual subject=blk(Trt*Rate Depth Trt*Rate*Depth Date Trt*Date Depth*Date) type=cs;
weight RL_wt;
run;

Thank you 

Caroline

[jiltao]

Your subject= effect in your second RANDOM statement is too convoluted. It needs to be simplified.

You also need the BY statement. So try the following program --

Proc glimmix data=one;
by exp;
class blk Trt Rate Depth Date;
model RL=Trt|Rate|Depth|Date/dist=lognormal ddfm=kr;
random intercept Trt*Rate Depth Date Trt*Date Depth*Date /subject=blk;
random Date/residual subject=blk(Trt*Rate*Depth) type=cs;
weight RL_wt;
run;

The program runs. However, you will see that the G matrix is npd. You might need to remove some of the random effects in your first RANDOM statement.

[palolix]

Thank you Jiltao,

It seems like I will not be able to combine both experiments and will have to analyze them separately. 

I completely agree that the random statements are too convoluted, but I'm afraid that if I start removing interactions it will no longer fit the original model with the strip factors and so on. When I work with proc glm, I pool error terms if possible, to simplify models, but I don't know how to do something analogous in glimmix. 

Thank you!

Caroline

[jiltao]

If you want to combine EXP in your analysis, rather than doing the analysis separately, you would need to include EXP in your model, such as in the CLASS, MODEL and RANDOM statements.

The strip factors and so on are modeled in your first RANDOM statement (the G-side random effects). The second RANDOM statement (the R-side random effect) only models the correlations in the residuals across DATE for a subject, where the subject is identified by blk(Trt*Rate*Depth) [or blk(Exp*Trt*Rate*Depth) if you want to include EXP in your analysis]. I am actually surprised that SAS did not complain your unusual subject specification there.

[palolix]

Thank you againJiltao!

Thats very good to know about the G- and the R-side random effects (I changed that already in the model).  I don't want Exp as another factor included in the model, that's why I'm trying to do a weighted analysis, so I can combine the data of the two experiments (which are replications).

If I include date in the R-side random statement, should I then remove it from the G-side random statement?

Thank you!

Caroline

[jiltao]

Date in the R-side random statement just identifies repeated effect. In your case (with TYPE=CS) the results should be the same with or without the DATE in the R-side random effect. So no, no need to remove DATE from the first (G-side) RANDOM statement. But of course, if you wanted not to model the DATE random effect for each BLK, you can certainly remove that, but that has nothing to do with the syntax in the R-side random effect specification.

[palolix]

Thank you Jiltao,

You were all right, with 4-way interactions the model was too overwhelmed and that was probably the reason of the warnings (also SAS was taking way too long to deliver results). If I instead analyze it by date, the model is simpler, and I get clean results without any warning. 

The last question: Is the random statement correct for a RCBD that includes Depth as a strip plot factor?

Proc glimmix data=one;
class blk Trt Rate Depth;
model RL=Trt|Rate|Depth/dist=lognormal ddfm=kr ;
random intercept Trt*Rate Depth Trt*Rate*Depth/subject=blk;
weight RL_wt;

run;

Thanks a lot!

Caroline

[palolix]

Thank you Jiltao for your quick reply!

Isn't the second random statement necessary to indicate SAS to treat date as repeated measures? 

Even if I remove it, I still get the same warning and no results.

Thank you,

Caroline

[StatsMan]

Once you understand what the RANDOM statements do in GLIMMIX, you will have a better understanding of whether or not you have fit the model you want. With a RANDOM statement like

     random int / subject=blk;

you are correlating the errors (the elements of the V matrix in your mixed model) for observations one the same level of BLK. This gives you a blocked covariance structure in V, where you have that common covariance among all the rows and columns corresponding to the same level of BLK. Add TRT to the RANDOM statement

     random int trt / subject=blk;

Now you have a second level of covariance. You still have the same covariance for all observations from the same level of BLK. On top of that, you have a common covariance for all observations with the same levels of TRT and BLK. As you add more terms to the RANDOM statement, you add more block structure to the V matrix. This discussion assumes that all terms (like TRT) are on the CLASS statement.  

Add the V option to the RANDOM statement (this works in PROC MIXED, too), and you can see the block structure for the first subject in your data (the first level of BLK). V= allows you to print the blocks of V for other subjects, using a numerical ordering of the subjects in your data. V=3 prints the block of V for the 3rd subject. Try the V option with different RANDOM statements to see how the results change.

As for the RANDOM RESIDUAL statement, if you have 

     random int / subject=blk;

and

     random date / residual subject=blk type=cs;

then those two statements fit the same covariance structure to the observations from the same level of BLK. You cannot do that. You get the message about the confounding message you have seen in your output in this case. The second RANDOM statement here is redundant and you can and should remove it.

[palolix]

Hi StatsMan,

Thank you very much for the explanation. I will definitely try the V option. I'm still a bit confused because I thought that the second random statement with the residual option was the Glimmix way of dealing with repeated measures.

Thank you

Caroline

[StatsMan]

RANDOM/RESIDUAL is the way to fit R-side covariance structures in GLIMMIX. Those are most commonly used in repeated measures models. However ...

These two RANDOM statements 

     random int / subject=blk;

and

     random date / residual subject=blk type=cs;

fit the same covariance structure in V when no link function is present. You can use V to see that. The first statement does this through G-side random effects, while the second sets this up through the R-side. This only works when there is no link function (ie, you are using the identity link) and when you have TYPE=CS on the second RANDOM statement. Add a link function or use a different TYPE= (like AR(1)) and the two statements are not equivalent. This holds true in PROC MIXED as well, with a RANDOM statement and a REPEATED statement using TYPE=CS.

Does that help?

[palolix]

Thank you very much for the clarification StatsMan!