Mixed modelling: discrete vs continuous time

[Irene Bobevski]

Hello,

I have a mixed model (Stata 18) examining the effect of depression diagnosis (binary moderator) on social support (continuous independent variable) over four time points (late pregnancy, 6 months postpartum, 12 months postpartum, and when the child is aged 4 years) in a sample of 675 women. Some covariates are also included. I have created an interaction term between depression and time. I have run two separate models, one defining time as continuous, and the other with discrete time points. Time follows a non-linear trend, so I have used quadratic time in the continuous time model.

The two models lead to different conclusions, so I would appreciate your help to determine which model is more appropriate for my data.

Here is the output from the two models and the marginsplot. By the way, I did try to include a quadratic time term in the random part of the model as well in the continuous time model, but this model did not converge.

Continuous time model:

. mixed sse_task i.depressed1##c.time_c##c.time_c ///
>       c.epds_c ///
>           i.parity c.agec i.Uni ///
>           i.w1_empl_bin ///
>           i.ethnicity_mo_bin i.ethnicity_fa_bin ///
>           c.sles_events_c  ///
>           ecr_anx_avg_w1to4_c ecr_avoid_avg_w1to4_c  ///
>       || mpewsid: time_c, covariance(unstructured) reml        
 
Performing EM optimization ...
 
Performing gradient-based optimization:
Iteration 0:  Log restricted-likelihood = -4770.1932
Iteration 1:  Log restricted-likelihood = -4767.5348
Iteration 2:  Log restricted-likelihood = -4767.5289
Iteration 3:  Log restricted-likelihood = -4767.5289
 
Computing standard errors ...
 
Mixed-effects REML regression                        Number of obs    =  1,855
Group variable: mpewsid                              Number of groups =    675
                                                     Obs per group:
                                                                  min =      1
                                                                  avg =    2.7
                                                                  max =      4
                                                     Wald chi2(15)    = 390.69
Log restricted-likelihood = -4767.5289               Prob > chi2      = 0.0000
 
----------------------------------------------------------------------------------------------
                    sse_task | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-----------------------------+----------------------------------------------------------------
                1.depressed1 |   .2865289   .3794629     0.76   0.450    -.4572046    1.030262
                      time_c |  -1.546026   .1987641    -7.78   0.000    -1.935596   -1.156455
                             |
         depressed1#c.time_c |
                          1  |  -.6075963   .4814766    -1.26   0.207    -1.551273    .3360804
                             |
           c.time_c#c.time_c |   .3762828   .0669686     5.62   0.000     .2450268    .5075388
                             |
depressed1#c.time_c#c.time_c |
                          1  |   .0417941   .1595366     0.26   0.793    -.2708918    .3544801
                             |
                      epds_c |  -.1388265   .0184937    -7.51   0.000    -.1750735   -.1025795
                             |
                      parity |
         1. Not Nulliparous  |   .2275241   .2618764     0.87   0.385    -.2857442    .7407925
                        agec |   .0338147   .0264577     1.28   0.201    -.0180414    .0856709
                             |
                         Uni |
              1. University  |  -.0600735   .2522428    -0.24   0.812    -.5544603    .4343134
               1.w1_empl_bin |  -.0888899   .3328792    -0.27   0.789    -.7413212    .5635415
                             |
            ethnicity_mo_bin |
                   1. Other  |   .0570491   .3783086     0.15   0.880     -.684422    .7985203
                             |
            ethnicity_fa_bin |
              1. Other/ETSI  |    .419161   .2723208     1.54   0.124    -.1145779    .9528999
               sles_events_c |  -.0595556   .0432661    -1.38   0.169    -.1443556    .0252444
         ecr_anx_avg_w1to4_c |  -.0445726    .023182    -1.92   0.055    -.0900086    .0008633
       ecr_avoid_avg_w1to4_c |   -.261847   .0301575    -8.68   0.000    -.3209546   -.2027394
                       _cons |    14.1528   .3740304    37.84   0.000     13.41971    14.88589
----------------------------------------------------------------------------------------------
 
------------------------------------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
mpewsid: Unstructured        |
                 var(time_c) |   .3016405   .1238084      .1349316    .6743193
                  var(_cons) |    6.49939   .6623925      5.322571    7.936403
           cov(time_c,_cons) |  -.3605016    .231909     -.8150348    .0940317
-----------------------------+------------------------------------------------
               var(Residual) |   5.798138   .3019169      5.235586    6.421134
------------------------------------------------------------------------------
LR test vs. linear model: chi2(3) = 400.36                Prob > chi2 = 0.0000
 
Note: LR test is conservative and provided only for reference.
 
.
end of do-file

Discrete time model:

. mixed sse_task i.depressed1##i.time_c ///
>       c.epds_c ///
>           i.parity c.agec i.Uni ///
>           i.w1_empl_bin ///
>           i.ethnicity_mo_bin i.ethnicity_fa_bin ///
>           c.sles_events_c  ///
>           ecr_anx_avg_w1to4_c ecr_avoid_avg_w1to4_c ///
>       || mpewsid: time_c, covariance(unstructured) reml
 
Performing EM optimization ...
 
Performing gradient-based optimization:
Iteration 0:  Log restricted-likelihood = -4766.5115
Iteration 1:  Log restricted-likelihood = -4763.8273
Iteration 2:  Log restricted-likelihood =  -4763.821
Iteration 3:  Log restricted-likelihood =  -4763.821
 
Computing standard errors ...
 
Mixed-effects REML regression                        Number of obs    =  1,855
Group variable: mpewsid                              Number of groups =    675
                                                     Obs per group:
                                                                  min =      1
                                                                  avg =    2.7
                                                                  max =      4
                                                     Wald chi2(17)    = 391.93
Log restricted-likelihood =  -4763.821               Prob > chi2      = 0.0000
 
---------------------------------------------------------------------------------------
             sse_task | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
----------------------+----------------------------------------------------------------
         1.depressed1 |   .2949209   .3827446     0.77   0.441    -.4552446    1.045086
                      |
               time_c |
                   1  |  -1.273449   .1761843    -7.23   0.000    -1.618764   -.9281342
                   2  |  -1.520477   .1821599    -8.35   0.000    -1.877504   -1.163451
                   3  |  -1.307527   .2046438    -6.39   0.000    -1.708622   -.9064327
                      |
    depressed1#time_c |
                 1 1  |  -.6201199   .4151962    -1.49   0.135    -1.433889    .1936496
                 1 2  |  -.9784331   .4578381    -2.14   0.033    -1.875779   -.0810869
                 1 3  |  -1.477005   .4838893    -3.05   0.002    -2.425411   -.5285996
                      |
               epds_c |  -.1384583   .0185332    -7.47   0.000    -.1747826    -.102134
                      |
               parity |
  1. Not Nulliparous  |   .2303759   .2620228     0.88   0.379    -.2831794    .7439312
                 agec |   .0335779   .0264718     1.27   0.205    -.0183058    .0854616
                      |
                  Uni |
       1. University  |  -.0594466   .2523863    -0.24   0.814    -.5541146    .4352214
        1.w1_empl_bin |  -.0888475   .3330563    -0.27   0.790    -.7416259    .5639309
                      |
     ethnicity_mo_bin |
            1. Other  |   .0557079   .3785079     0.15   0.883     -.686154    .7975697
                      |
     ethnicity_fa_bin |
       1. Other/ETSI  |   .4190539    .272474     1.54   0.124    -.1149853    .9530931
        sles_events_c |  -.0543912   .0434807    -1.25   0.211    -.1396118    .0308293
  ecr_anx_avg_w1to4_c |  -.0445926   .0231945    -1.92   0.055    -.0900529    .0008678
ecr_avoid_avg_w1to4_c |  -.2614904   .0301753    -8.67   0.000    -.3206328   -.2023479
                _cons |   14.17102   .3747369    37.82   0.000     13.43655    14.90549
---------------------------------------------------------------------------------------
 
------------------------------------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
mpewsid: Unstructured        |
                 var(time_c) |   .2990392   .1237861      .1328546    .6731003
                  var(_cons) |   6.494387   .6626208      5.317289    7.932061
           cov(time_c,_cons) |  -.3536033   .2318515     -.8080239    .1008173
-----------------------------+------------------------------------------------
               var(Residual) |   5.801771     .30254      5.238099      6.4261
------------------------------------------------------------------------------
LR test vs. linear model: chi2(3) = 400.09                Prob > chi2 = 0.0000
 
Note: LR test is conservative and provided only for reference.
 
.
end of do-file

Marginsplot:


From the continuous time model I would conclude that there is no difference in the time trend between participants with and without depression, based on the non-significant interaction term of depression and quadratic time. But from the discrete time model the conclusion would be that compared to the baseline time, at the last two time points there is a difference between the two groups, with participants with depression having lower support rating than those without depression. This is based on the significant interaction of depression with the last two time points. Visually, the graph shows a relationship as described in my interpretation of the discrete time model, although the differences between the groups with and without depression at the last two time points are small.

Which model/interpretation would be more appropriate for this data, given the number of time points and the very uneven spacing between them?

[Clyde Schechter]

Well, I don't think a model using "continuous time" with a variable time_c that appears to have numerical values 0, 1, 2, 3 makes much sense when the actual times in question are 0 (or some small negative), 6, 12, and 48 months, regardless of whether a quadratic term is included or not.

If you did use the actual time values (0, 6, 12, and 48) and treated it as a continuous variable and included both linear and quadratic terms in both the fixed and random effects part, I think you might have a better chance of getting convergence if you centered that continuous time variable.

That said, looking at the output from the discrete time model, the results do not look at all quadratic. Your expected outcome variable, depressed or not, is 14 and change at time 0. It then drops down to about 12.8 give or take at time 6 months, descends a bit farther to 12.5-12.6 at 12 months and then rises a tad to around 12.8 at 48 months. These data are completely lacking the symmetry of a parabola around its vertex. So while the relationship to time is non-linear, it seems like a quadratic model would not describe it at all well. If any (semi-)parametric relationship were to be used for this, it would need to be more like a step function or perhaps a piecewise linear function. So the discrete model is probably a better bet, having the advantage of being a very close fit to the data and being very easy to describe to almost any audience.

Your margins plot was not displayable, at least not on my setup. This is often true of images posted on this forum, a reason why they are discouraged. Posting the table of -margins- output would have been more helpful.

[Irene Bobevski]

Thanks Clyde. I did try with irregular time intervals (0, 9, 15, 51 - as there are 9 months between late pregnancy 28 weeks and 6mths postpartum)) for the continuous time, as you suggested, but it still does not converge if the a quadratic time term is included in the random part:

mixed sse_task i.depressed1##c.time_irreg##c.time_irreg_c ///
>       c.epds_c ///
>           i.parity c.agec i.Uni ///
>           i.w1_empl_bin ///
>           i.ethnicity_mo_bin i.ethnicity_fa_bin ///
>           c.sles_events_c  ///
>           ecr_anx_avg_w1to4_c ecr_avoid_avg_w1to4_c  ///
>       || mpewsid: time_irreg_c##time_irreg_c, /// 
>           covariance(unstructured) reml   
......
Iteration 295: Log restricted-likelihood = -4765.5074  (backed up)
Iteration 296: Log restricted-likelihood = -4765.5074  (backed up)
Iteration 297: Log restricted-likelihood = -4765.5074  (backed up)
Iteration 298: Log restricted-likelihood = -4765.5074  (backed up)
Iteration 299: Log restricted-likelihood = -4765.5074  (backed up)
Iteration 300: Log restricted-likelihood = -4765.5074  (backed up)
convergence not achieved

Here is the margins table from the model with discrete time:

margins depressed1, at(time_irreg_c=(0 9 15 51)) atmeans


--

--------------------------------------------------------------------------------
               |            Delta-method
               |     Margin   std. err.      z    P>|z|     [95% conf. interval]
---------------+----------------------------------------------------------------
_at#depressed1 |
          1 0  |   14.24368   .1538221    92.60   0.000     13.94219    14.54516
          1 1  |     14.536    .351186    41.39   0.000     13.84769    15.22431
          2 0  |   12.97363   .1739217    74.59   0.000     12.63275    13.31451
          2 1  |    12.6449    .376004    33.63   0.000     11.90795    13.38186
          3 0  |   12.73601    .173371    73.46   0.000      12.3962    13.07581
          3 1  |   12.03612   .4060123    29.64   0.000     11.24035    12.83189
          4 0  |   12.94705   .1919876    67.44   0.000     12.57077    13.32334
          4 1  |   11.78416   .4265181    27.63   0.000      10.9482    12.62012
--------------------------------------------------------------------------------

I have tried to include the plot again as an attachment.

[Joseph Coveney]

I did try with irregular time intervals . . . for the continuous time, as you suggested, but it still does not converge if the a quadratic time term is included in the random part:

From your margins it doesn't look as if you centered the time variable. I second Clyde's other suggestion, too, and go with discrete time.

In lieu of trying to fit a model with discrete time in the fixed-effects equation and linear time as the random slope, I recommend fitting an unstructured residual covariance matrix. Something like the following.

mixed sse_task i.depressed1##i.time_irreg_c ///
    c.(epds_c agec sles_events_c ecr_anx_avg_w1to4_c ecr_avoid_avg_w1to4_c) ///
    i.(parity Uni w1_empl_bin ethnicity_??_bin) || mpewsid: , ///
        noconstant residuals(unstructured, t(time_irreg_c)) ml

I gather that these women's ages were recorded just once, at enrollment. Likewise for the diagnosis of depression.

[Irene Bobevski]

Thanks for your suggestions, Joseph.

From your margins it doesn't look as if you centered the time variable.

What do you mean about not centering my time variable? I though that I had centered it around the first time point: 0, 9 months, 15 monts, 51 months.

I gather that these women's ages were recorded just once, at enrollment. Likewise for the diagnosis of depression.

Yes, depression, age, and other demographic variables were only measured once at recruitment. Two of the continous covariates, epds_c and sles_events_c, are time variant, as they were mesured at each time point.

I second Clyde's other suggestion, too, and go with discrete time.

I agree with Clyde and you, the discrete time approach makes more sense with this data. I run it with unstructured residual covariance matrix, as you suggested. Does the model below look correct?

. mixed sse_task i.depressed1##i.time_irreg_c ///
>       c.epds_c ///
>           i.parity c.agec i.Uni ///
>           i.w1_empl_bin ///
>           i.ethnicity_mo_bin i.ethnicity_fa_bin ///
>           c.sles_events_c  ///
>           ecr_anx_avg_w1to4_c ecr_avoid_avg_w1to4_c  ///
>       || mpewsid: , /// 
>           noconstant ///
>           residuals(unstructured, t(time_irreg_c)) reml   

Performing gradient-based optimization: 
Iteration 0:  Log restricted-likelihood = -4963.8649  
Iteration 1:  Log restricted-likelihood = -4883.9944  
Iteration 2:  Log restricted-likelihood = -4845.2537  
Iteration 3:  Log restricted-likelihood = -4753.6429  
Iteration 4:  Log restricted-likelihood = -4744.6715  
Iteration 5:  Log restricted-likelihood = -4744.6073  
Iteration 6:  Log restricted-likelihood = -4744.6073  

Computing standard errors ...

Mixed-effects REML regression                        Number of obs    =  1,855
Group variable: mpewsid                              Number of groups =    675
                                                     Obs per group:
                                                                  min =      1
                                                                  avg =    2.7
                                                                  max =      4
                                                     Wald chi2(17)    = 385.62
Log restricted-likelihood = -4744.6073               Prob > chi2      = 0.0000

-----------------------------------------------------------------------------------------
               sse_task | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
------------------------+----------------------------------------------------------------
           1.depressed1 |   .2544759   .3610448     0.70   0.481    -.4531589    .9621108
                        |
           time_irreg_c |
                     9  |  -1.269773   .1777805    -7.14   0.000    -1.618216   -.9213295
                    15  |  -1.486523   .1780766    -8.35   0.000    -1.835547     -1.1375
                    51  |  -1.293413    .205566    -6.29   0.000    -1.696315   -.8905108
                        |
depressed1#time_irreg_c |
                  1  9  |  -.5919058   .4181723    -1.42   0.157    -1.411508    .2276969
                  1 15  |  -1.041509   .4461419    -2.33   0.020    -1.915931   -.1670865
                  1 51  |  -1.429514   .4874679    -2.93   0.003    -2.384934   -.4740946
                        |
                 epds_c |  -.1343405   .0183425    -7.32   0.000    -.1702911   -.0983899
                        |
                 parity |
    1. Not Nulliparous  |   .1020109   .2561787     0.40   0.690    -.4000901    .6041118
                   agec |   .0188553   .0259932     0.73   0.468    -.0320905    .0698011
                        |
                    Uni |
         1. University  |  -.1038651   .2471047    -0.42   0.674    -.5881813    .3804512
          1.w1_empl_bin |  -.0173155   .3254857    -0.05   0.958    -.6552557    .6206247
                        |
       ethnicity_mo_bin |
              1. Other  |   .0788357   .3685585     0.21   0.831    -.6435257    .8011971
                        |
       ethnicity_fa_bin |
         1. Other/ETSI  |   .4008662    .266423     1.50   0.132    -.1213133    .9230458
          sles_events_c |  -.0516442   .0424842    -1.22   0.224    -.1349117    .0316233
    ecr_anx_avg_w1to4_c |  -.0454129   .0227916    -1.99   0.046    -.0900836   -.0007422
  ecr_avoid_avg_w1to4_c |  -.2495812   .0296356    -8.42   0.000     -.307666   -.1914964
                  _cons |   14.18797   .3646515    38.91   0.000     13.47326    14.90267
-----------------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
mpewsid:             (empty) |
-----------------------------+------------------------------------------------
Residual: Unstructured       |
                     var(e0) |   10.88383   .6121461      9.747807    12.15224
                     var(e9) |   14.68343   .9909199      12.86423     16.7599
                    var(e15) |   12.53711    .870972      10.94115    14.36586
                    var(e51) |   11.89054   .8646098      10.31116    13.71185
                  cov(e0,e9) |   6.879468   .6360087      5.632914    8.126022
                 cov(e0,e15) |   5.588027   .6076663      4.397023    6.779032
                 cov(e0,e51) |   4.122202   .6017744      2.942746    5.301658
                 cov(e9,e15) |   8.647218   .7852754      7.108107    10.18633
                 cov(e9,e51) |   5.640203   .7617521      4.147196    7.133209
                cov(e15,e51) |   5.900792   .7039363      4.521102    7.280482
------------------------------------------------------------------------------
LR test vs. linear model: chi2(9) = 438.52                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.

. 
end of do-file

.

[Joseph Coveney]

What do you mean about not centering my time variable?

You can get a better idea by Googling multilevel modeling centering and taking a look at a few hits there. There is also some discussion about the relative merits to subtracting the grand mean (median) from everyone's values versus participantwise centering.

. . . the discrete time approach makes more sense with this data. I run it with unstructured residual covariance matrix, as you suggested. Does the model below look correct?

Well, the algorithm converged and you get coefficients and their standard errors for everything. So in that sense it appears to give what the command asks for. I don't feel qualified to speak to exogeneity of the predictors, whether the model adequately represents the data-generating process and is suitable for the intended purpose and so on.

A couple of minor observations:

1. Is there a particular reason why you used REML here? I would have been concerned about its limiting the likelihood-ratio testing / profile likelihood confidence intervals possible postestimation and I'm not sure of any compensatory advantage that REML confers in your case.

2. There is some temporal trail off evident in the correlation as expected

mata:
Cov = invvech((10.88383, 6.879468, 5.588027, 4.122202,
      14.68343, 8.647218, 5.640203,
      12.53711, 5.900792,
      11.89054)')

corr(Cov)

end

but the residual variance and covariance do look kind of similar between time points and you might be able to get away with a simpler residual covariance structure (exchangeable) if you wanted to keep things simple.

[Irene Bobevski]

Thanks for you help, Joseph. Regarding REML, I have used it because I thought that it offers an advantage with imbalanced designs, as we have different number of participants at different time points due to some missing data. Is this correct?

[Joseph Coveney]

Regarding REML, I have used it because I thought that it offers an advantage with imbalanced designs . . .

The only circumstance that I am aware of with that consideration is when one wants to fit a small-sample repeated-measures ANOVA or MANOVA with missing data (unbalanced datasets). In that case you cannot use conventional approaches, such as anova or manova unless one omits incomplete records by listwise deletion. But one can use maximum-likelihood mixed models, such as with mixed, but one needs to accommodate the small sample size with adjusted degrees of freedom methods, namely, Kenward-Roger,and that method requires that the model be fitted via REML.

Other than that, REML's benefit is for those interested in the variance components in cases where the sample size is small relative to the number of fixed effects. REML is supposed to have less downwardly biased estimates of the variance components in those cases. In your case, though, you have fewer than 20 fixed effects parameters to estimate and a large sample size (675), and so the advantage of REML is liable to be nil.

On another point, I don't want to put words into his mouth, but it struck me after posting is that I suspect that Clyde's recommendation about centering the time variable has more to do with alleviating the near-collinearity that often arises when forming quadratic polynomial predictors from continuous variables (the c.time##c.time expression). This near collinearity can result in numerical instability in the fitting algorithm and it will have difficulty in converging. You mentioned convergence failure when attempting to put the quadratic predictors in the random effects equation, and centering can help in just such cases. I'm not sure whether the Google search brought that point up.

[Clyde Schechter]

I don't want to put words into his mouth, but it struck me after posting is that I suspect that Clyde's recommendation about centering the time variable has more to do with alleviating the near-collinearity that often arises when forming quadratic polynomial predictors from continuous variables (the c.time##c.time expression).....

Yes, that's precisely what I had in mind.

[Irene Bobevski]

Thanks for your previous help with my model. I have a further question about my mixed effects model with discrete time points discussed previously in this post. Would it make sence to estimate random slopes for the discrete time? I know that there are no slopes with discrete time points, but I was interested in the variances of the random time ponts as an indicator of whether the effect of time at each point may be different for different groups of people (i.e. if the 95% CI of the variance is greater than 0). So would this model make sense?

mixed sse_task i.time_irreg_c##i.depressed1 ///
      c.epds_c ///
      c.agec i.Uni i.parity ///
      i.w1_empl_bin ///
      i.ethnicity_mo_bin i.ethnicity_fa_bin ///
      c.sles_events_c  ///
      ecr_anx_avg_w1to4_c ecr_avoid_avg_w1to4_c ///
      || mpewsid: i.time_irreg_c

[Clyde Schechter]

From a purely statistical perspective, it is possible to incorporate those random slopes. The values of those variances will tell you the extent to which different mpewsid's (? = study participants) vary in the time trajectory of outcome sse_task conditional on the values of your right hand side variables. Fine.

As to whether it makes sense to do so, that is a substantive question. What would those random slopes signify in your project? How would you explain and interpret those results? I do not have a sufficient understanding of your research context to answer those questions.

[Irene Bobevski]

We are interested in whether there is variation between individuals in how sse_task (partner support with tasks from pregnancy to 4 years postpartum) changes over time. This will indicate that there may be subgroups with different patterns of change over time.

[Clyde Schechter]

I don't see how random slopes will help you identify subgroups with different patterns of change over time. If you have some hypotheses about what those subgroups are, you can create a variable that defines them, and then add interactions between the subgroup indicators and the time indicators. Random slopes quantify the extent of individual (mpewsid)-level variation in how the outcome level is associated with each time period.

[Erik Ruzek]

Agreed with Clyde. To get subgroups without a priori ideas about group composition (i.e., some variable in your dataset that helps identifies the groups), you would have to use mixture modeling to identify the latent (unobserved) groups. See this article. These models can potentially be estimated with gsem, though I have not seen a worked example to point you to. The much earlier Stata command gllamm could almost certainly estimate such a model, but again, I do not have a good worked example. The other Stata possibility is the traj plugin by Jones and Nagin.