Programming the statistical procedures from SAS

proc mixed estimated V matrix for different subjects

Reply
Occasional Contributor RVS
Occasional Contributor
Posts: 5

proc mixed estimated V matrix for different subjects

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.

Respected Advisor
Posts: 2,655

Re: proc mixed estimated V matrix for different subjects

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

Occasional Contributor RVS
Occasional Contributor
Posts: 5

Re: proc mixed estimated V matrix for different subjects

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

Occasional Contributor RVS
Occasional Contributor
Posts: 5

Re: proc mixed estimated V matrix for different subjects

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

Respected Advisor
Posts: 2,655

Re: proc mixed estimated V matrix for different subjects

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

Ask a Question
Discussion stats
  • 4 replies
  • 419 views
  • 0 likes
  • 2 in conversation