PROC MIXED SLICE STATEMENT

[beahern]

Hi,

I am comparing 3 treatments on a continuous measure administered at 3 time points (1, 2, 3) via a mixed model using the following code:

proc mixed data=dataset METHOD=ML COVTEST EMPIRICAL;
  class RECORD_ID time(ref='1') treatment(ref='CT');
  model ISI_Total = TREATMENT TIME age female TREATMENT*TIME / SOLUTION DDFM=BW;
  RANDOM intercept / TYPE=UN SUBJECT = record_id(TREATMENT);
  SLICE TREATMENT*TIME / SLICEBY TIME;
  LSMEANS TIME TREATMENT*TIME / DIFF ;
  TITLE1 "Mixed Model: ISI_Total predicted by fixed Treat, Time and Treatment*Time, random intercept";
  TITLE2 "CT is the reference treatment";
run;

The SLICE statement generates an F-test for TREATMENT within each TIME that I need, but I also need to generate an F-test for change between time points (change from TIME=1 to TIME=2, TIME=1 to TIME=3, TIME=2 to TIME=3) by treatment within the same mixed model. I am having a hard time modifying my code to generate these F-tests within the SLICE statement.  I would appreciate any guidance. 

[ChrisNZ]

Moved to stats community.

[SteveDenham]

I think this will be a case where GLIMMIX offers what you need but MIXED doesn't without writing LSMESTIMATE statements.  Consider trying this:

proc glimmix data=dataset METHOD=MSPL  EMPIRICAL;
  class RECORD_ID time(ref='1') treatment(ref='CT');
  model ISI_Total = TREATMENT TIME age female TREATMENT*TIME / SOLUTION DDFM=BW;
  RANDOM intercept / TYPE=UN SUBJECT = record_id(TREATMENT);
    LSMEANS TIME  / DIFF   ;
    LSMEANS TREATMENT*TIME/SLICEDIFF = TREATMENT SLICEDIFFTYPE=ALL;
COVTEST GLM/ESTIMATES CLASSICAL CL(TYPE=PLR);
  TITLE1 "Mixed Model: ISI_Total predicted by fixed Treat, Time and Treatment*Time, random intercept";
  TITLE2 "CT is the reference treatment";
run;

I am a bit curious about the selection of type=UN for the random intercepts.  If you have 10 subjects, that means estimating 45 parameters, and this grows like a weed with more subjects (100 subjects-->4950 parameters).  Unless you have a lot of data, you are going to get a 'G matrix is not positive definite' message in the output. Instead, try type=UN(1), which would give one variance component for each subject.  If that still leads to a problem, try dropping type= altogether, which will give one estimate for the variance of the intercepts. 

SteveDenham

[StatsMan]

@SteveDenham is correct that the SLICEDIFF= option should give you what you are looking for here.

One additional point. With a RANDOM statement like

random intercept / type=un subject=id(trt);

the TYPE= speicifcation does nothing for you here. You get two variance components in the output from GLIMMIX, one for the ID(TRT) effect and one for the residual. You can drop the TYPE=UN (or change the covariance specification to just about anything else) and the model results will not change. This behavior only works if we are talking about g-side covariance structures. 

If you change the RANDOM statement to

random id(trt) / type=un;

then you get the behavior that is described in the previous reply, a large number of variance/covariance parameters and most likely a convergence issue with the model. 

Why the differenece? The first RANDOM statement allows you to process the V (ZGZ' + R) matrix in blocks, one block per subject, and apply the common covariance across all observations from the same subject. The second RANDOM statement does not block the V matrix and applies a unique variance to each ID(TRT) level and a unique covariance across each pair of ID(TRT) levels. This can be a helpful structure in some agricultural studies, where you want to correlate observations from adjacent plots, but as Steve points out you really need a small number of levels to ID(TRT) to make this problem approachable. 

[beahern]

Thank you very much for taking the time to explain. This is very helpful.

[SteveDenham]

THANKS!!! @StatsMan , I have always been unsure about the difference between the two, and made a really wrong assumption that the two would result in the same number of parameters.  Now, I won't run away from UN and CHOL as fast as in the past.

SteveDenham

[beahern]

Hi @SteveDenham,

Thank you for your help with this problem.  We are trying to replicate results from a published study (tables 2, 3, 4 in enclosed article) using a different sample, and unfortunately, not enough details are provided about the SAS code. 

The code you suggested produces pairwise comparisons rather than an overall F-test. Within PROC MIXED, the SLICE statement produces the F-test I need for treatment on each value of TIME (time=1 and 2 included below). I also need a similar test for differences in times (Change in TIME1 - TIME 2, for example).  Apparently all these comparisons were run within the same mixed model.

[SteveDenham]

You can get those F tests using a SLICE statement in GLIMMIX.  It just isn't the best way to extract the pairwise comparisons.  So we would now have this (incorporating @StatsMan 's comments as well).

proc glimmix data=dataset METHOD=MSPL  EMPIRICAL;
  class RECORD_ID time(ref='1') treatment(ref='CT');
  model ISI_Total = TREATMENT TIME age female TREATMENT*TIME / SOLUTION DDFM=BW;
  RANDOM intercept /  SUBJECT = record_id(TREATMENT);
    LSMEANS TIME  / DIFF   ;
    SLICE TREATMENT*TIME/SLICEBY = TREATMENT DIFF=ALL;
COVTEST GLM/ESTIMATES CLASSICAL CL(TYPE=PLR);
  TITLE1 "Mixed Model: ISI_Total predicted by fixed Treat, Time and Treatment*Time, random intercept";
  TITLE2 "CT is the reference treatment";
run;

I believe you should consider modeling the covariance structure for TIME with something like:

RANDOM TIME/RESIDUAL SUBJECT=record_id(TREATMENT) type= <pick something appropriate to the number of observations and their spacing in time >;

as the measurements at each time point are very likely not independent of one another.

SteveDenham

.

[StatsMan]

I think the differences you see in the associated pdf can be obtained from differences in the lsmeans, something you can get with the PDIFF option on the LSMEANS statement. If you want differences of differences of lsmeans, then as Steve mentioned earlier, you will need to get familiar with the LSMESTIMATE statement (or the ESTIMATE statement if you prefer).