How to interpret results from a model with quadratic trends

[Urmeli Joost]

Hi all,

I am using Stata 14. I am trying to analyze longitudinal data from a cohort study where the same subjects are measured at three different time points (ages 15, 18 and 25 years). To allow for departures from linearity, I am considering a model with quadratic trends. Unfortunately I cannot interpret the results. Here is the code and results:

WHR - waist to hip ratio
kohort - birth cohort
aeg - time
sugu - sex
kood - ID

mixed WHR i.kohort##c.aeg##c.aeg if sugu ==2|| kood: , noconst residuals(unstructured, t(aeg)) stddev reml

Obtaining starting values by EM:

Performing gradient-based optimization:

Iteration 0:   log restricted-likelihood =  2621.9848  (not concave)
Iteration 1:   log restricted-likelihood =   2817.847  
Iteration 2:   log restricted-likelihood =  2854.2901  
Iteration 3:   log restricted-likelihood =   2858.196  
Iteration 4:   log restricted-likelihood =  2858.2035  
Iteration 5:   log restricted-likelihood =  2858.2035  

Computing standard errors:

Mixed-effects REML regression                   Number of obs     =      1,602
Group variable: kood                            Number of groups  =        655

                                                Obs per group:
                                                              min =          1
                                                              avg =        2.4
                                                              max =          3

                                                Wald chi2(5)      =     248.23
Log restricted-likelihood =  2858.2035          Prob > chi2       =     0.0000

------------------------------------------------------------------------------------
               WHR |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------------+----------------------------------------------------------------
          2.kohort |    .134359     .05627     2.39   0.017     .0240719    .2446461
               aeg |   .0238618   .0041599     5.74   0.000     .0157085    .0320152
                   |
      kohort#c.aeg |
                2  |  -.0170974   .0058537    -2.92   0.003    -.0285705   -.0056243
                   |
       c.aeg#c.aeg |  -.0005523    .000103    -5.36   0.000     -.000754   -.0003505
                   |
kohort#c.aeg#c.aeg |
                2  |   .0004883    .000145     3.37   0.001     .0002042    .0007725
                   |
             _cons |   .4994664   .0400719    12.46   0.000      .420927    .5780058
------------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
kood:                (empty) |
-----------------------------+------------------------------------------------
Residual: Unstructured       |
                     sd(e15) |   .0370384   .0010755      .0349892    .0392075
                     sd(e18) |   .0444294   .0013792      .0418069    .0472165
                     sd(e25) |    .056315   .0017583      .0529721    .0598688
               corr(e15,e18) |   .5656019   .0306453      .5025527    .6226668
               corr(e15,e25) |   .4291414   .0360977      .3558472    .4971886
               corr(e18,e25) |   .6089058   .0305249      .5456038    .6652989
------------------------------------------------------------------------------
LR test vs. linear model: chi2(5) = 472.44                Prob > chi2 = 0.0000

Note: The reported degrees of freedom assumes the null hypothesis is not on the boundary of the
      parameter space.  If this is not true, then the reported test is conservative.

If anyone can help or point to a good reference (have not found anything appropriate this far), I would be eternally grateful!

Urmeli

[Clyde Schechter]

Before interpreting the model, I would like to know a bit more about how this particular model was arrived at. If I understand correctly, you have only 3 time points, and they are exactly the same for all participants: 15, 18, and 25. While it is perfectly legal to model this as a continuous variable and fit parabolas through the three time points, it is unusual to do so. You can always exactly fit a parabola through three points, but more generally, quadratic models are not usually good descriptions of real world data generating processes. Is there a strong theoretical basis for thinking that the relationship of outcome to time more generally should be quadratic? If you are just thinking that you expect the outcome to go up and then decline (or the other way around), why not just model the three timepoints as a discrete variable? The results would be simpler to interpret and not claim anything about other points in time that were not observed. Also legal, but unusual, is the -noconst- option in the random effects model. Is there some specific reason to believe that in each cohort, all the participants are starting out from the same baseline value of WHR? I ask, because that is a constraint that your model is imposing.

Correction: You are not imposing a constraint that all participants in a cohort are starting from the same value of WHR. That is a mis-statement on my part. But you are instead specifically failing to capture that person-level baseline variance separately and are ascribing it entirely to noise. This seems strange for a variable like WHR, where, I would think there are many person-level variables that affect it, and whose effects you would want to dissect away from noise. So is there some specific reason you wanted to do that? It looks like you are trying to do a variant of difference-in-differences estimation here, and doing that requires a fixed or random effect at the individual (kood) level so you can clearly and completely separate between-person variation at baseline from within-person effects of cohort membership.

[Urmeli Joost]

Hi Clyde,

Thank you for your response, I really appreciate it. When I started to visually examine the data (using the "xtgraph" command) I realized that the changes in the mean response over time were not linear (in one or both of the birth cohorts, depending on the variable of interest) (see examples below). To be honest, this was the first time I have encountered this so I am actually not that confident about my models (i.e., the points that you have brought out seem pretty reasonable).

As far as I understood, I had two options:
a) fit a model with quadratic trends;
b) use linear splines.

As you can see, I ended up choosing the quadratic trends. I am not actually sure how to approach this, considering I have only 3 data points.

All feedback and thoughts are welcome. And directions to any good reasources/materials if my question is actually too basic (without me realizing it).

Thank you again,
Urmeli

[Clyde Schechter]

Well, if I were doing this, I would not try to specify a quadratic, or any other parametric model for the trajectory here. Even if it were conceptually correct for WHR,it might not be for Wmaxkg or other outcomes you might look at. Just using three discrete timepoints, when that is the way the data are structured, is much simpler and is free of any assumptions about the shape of the trajectory. (With only three time points, the discrete timepoint approach is equivalent to linear splines, but simpler to do and interpret.) So I would be running it as:

mixed WHR i.kohort##i.aeg if sugu ==2|| kood:, residuals(unstructured, t(aeg)) stddev reml
margins aeg#kohort
marginsplot

Note that I have also taken out the -noconst- option from the random effect as I don't see any rationale for it and you have not explained why you included it.

-margins aeg- will report the expected values of WHR in each cohort at each time period, and -marginsplot- will graph them for you.

[Urmeli Joost]

Thank you, Clyde, for the recommendation. It is highly appreciated!

Urmeli