proc nlmixed - test if the interaction term is significant at a specific value of another variable

[kc]

I have the following statement in proc mixed:

model var = var0 trt time time2 time*trt time2*trt Category time*Category time2*Category trt*Category time*trt*Category

I have the following statement in proc nlmixed based on the model statement from proc mixed above:

mu1 = beta0 + beta1*var0 + beta2*(trt=1) + beta3*time + beta4*time2 +

      beta5*time*(trt=1) + beta6*time2*(trt=1) + beta7*(Category=2) + beta8*(Category=1) + beta9*time*(Category=2) +           

      beta10*time*(Category=1) + beta11*time2*(Category=2) + beta12*time2*(Category=1) + beta13*(trt=1)*(Category=2) +

      beta14*(trt=1)*(Category=1) + beta15*time*(trt=1)*(Category=2) + beta16*time*(trt=1)*(Category=1) + a ;

There are 11 terms in the original linear mixed model.

var – dependent variable

var0 – continuous variable

trt – treatment groups (0, 1) – reference group is 0

Category – 3 categories (0, 1, 2) – reference group is 0

time – linear time

time2 – quadratic time

beta0 to beta16 – coefficients from linear mixed model

What I need is the overall p-value for the interaction 'trt*Category' at time=2 from the proc nlmixed procedure.

Any help is greatly appreciated.

Thanks!

[StatDave · Accepted solution]

You would do that with a two row CONTRAST statement in NLMIXED. Since there are three categories and two treatments, the interaction (regardless of time) is the cat 1 - cat 2 difference in trt 1 minus the cat 1 - cat 2 difference in trt 0 AND the cat 1 - cat 3 difference in trt 1 minus the cat 1 - cat 3 difference in trt 0. That is, the test of the trt*cat interaction is a joint test of two difference in differences. Instead of trying to directly use the model specification to write out this joint contrast in terms of the model parameters, you can get some help by using the LSMEANS and LSMESTIMATE statements in MIXED to tell you what the contrast coefficients should be. The LSMEANS statement below gives the LS-means for the six trt*cat combinations at time 2.

proc mixed;
class trt cat/ref=first;
model xu5=var0 trt|time|time cat|time|time trt|cat|time / solution;
lsmeans trt*cat / at time=2; 
run;

Using the order of the LS-means as shown, you can easily write the two difference in difference contrasts for the interaction. The E option shows the coefficients you specified as applied the LS-means. The ELSM option shows the resulting coefficients as applied to the model parameters. Note that trt*cat*time parameters are involved using the desired time=2 setting specified with the AT option. The JOINT option gives the overall test of trt*cat in time 2.

lsmestimate trt*cat 
   't1c1-t1c2 - (t0c1-t0c2)' 1 -1 0 -1 1 0 ,
   't1c1-t1c3 - (t0c1-t0c3)' 1 0 -1 -1 0 1 
   / e elsm at time=2 joint;

You can now use the ELSM coefficients to write an equivalent CONTRAST statement.

contrast 'trt*cat' 
  trt*cat 1 -1 0 -1 1 0 trt*cat*time 2 -2 0 -2 2 0, 
  trt*cat 1 0 -1 -1 0 1 trt*cat*time 2 0 -2 -2 0 2 / e;

Now you need to translate this into an equivalent CONTRAST statement in NLMIXED. A complication is that your NLMIXED model is written in full-rank reference coding as opposed to the nonfull-rank coding that MIXED uses for CLASS effects. But the parameters in the reference coded model just throw out all the zero parameters corresponding to reference levels. This means that the contrast above just needs to use the first two coefficients in each of trt*cat and trt*cat*time corresponding to the nonzero parameter estimates in those effects. So, the equivalent contrast with reference coding is simply

contrast 'trt*cat' 
  trt*cat 1 -1 trt*cat*time 2 -2 , 
  trt*cat 1 0  trt*cat*time 2  0 ;

So, finally back to your NLMIXED model. Rather than use time2, I've written the quadratic term as time*time and put the parameters for each model effect on separate lines. The above contrast is then written as the coefficients times the corresponding parameters (betas).

mu1 = beta0 + beta1*var0 + 
      beta2*(trt=1) + 
      beta3*time + beta4*time*time +
      beta5*time*(trt=1) + beta6*time*time*(trt=1) + 
      beta7*(Category=2) + beta8*(Category=1) + 
      beta9*time*(Category=2) + beta10*time*(Category=1) + 
      beta11*time*time*(Category=2) + beta12*time*time*(Category=1) + 
      beta13*(trt=1)*(Category=2) + beta14*(trt=1)*(Category=1) + 
      beta15*time*(trt=1)*(Category=2) + beta16*time*(trt=1)*(Category=1) ;

contrast 'trt*cat at t=2' -beta13+beta14+2*beta16-2*beta15 , beta14+2*beta16;

[kc]

Thank you so much for the detailed response @StatDave. 

[kc]

@StatDave:

Using the same equation below, could you please check if the following estimate statements are correct when estimating average estimate of the outcome when

1. treatment group = 1, time = 2, and category = 0

2. treatment group = 0, time = 2, and category = 0 and 

3. Difference of 1 - 2 above (mean difference in outcome between trt 1 - trt 0 in category = 0 at time = 2)

mu1 @ trt=1, category=0 and time=2
    = beta0 + beta1*var0 + 
      beta2*1 + 
      beta3*2+ beta4*2*2+
      beta5*2*1 + beta6*2*2*1;
mu1 @ trt=0, category=0 and time=2
    = beta0 + beta1*var0 + 
      beta3*2 + beta4*2*2 ;
Difference between treatment group 1 vs. 0 at category = 0 and time = 2 :
Diff = beta2*1 + beta5*2*1 + beta6*2*2*1;

[StatDave]

That just needs straightforward plug-in to the specified model, so yes, that looks right.