Programming the statistical procedures from SAS

Using the glimmix radial smoother

Reply
Occasional Contributor
Posts: 12

Using the glimmix radial smoother

Dear Everybody,

I hope you can provide some help about my issue.

My issue is about some particular mean comparison using the glimmix procedure (using the radial smoother).

My subject is about a normally distributed continuous outcome.

I have more than 10 post baseline measurements.

 

I planned first to use a random coefficient model to estimate an intercept and a slope by subject using the mixed procedure to take advantage from the shrinkage property of the model but the kinetic is not linear.

I discovered the radial smoother of the proc Glimmix and the model fitting is very good.

The normality of the residuals is also very good.

With this model, I can easily compare two treatment means (treatment a vs treatment b) at a given time-point: (Ta-Tb) at day 11

But what I want is the compare two mean differences: (Ta-Tb) at day 39 vs (Ta-Tb) at day 11

This comparison may be estimated easily with the mixed procedure with the estimate statement.

I prefer the radial smoother because I can detect smaller differences (high power: standard errors of 0.15 vs 0.60) and and the model properties are very good.

The other issue with the mixed procedure is that I cannot take benefit from the shrinkage property of the mixed model in this context.

Do you think that I can do this particular comparison with the radial smoother with time as a continuous covariate?

I use this syntax to estimate the mean difference at a given time point:

PROC GLIMMIX DATA=A PLOTS=ALL;

    CLASS USUBJID FORMULE;

MODEL Y= FORMULE LOGTIME FORMULE*LOGTIME;

RANDOM LOGTIME /SUBJECT=USUBJID TYPE=RSMOOTH KNOTMETHOD=EQUAL(10);

    LSMEANS FORMULE/AT LOGTIME=1.945910149 DIFF;

RUN;

Thanks in advance for your feedback.

Samir

Respected Advisor
Posts: 2,655

Re: Using the glimmix radial smoother

I do not know what version of SAS/STAT you are using, but if it is late 9.3 or after, you should look at the use of the LSMESTIMATE statement.  It is specifically designed to create tests of the sort you want to look at.

Suppose you have two treatments and two days (example only).  Then something like

LSMESTIMATE treatment*time 'Difference in differences' 1 -1 -1 1;

Will give the difference between the treatment differences over time.

Steve Denham

Valued Guide
Valued Guide
Posts: 684

Re: Using the glimmix radial smoother

LSMESTIMATE is just for means, not continuous variables (i.e., with class variables). So, unfortunately, you cannot use it with the example (time is continuous).

The way I understand the question, you want a contrast of fixed effects. It doesn't matter if you have a radial smooth random effect, a simple variance component, or other random-effect structure. You can certainly get the desired contrast with an ESTIMATE statement. It gets tricky.  Here is a different example I put together quickly with treat with two levels (factor) and time as continuous. I show how to get the expected value for each treatment at time 1 and at time 3. Then the difference of treatments at each time (just subtract one from the other, some things cancel). Then the difference of the difference to get the difference of the two treatments at time 1 and at time 3 (terms cancel again).  All of this can also be done using the newer nonpositional syntax in the estimate statement.

proc glimmix data=b;

class indiv treat;

model y = treat time treat*time / s ;

random int / sub=indiv;

estimate 'trt1 time 1' int 1 treat 1 0 time 1 treat*time 1 0;

estimate 'trt2 time 1' int 1 treat 0 1 time 1 treat*time 0 1;

estimate 'trt1 time 3' int 1 treat 1 0 time 3 treat*time 3 0;

estimate 'trt2 time 3' int 1 treat 0 1 time 3 treat*time 0 3;

estimate 'diff @ time 1' treat 1 -1 treat*time 1 -1; *just do the subtractions;

estimate 'diff @ time 3' treat 1 -1 treat*time 3 -3; *time main effect cancels out;

estimate 'diff of differneces' treat*time -2 2; *just subtract the above two. some things cancel and coefficients change;

run;

Hope this gives you guidance. If you have more than two levels to your treatment factor, things will change. Always look at the solution table (option /s in model statement) to see the order the procedure is using for the parameters.

Respected Advisor
Posts: 2,655

Re: Using the glimmix radial smoother

Good catch on that, Larry.

Here is another approach--that combines a categorical time variable with the radial smoother, and enables the use of the LSMESTIMATE statement.

PROC GLIMMIX DATA=A PLOTS=ALL;

    CLASS USUBJID FORMULE TIME;

LOGTIME=LOG(TIME);

MODEL Y= FORMULE TIME FORMULE*TIME;

RANDOM LOGTIME /SUBJECT=USUBJID TYPE=RSMOOTH KNOTMETHOD=EQUAL(10);

   LSMESTIMATE FORMULE*TIME "Differences of differnces,assuming only two FORMULE and two time points" 1 -1 -1 1;

RUN;

Note that under all of these approaches, only the covariance structure has the radial smoother applied. Yet another approach might be to apply a spline basis to the fixed effects. I am curious as to behavior and performance on real data.

Steve Denham

Occasional Contributor
Posts: 12

Re: Using the glimmix radial smoother

Great Thanks Larry and Steve,

The solution of steve is correct.

I tried it and the difference of differences is correct.

Neverthless, this result disapointed me.

the "estimate 'diff of differneces' treat*time -2 2;" don't depend of the time location but only on the time interval lenght.


this estimated difference is always the  same for a fixed time interval !!


the difference of differences between day 1 and day 2 is equal to one obtained between day 2 and day 3.

This result is only valid for linear kinetics.

I ask myself now if the estimate of difference at a given time estimated by the lsmeans (using the "AT" option) or by the estimate statement are correct.

The solution of steve is also correct but I have a high number of time levels (189), the associated standard errors are very high compared to those obtained with the continuous time covariate (assuming a comparable effcet size).

Thanks a lot,

Samir

Occasional Contributor
Posts: 12

Re: Using the glimmix radial smoother

Hi,

I finally understood why my difference doen't change whith shifting the time interval..

My mean model assume a linear relationship between time and response even the splin kinetics of each subject.

I changed the structure of the mean model by using the spline function in the effect statement of the procedure.

I get now an overall kinetic close to the individual kinetic.

The problem of the estimation of the the difference of difference is still not resolved.

My new model is now:

PROC GLIMMIX DATA=DERIVED PLOTS=STUDENTPANEL;

    CLASS SUJET TRT JOUR ACTIVITY;

  EFFECT spl=spline(TIME/KNOTMETHOD=PERCENTILES(100));

    MODEL RESP= TRT TRT*spl/ S DDFM=BW;

    RANDOM INT/SUBJECT=SUJET;

    RANDOM TIME /SUBJECT=SUJET(TRT) TYPE=RSMOOTH KNOTMETHOD=KDTREE(BUCKET=100 KNOTINFO);

  OUTPUT OUT=PREDJ9 PRED=PRED;

    ESTIMATE '0' TRT 1 -1 TRT*SPL [1,1 500] [-1,2 500];

RUN;

Thanks in advance for your feedback.

Samir

Valued Guide
Valued Guide
Posts: 684

Re: Using the glimmix radial smoother

Glad you figured out why you get the result with the linear model for time.  The same process I showed with the linear model works with a spline. Define the estimate statements for each treatment at each relevant time. Then get the difference of the two treatments at each of two single times. Then get this difference of these two differences, which is what you want. You can skip the other estimates once you have the one you really want. It is worth starting with all the different estimate statements so that you can check your work by looking at your data.

Occasional Contributor
Posts: 12

Re: Using the glimmix radial smoother

Hi,

I finally resolved my problem with this syntax (with the model above):

    ESTIMATE 'ESTIMATE AT TIME 500' TRT 1 -1 TRT*SPL [1,1 500] [-1,2 500];

    ESTIMATE 'ESTIMATE AT TIME 1500' TRT 1 -1 TRT*SPL [1,1 1500] [-1,2 1500];

    ESTIMATE 'ESTIMATE OF THE DIFFERENCE OF DIFFERENCES' TRT*SPL [1,1 500] [-1,2 500] [-1,1 1500] [1,2 1500];

Many thanks for your help.

Samir

Respected Advisor
Posts: 2,655

Re: Using the glimmix radial smoother

Great work!

I would like the next version of SAS/STAT to be able to implement something that readily incorporates what is going on in the EFFECT statement into the equivalent of LSMESTIMATE, mostly because I'm too lazy to do a good job writing ESTIMATE statements and thus get frustrated easily. Smiley Wink

Steve Denham

PROC Star
PROC Star
Posts: 188

Re: Using the glimmix radial smoother

It seems to me that

    RANDOM INT/SUBJECT=SUJET;

is already specified by

    RANDOM TIME /SUBJECT=SUJET(TRT) TYPE=RSMOOTH KNOTMETHOD=KDTREE(BUCKET=100 KNOTINFO);

assuming that SUBJECT=SUJET produces the same model structure as SUBJECT=SUJET(TRT)--in other words, each value of SUJET is unique. And there's a double-specification of splines, one in the MODEL statement as a fixed effect and one in RANDOM as a random effect that could be a concern. I don't know what to think about that.

Example 40.6 in the GLIMMIX procedure documentation suggests another way to allow the smoothing parameters to vary with TRT levels. In your case, code might be something like

PROC GLIMMIX DATA=DERIVED PLOTS=STUDENTPANEL;

    CLASS SUJET TRT JOUR ACTIVITY;

    MODEL RESP= TRT TIME / S DDFM=BW;

    RANDOM INT/SUBJECT=SUJET(TRT);

    RANDOM TIME /SUBJECT=TRT TYPE=RSMOOTH KNOTMETHOD=KDTREE(BUCKET=100 KNOTINFO);

RUN;

or the two RANDOM statements could be replaced by

    RANDOM TIME /SUBJECT=SUBJET(TRT) GROUP=TRT TYPE=RSMOOTH KNOTMETHOD=KDTREE(BUCKET=100 KNOTINFO);

Best,

Susan

Valued Guide
Valued Guide
Posts: 684

Re: Using the glimmix radial smoother

If SUJET is coded as you indicate, your first comment could be right. But I can't tell how the OP codes this variable. SUJET may not be unique across all observations. But I think (?) that the first random statement may not be redundant with the radial smoother; the former is giving a constant covariance within SUJET, and the radial smoother is adding a time-dependent modification to that covariance. A bit like the AR(1)+RE structure on page 181 of Littell et al. (2006; SAS for Mixed Models, 2nd edition). This comes from combining CS and and AR(1) (see also page 177 in this book). I would have to play around with this, probably with an example, to really figure this out.

As far as double splines, I would normally agree with you. However, Walt Stroup does the same thing in his 2013 book (see section 15.3, and 15.4.5.1). A bit like a random coefficient model. So, I see justification for this approach. However, the spline for the fixed effect is calculated differently from the radial spline in the random statement. The knots are at different places and the basis function is different. If I were doing it, I would try to line these up directly. But this would require some specialized coding.

Example 40.6 is a different, and definitely a good, way to approach the problem. I have used this approach for other purposes. Basically, one then uses BLUPs to get Y_hat for the interaction terms. This would make the contrasts of interest trickier (but achievable with some work).

PROC Star
PROC Star
Posts: 188

Re: Using the glimmix radial smoother

Chapter 15 in Walt's book is the one chapter I haven't read, so thank you for pointing out his examples. Walt is very clever with modeling. I really ought to work through his book again--or at least all the sections I've tabbed for a second look.

At the Kansas meeting in April he said they were working on a new edition of SAS for Mixed Models. Something to look forward to!

For those following this discussion, this paper

Littell RC et al, 2000, Modelling covariance structure in the analysis of repeated measures data, Statistics in Medicine 19:1793-1819

has a few sentences (p 1801) about what the AR(1)+RE structure is doing, which is illustrated with Fig 5.4 in Littell et al (2006, p 183).

Susan

Respected Advisor
Posts: 2,655

Re: Using the glimmix radial smoother

How timely is all of this?.  Right now, on my other machine, I am using the following to model tumor volume growth rates.  The data are "nearly evenly" spaced over 63 days. What happens that makes this tricky is that animals are removed from study for humane reasons when the exoplanted tumor reaches a certain size--so in this case, the control animals are censored at around day 40, and other groups have sporadic censoring.  Certainly not a missing at random situation.

/* Fixed effect spline approach */

proc glimmix data=dat.longexp1 chol;

t2=time

val2=value/100;

class grp_no time anml_no;

nloptions maxiter=5000 tech nrridg;

effect spl=spline(time);

model val2=grp_no grp_no*spl/noint ddfm=kr2;

random time/residual type=chol subject=anml_no;

lsmeans grp_no/at time=1;

<repeat this part for all 19 time points>

run;

I realize I end up estimating 133 variance-covariance parameters here, but I have 75 subjects in 5 groups, so this should be tractable.

THE KEY THING I HAVE LEARNED FROM THIS PART:

When fitting fixed effect splines, GLIMMIX has a hard time getting started if you specify a spline main effect and a spline by treatment interaction.  All that is needed is the spline by treatment term.  In this case, it is like fitting separate intercept and 'slope' fixed effects.

I have also run this with a fixed effect of time (in a factorial design) with the radial smoother.

/* Random smoother approach */

 

proc glimmix data=dat.longexp1;

t2=time;

class grp_no time anml_no;

effect spl=spline(time);

model value=grp_no|time/ ddfm=kr2;

random t2/ type=rsmooth subject=anml_no group=grp_no;

covtest homogeneity;

lsmeans grp_no*time;

run;

and , am I missing anything obvious here?

Steve Denham

Occasional Contributor
Posts: 12

Re: Using the glimmix radial smoother

Dear Steve,

In your first model, the high number of covariance parameters are due to the fact that the time is entered as class variable.

I suggest you this model with time as a continuous covariate:

proc glimmix data=dat.longexp1;

class grp_no anml_no;

effect spl=spline(time);

model val2=grp_no grp_no*spl/noint ddfm=kr2;

random time/subject=anml_no TYPE=RSMOOTH KNOTMETHOD=KDTREE(BUCKET=100 KNOTINFO) ;

lsmeans grp_no/at time=1;

<repeat this part for all 19 time points>

run;

I hope this will work better.

Samir

PROC Star
PROC Star
Posts: 188

Re: Using the glimmix radial smoother

Hi Steve,

I tried your two models on the cow data from the SAS example (40.6 in GLIMMIX) and, apart from me not being patient enough to wait out the CHOL covariance structure, they worked just fine. I can't think of anything amiss. Two distinct approaches to fitting a visually similar model, what are your thoughts about how to choose between them?

Susan

Ask a Question
Discussion stats
  • 17 replies
  • 786 views
  • 4 likes
  • 4 in conversation