topic Re: Using the glimmix radial smoother in Statistical Procedures

Using the glimmix radial smoother

Sam28041977 — Tue, 26 May 2015 08:36:37 GMT

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

Re: Using the glimmix radial smoother

SteveDenham — Thu, 28 May 2015 18:14:07 GMT

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

Re: Using the glimmix radial smoother

lvm — Thu, 28 May 2015 21:37:32 GMT

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.

Re: Using the glimmix radial smoother

SteveDenham — Fri, 29 May 2015 13:19:14 GMT

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

Re: Using the glimmix radial smoother

Sam28041977 — Sat, 30 May 2015 18:30:53 GMT

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

Re: Using the glimmix radial smoother

Sam28041977 — Sun, 31 May 2015 03:17:40 GMT

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

Re: Using the glimmix radial smoother

lvm — Sun, 31 May 2015 03:34:41 GMT

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.

Re: Using the glimmix radial smoother

Sam28041977 — Sun, 31 May 2015 03:40:56 GMT

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

Re: Using the glimmix radial smoother

SteveDenham — Mon, 01 Jun 2015 13:15:12 GMT

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.

Steve Denham

Re: Using the glimmix radial smoother

sld — Tue, 02 Jun 2015 19:57:14 GMT

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

Re: Using the glimmix radial smoother

lvm — Tue, 02 Jun 2015 21:03:51 GMT

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).

Re: Using the glimmix radial smoother

sld — Wed, 03 Jun 2015 20:11:49 GMT

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

Re: Using the glimmix radial smoother

SteveDenham — Fri, 05 Jun 2015 17:03:38 GMT

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;

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

Re: Using the glimmix radial smoother

Sam28041977 — Fri, 05 Jun 2015 22:19:54 GMT

Dear SLD,

The first random statement allow me to do within subjects comparison of treatments.

Every subject in the study receive all the treatments.

The second random statement capture the correlation from repeated measurements within the same treatment.

I hope this is clear.

Samir

Re: Using the glimmix radial smoother

Sam28041977 — Fri, 05 Jun 2015 22:26:21 GMT

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;

run;

I hope this will work better.

Samir

Re: Using the glimmix radial smoother

sld — Fri, 05 Jun 2015 22:28:18 GMT

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

Re: Using the glimmix radial smoother

sld — Sat, 06 Jun 2015 17:36:26 GMT

Hello Samir,

Yes, that clears up my concern. Thank you. SUJET(TRT) is read as "sujet nested within trt" so my guess was that SUJET was the experimental unit assigned to a single level of TRT, rather than SUJET being a block containing smaller experimental units to which levels of TRT were assigned (in your words, each subject receives all treatments). Because both SUJET and TRT were defined in your model, I think that GLIMMIX would interpret SUJET(TRT) as SUJET*TRT which makes sense.

Susan

Re: Using the glimmix radial smoother

SteveDenham — Mon, 08 Jun 2015 13:23:43 GMT

Hi Susan,

It seems like I am hijacking Samir's thread, but it still all boils down to using splines in a mixed model setting.

As far as selecting which to go with--the missing data means the two approaches give very different results. The group that is early right censored gives a continually increasing estimate under the "only random effect spline" which is probably closer to the biology. The fixed-effect spline/random spline model yields values that decrease once the data starts to be missing en masse--which is not unexpected, as the least squares estimates will regress toward the mean under the MAR assumption that usually applies to REML.

Now I am waiting for some NLMIXED code that incorporates a general log likelihood function for a 3 parameter Gompertz distribution/growth model.

In the end, I'll probably end up selecting the model that minimizes the standard error of predicted values across the range of non-missing data.

Steve