topic Re: proc mixed estimated V matrix for different subjects in Statistical Procedures

proc mixed estimated V matrix for different subjects

RVS — Wed, 24 Jun 2015 20:34:55 GMT

Hello SAS community,

Following is the model that I'm interested to extract blocks of the V matrix using hold option with known or predefined variance parameters. I would like to get V matrix that contains G-side and R-side structure at plot level, i.e.a matrix of 240 x 240 in the following example.

I was able to get V matrix for location as subject, but I'm not sure how to change the subject to a plot level to get V with both G and R structures.
*Generate numbers;
data one;
do har_y = 1 to 2; * two harvest years;
   do loc = 1 to 5; * five locations ;
    do blk = 1 to 4; * four blocks;
     do cult = 1 to 6;*six cultivars;
       yield= 1; * dummy value
         output;
      end;
     end;
    end;
   end;
*Sort
proc print; run;
proc sort data=one out=two;by har_y loc blk cult;
run;

proc mixed lognote data=two;
class cult loc har_y blk;
model yield=cult har_y cult*har_y;
random har_y/subject=loc type=csh;
random har_y/subject=loc*cult type=csh;
random har_y/subject=loc*blk type=csh v ;
repeated har_y/subject=loc*blk*cult type=csh r;
parms(6.1387) (3.9508) ( 0.7087) (0.09453) (0.1454) (0.8500) (0.01970) (0.04729) (0) (0.8706) (0.7041) (0.5151)/noiter hold= 1 2 3 4 5 6 7 8 9 10 11 12; *parameters used to get V matrix
ods output Tests3 = vcfixed CovParms=vcparms FitStatistics= vcfitstat v=v1 r=r1 ;
run;

Thank you in advance for your time and consideration.

Re: proc mixed estimated V matrix for different subjects

SteveDenham — Fri, 26 Jun 2015 17:11:42 GMT

This isn't going to help at all, but I am surprised that you don't get a message in the output that the G matrix is not positive definite.

It appears that the repeated statement has a subject of loc*blk*cult, which I assume is at the plot level. With what you have, I guess I would have approached this somewhat differently:

proc mixed lognote data=two;

class cult loc har_y blk;

model yield=cult har_y cult*har_y;

random intercept cult blk/subject=loc type=csh v;

repeated har_y/subject=loc*blk*cult type=csh r;

Note that this fits 3 random effects (loc, loc*cult and loc*blk) each with different values and constant covariances between all effects. Putting harvest year in as both a random and repeated effect seems hard to justify in my opinion--I would see it as random if measurements were on different plots in the two years, and repeated if on the same plot.

See if this helps any.

Steve Denham

Re: proc mixed estimated V matrix for different subjects

RVS — Fri, 26 Jun 2015 19:38:38 GMT

Thanks, Steve for your suggestion on the model fitting. I was trying different variance-covariance structures on the field trials data set that consists of multiple harvest.Yes,I had a problem with not positive definite G matrix (only for few trials), but I resolved that by providing different initial parms. The three random statements in my original question was to model har_y at three random subject levels, i.e. loc, loc*cult and loc*blk.

RVS

Re: proc mixed estimated V matrix for different subjects

RVS — Fri, 26 Jun 2015 20:39:25 GMT

Hi Steve,

I did replicate proc mixed code into proc glimmix, this provided a V matrix consists of 240 x 240 to the plot level. Do you think it is appropriate?

proc glimmix data=two;

class cult loc har_y blk;

model yield=cult har_y cult*har_y;

random har_y/subject=loc type=un;

random har_y/subject=loc*cult type=un;

random har_y/subject=loc*blk type=un;

random har_y/ residual subject=loc*blk*cult type=un v;

ods output Tests3 = vcfixed CovParms=vcparms FitStatistics= vcfitstat v=glimmix ;

run;

Thank you for your time and consideration.

RVS

Re: proc mixed estimated V matrix for different subjects

SteveDenham — Mon, 29 Jun 2015 13:04:03 GMT

If you are looking to fit a random slopes model (in har_y), this looks good. I still think of it as a random intercepts model, with a single G side random statement:

random intercept cult blk/subject=loc type=un;

I end up thinking this because of the R side factors including all of these variance components in an unstructured covariance structure matrix. My experience would be that the output would throw a message that the G matrix was not positive definite--basically saying that the R part was explaining all of the variability. However, with sufficient data, it could be done.

Steve Denham