Mixed regression/growth model with a quadratic interaction term: how to obtain and interprete simple slopes

[Irene Bobevski]

Hello Statalisters,

I am running a growth model using mixed with long formatted data. I am examining the effect of IVF treatment (binary variable ivf = 0 or 1) on anxiety over four time points. The research question is whether women with and without IVF will have different anxiety trajectories over time. The sample size is N=764. I am using study_wave (with values 0, 1, 2, 3) as a continuous variable indicating the time points. Data was collected every six months. The outcome variable is anxiety (a continuous measure of anxiety). There are also two covariates in the model. Prior investigations showed that study_wave is not linear, so I have included its quadratic term, as well as an interaction of the quadratic term with ivf: i.ivf##c.study_wave##c.study_wave in order to compare the trajectories over time for ivf=0 (no IVF) and ivf=1 (IVF treatment received).

Here is the code and output, as well as the graph using marginsplot.

. mixed anxiety i.ivf##c.study_wave##c.study_wave i.Depressed ib(4).mod_4|| ///
>       mpewsid: c.study_wave, covariance(unstructured) reml        

Performing EM optimization ...

Performing gradient-based optimization:
Iteration 0:  Log restricted-likelihood = -8954.5523  
Iteration 1:  Log restricted-likelihood = -8953.0468  
Iteration 2:  Log restricted-likelihood = -8953.0439  
Iteration 3:  Log restricted-likelihood = -8953.0439  

Computing standard errors ...

Mixed-effects REML regression                        Number of obs    =  2,511
Group variable: mpewsid                              Number of groups =    764
                                                     Obs per group:
                                                                  min =      1
                                                                  avg =    3.3
                                                                  max =      4
                                                     Wald chi2(9)     = 166.49
Log restricted-likelihood = -8953.0439               Prob > chi2      = 0.0000

-----------------------------------------------------------------------------------------------
                      anxiety | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
------------------------------+----------------------------------------------------------------
                          ivf |
              1. IVF or GIFT  |  -3.755978    1.69605    -2.21   0.027    -7.080175   -.4317808
                   study_wave |  -1.545982   .4163235    -3.71   0.000    -2.361961   -.7300031
                              |
             ivf#c.study_wave |
              1. IVF or GIFT  |   4.937875   1.776305     2.78   0.005     1.456381    8.419369
                              |
    c.study_wave#c.study_wave |   .3376616   .1371033     2.46   0.014      .068944    .6063792
                              |
ivf#c.study_wave#c.study_wave |
              1. IVF or GIFT  |  -1.410498   .5727671    -2.46   0.014    -2.533101   -.2878951
                              |
                    Depressed |
                1. Depressed  |   8.884062   .7758624    11.45   0.000       7.3634    10.40472
                              |
                        mod_4 |
                      1. NVD  |  -1.313946   1.074189    -1.22   0.221    -3.419317    .7914254
          2.  NVD - assisted  |  -2.590648   1.165324    -2.22   0.026    -4.874642   -.3066545
             3. Emergency CS  |  -1.680633   1.149658    -1.46   0.144     -3.93392    .5726549
                              |
                        _cons |   34.92902   1.024893    34.08   0.000     32.92026    36.93777
-----------------------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
mpewsid: Unstructured        |
             var(study_wave) |   3.128821   .8150812      1.877744    5.213449
                  var(_cons) |   68.30227   5.326396      58.62139    79.58187
       cov(study_wave,_cons) |   -7.12071   1.719405     -10.49068   -3.750739
-----------------------------+------------------------------------------------
               var(Residual) |   41.56179   1.745478      38.27773    45.12761
------------------------------------------------------------------------------
LR test vs. linear model: chi2(3) = 738.77                Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.
  
.
end of do-file

margins ivf, at(study_wave=(0 1 2 3)) atmeans
marginsplot, x(study_wave) noci

I have searched for answers to my questions, but I am still not clear how to interpret a quadratic interaction term. So here are the questions.
Question 1: I have used margins ivf, dydx(study_wave) atmeans to obtain what I understand are simple slopes for ivf=0 and ivf=1. My understanding is that these are the simple slopes based both the quadratic and linear trends. Is this correct?

Here is the code and output:

. margins ivf, dydx(study_wave) atmeans    

Conditional marginal effects                             Number of obs = 2,511

Expression: Linear prediction, fixed portion, predict()
dy/dx wrt:  study_wave
At: 0.ivf       = .9438471 (mean)
    1.ivf       = .0561529 (mean)
    study_wave  = 1.326961 (mean)
    0.Depressed = .8191955 (mean)
    1.Depressed = .1808045 (mean)
    1.mod_4     = .4532059 (mean)
    2.mod_4     =  .215452 (mean)
    3.mod_4     = .2345679 (mean)
    4.mod_4     = .0967742 (mean)

---------------------------------------------------------------------------------
                |            Delta-method
                |      dy/dx   std. err.      z    P>|z|     [95% conf. interval]
----------------+----------------------------------------------------------------
study_wave      |
            ivf |
         0. No  |  -.6498543   .1468912    -4.42   0.000    -.9377558   -.3619528
1. IVF or GIFT  |   .5446678   .5929341     0.92   0.358    -.6174616    1.706797
---------------------------------------------------------------------------------

.
end of do-file

Question 2: I think that it does not make sense to interpret the overall simple slope for ivf=1, as it is not linear, and different parts of the slope need to be considered separately. Is this correct?
However, for ivf=0 the slope seems approximately linear, so my interpretation is that overall there is a significant decrease in anxiety over time for ivf=0, given the negative significant slope. Is this correct?

Question 3: For the slope of ivf=1, my main interest is in this group and how it differs from ivf=0 at different time points, so I am using pwcompare to examine the four estimated mean differences between ivf=0 and ivf=1 at each time point. Is this correct?

. margins ivf, pwcompare(effects) at(study_wave=(0 1 2 3)) atmeans

So, from the above results my conclusion is that the anxiety trajectories for ivf=0 and ivf=1 follow a different pattern over time. Overall ivf=0 decreased over time. However, ivf=1 started off lower on anxiety than ivf=0 at study_wave 0, but then increased to the same level as ivf=0 between study_wave 1 and 2, and then again was lower than ivf=0 at study_wave 3. The mean differences were significant at study_wave 0 and 3 (from the above code).

[Clyde Schechter]

Before we even try to interpret the results, you have to fix the model. The model you have created includes a random slope at the mpewsid level for the linear term in wave but not for the quadratic. So you are constraining the quadratic coefficient to be constant across participants while having the linear coefficient and the constant to vary among study subjects. This produces a rather strange kind of variation across subjects in the relationship between. To get a flavor of the kind of variation this corresponds to, look at the graph produced by this code (with three quadratics all having the same quadratic coefficient but varying linear and constant terms):

graph twoway function y = x^2 + x, range(0 3) || function y = 2*x^2 + x + .5, range(0 3) || function y = -4*x^2 + x +2, range(0 3)

I think you will agree that this is not a plausible model of variation. So to pursue this analysis, re-do the mixed model adding a random coefficient for the quadratic term.

Beyond that, it doesn't make sense to talk about a "simple slope" in a quadratic model. Quadratic functions do not have "a slope." Their slope ranges from negative to positive infinity according to the value of x (wave, in your case). The marginal effect of the linear term that you are calculating is just not interpretable as a simple slope, or as a slope at all. In a quadratic function, the linear term coefficient is really about the horizontal location of the parabola, and the slope at any given value of x is at least as strongly affected by the quadratic coefficient.

Then there is the issue of whether it makes sense to fit a parabola when you only have four time points.

Frankly, I think you would be much better off modeling time as a discrete variable and just study the differences between the two study groups at each of the four time points.