piecewise growth curve

[carole fantini]

Hi everybody,
I try to model a growth curve for my data.
I try to observe differences in the trajectory of depression regarding the genetic status a participants (0-carier or 1-not carier). depression was measured at baseline (1), 2 month after result (2), 6 month (3), 1 (4) and 2 years (5) after results. I have read all is possible on LGC, and regarding the non linearity of my data, a piecewise model is the best solution.
I have tried tu run such a model but i'm clearly not sure about the syntax i have used. I'm interesting about the slope change between time intervals, and particularly if there is differences in the slope changes between carier and not carier of the gene.
Unfortunately, i have a troubling result.
the variables are :
depression : continuous
gene status : carier (1), non carier (0): categorical
Time : 1, 2, 3, 4 &5 : categorical

first i use mkspline to define 4 knots :

mkspline time1 2 time2 3 time3 4 time4 = Time, marginal

second, I ran the model with xtmixed

xtmixed dep porteur##time1 porteur##time2 porteur##time3 porteur##time4 || ID: Time, mle

I obtain the following estimates with warning message of colinearity.

. xtmixed dep porteur##time1 porteur##time2 porteur##time3 porteur##time4 || ID: Time, mle
note: 1.time2 omitted because of collinearity
note: 2.time2 omitted because of collinearity
note: 3.time2 omitted because of collinearity
note: 1.porteur#1.time2 omitted because of collinearity
note: 1.porteur#2.time2 omitted because of collinearity
note: 1.porteur#3.time2 omitted because of collinearity
note: 1.time3 omitted because of collinearity
note: 2.time3 omitted because of collinearity
note: 1.porteur#1.time3 omitted because of collinearity
note: 1.porteur#2.time3 omitted because of collinearity
note: 1.time4 omitted because of collinearity
note: 1.porteur#1.time4 omitted because of collinearity

Performing EM optimization:

Performing gradient-based optimization:

Iteration 0: log likelihood = -2278.8898
Iteration 1: log likelihood = -2277.0752
Iteration 2: log likelihood = -2277.0463
Iteration 3: log likelihood = -2277.0456
Iteration 4: log likelihood = -2277.0456

Computing standard errors:

Mixed-effects ML regression Number of obs = 631
Group variable: ID Number of groups = 193

Obs per group:
min = 1
avg = 3.3
max = 5

Wald chi2(9) = 51.69
Log likelihood = -2277.0456 Prob > chi2 = 0.0000

-------------------------------------------------------------------------------
dep | Coef. Std. Err. z P>|z| [95% Conf. Interval]
--------------+----------------------------------------------------------------
1.porteur | -1.000888 1.544808 -0.65 0.517 -4.028657 2.026881
|
time1 |
2 | -2.529787 1.084861 -2.33 0.020 -4.656075 -.4034987
3 | 1.761583 1.164203 1.51 0.130 -.5202119 4.043378
4 | 4.017521 1.201632 3.34 0.001 1.662365 6.372676
5 | -.4507877 1.278628 -0.35 0.724 -2.956852 2.055276
|
porteur#time1 |
1 2 | 4.780729 1.634244 2.93 0.003 1.577671 7.983788
1 3 | -3.651737 1.792905 -2.04 0.042 -7.165766 -.137707
1 4 | .4716759 1.868767 0.25 0.801 -3.19104 4.134392
1 5 | 1.41006 2.004701 0.70 0.482 -2.519083 5.339202
|
time2 |
1 | 0 (omitted)
2 | 0 (omitted)
3 | 0 (omitted)
|
porteur#time2 |
1 1 | 0 (omitted)
1 2 | 0 (omitted)
1 3 | 0 (omitted)
|
time3 |
1 | 0 (omitted)
2 | 0 (omitted)
|
porteur#time3 |
1 1 | 0 (omitted)
1 2 | 0 (omitted)
|
1.time4 | 0 (omitted)
|
porteur#time4 |
1 1 | 0 (omitted)
|
_cons | 12.6734 1.016518 12.47 0.000 10.68106 14.66574
-------------------------------------------------------------------------------

------------------------------------------------------------------------------
Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
-----------------------------+------------------------------------------------
ID: Independent |
sd(Time) | .2596632 1.007456 .0001294 521.1367
sd(_cons) | 7.785805 .5415579 6.793544 8.922994
-----------------------------+------------------------------------------------
sd(Residual) | 7.059303 .2545035 6.577701 7.576166
------------------------------------------------------------------------------
LR test vs. linear model: chi2(2) = 214.80 Prob > chi2 = 0.0000

I think that the time values in variable colons are curious. for example there is a 5 under time1 section. i don't understand this labelling.
second, I suppose that terms omitted are just redundant with the porteur#time2 for example, as time1-2 represent the slope for non cariers and porteur#time2 represent the slope for cariers.
what we see here is that slope changes between the previous knots are sometimes significant. for example, time1-2 is significant as well as porteur#time1-12. but that does not tell me if the two slopes changes are statistically and significantly different between each other.
further, graphing this is verry challenging and I don't have success to plot the two growth curves.
any help and feedback would be very helpful. I give you also some data...
thanks a lot
carole

[Clyde Schechter]

There are two mistakes, one of which is causing the colinearity, and the other you have not yet tripped over but it invalidates your model.

The first mistake is the you have treated your time1 through time4 variables as discrete, but they should be treated as continuous. So, it should be porteur##c.time1, etc.

The second mistake is that you are asking Stata to model random slopes on time, even though time is not included as a fixed effect. The result is that the random slopes will be calculated as if the mean random slope for time is zero--which really defeats the purpose of a growth model! You need to use the spline variables for random effects. So it should be:

mixed dep i.porteur##c.(time1 time2 time3 time4) || ID: time1 time2 time3 time4, mle

Not tested. Although you show example data, you did so in a screenshot, which the FAQ specifically asks people not to do. There is no way to import data from a screenshot to Stata. So I could not test this code. It is probably correct, although I have not thought through sufficiently whether a model with all of the time splines as random slopes is identifiable. If the estimation does not converge, consider dropping the random slope on one of the time spline variables.

Note that the -xtmixed- command was renamed -mixed- as of version 14. While the original name is still recognized, and is even treated slightly differently (reporting variance components as standard deviations instead of variances by default), at some point it may not be. Best to get in the habit of using the up-to-date name now.

In the future, when showing data examples, please use the -dataex- command to do so. If you are running version 15.1 or a fully updated version 14.2, it is already part of your official Stata installation. If not, run -ssc install dataex- to get it. Either way, run -help dataex- to read the simple instructions for using it. -dataex- will save you time; it is easier and quicker than typing out tables. It includes complete information about aspects of the data that are often critical to answering your question but cannot be seen from tabular displays or screenshots. It also makes it possible for those who want to help you to create a faithful representation of your example to try out their code, which in turn makes it more likely that their answer will actually work in your data.

[carole fantini]

Hi Clyde
Thank you so much for your feedback. I was also intriguing by the random structure with Time as continuous but I have follow the very excellent book of Sophia Rabe-Esketh with an example of piecewise HLM. but I probably have misunderstood something as she use the same syntax, I have reproduced with the data example she used. HLM is often very complex and subtle, particularly with growth curve.
I have used dataex to produce a useful file with data, joined to my message.
I don' t understand why we need to consider that time1, time2...are continuous because they took only one value, resulting from the mkspline command?
here is what i obtained with your syntax

 mixed dep i.porteur##c.(time1 time2 time3 time4) || ID: time1 time2 time3 time4, mle

Performing EM optimization: 

Performing gradient-based optimization: 

Iteration 0:   log likelihood =  -2297.211  
Iteration 1:   log likelihood = -2276.9154  (not concave)
Iteration 2:   log likelihood = -2276.3253  
Iteration 3:   log likelihood = -2275.5142  
Iteration 4:   log likelihood = -2275.4972  
Iteration 5:   log likelihood =  -2275.497  

Computing standard errors:

Mixed-effects ML regression                     Number of obs     =        631
Group variable: ID                              Number of groups  =        193

                                                Obs per group:
                                                              min =          1
                                                              avg =        3.3
                                                              max =          5

                                                Wald chi2(9)      =      53.55
Log likelihood =  -2275.497                     Prob > chi2       =     0.0000

---------------------------------------------------------------------------------
            dep |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
----------------+----------------------------------------------------------------
      1.porteur |  -5.732588   2.631107    -2.18   0.029    -10.88946   -.5757127
          time1 |  -2.488448   1.051253    -2.37   0.018    -4.548866    -.428029
          time2 |   6.776275    1.91673     3.54   0.000     3.019554      10.533
          time3 |  -2.016079   2.121262    -0.95   0.342    -6.173675    2.141518
          time4 |  -6.681047   2.237862    -2.99   0.003    -11.06718   -2.294918
                |
porteur#c.time1 |
             1  |   4.737828   1.583423     2.99   0.003     1.634376    7.841281
                |
porteur#c.time2 |
             1  |  -13.11178   2.918263    -4.49   0.000    -18.83147   -7.392086
                |
porteur#c.time3 |
             1  |   12.41117   3.295531     3.77   0.000     5.952046    18.87029
                |
porteur#c.time4 |
             1  |  -3.049441   3.536357    -0.86   0.389    -9.980574    3.881692
                |
          _cons |   15.15541   1.740362     8.71   0.000     11.74436    18.56646
---------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
ID: Independent              |
                  var(time1) |   4.71e-12          .             .           .
                  var(time2) |   6.66e-09          .             .           .
                  var(time3) |   4.847292          .             .           .
                  var(time4) |   2.05e-11          .             .           .
                  var(_cons) |   62.03209          .             .           .
-----------------------------+------------------------------------------------
               var(Residual) |   46.66183          .             .           .
------------------------------------------------------------------------------
LR test vs. linear model: chi2(5) = 217.90                Prob > chi2 = 0.0000

If i have well understood, the interpretation is that non carier at time 1 show a decrease (-2.49) of depression scores and there is a significant interaction as the estimate for carier (porteur in the table) at time 1 is significant. right?
is there a way to plot the graph with predicted scores in order to show the growth curves for carier and non carier? not sure that predict xb is adequat with a piecewise HLM (wihiwh is not linear by definition), nor the command lincom to compare coefficients for example time1 versus time 2 (even if these segments should be considered as linear)...

Thanks a lot Clyde

* Example generated by -dataex-. To install: ssc install dataex
clear
input double(dep porteur) byte(time1 time2 time3 time4) double Time
 1 1 1 0 0 0 1
 0 1 2 0 0 0 2
 0 1 3 1 0 0 3
 6 1 4 2 1 0 4
 0 1 5 3 2 1 5
 9 0 1 0 0 0 1
 0 0 2 0 0 0 2
 6 0 3 1 0 0 3
 6 0 4 2 1 0 4
 1 0 5 3 2 1 5
 1 1 1 0 0 0 1
 3 1 2 0 0 0 2
 . 1 3 1 0 0 3
 . 1 4 2 1 0 4
 . 1 5 3 2 1 5
33 0 1 0 0 0 1
24 0 2 0 0 0 2
 . 0 3 1 0 0 3
 . 0 4 2 1 0 4
 . 0 5 3 2 1 5
 4 0 1 0 0 0 1
15 0 2 0 0 0 2
 9 0 3 1 0 0 3
21 0 4 2 1 0 4
 . 0 5 3 2 1 5
 . 0 1 0 0 0 1
 . 0 2 0 0 0 2
 . 0 3 1 0 0 3
 . 0 4 2 1 0 4
 . 0 5 3 2 1 5
 3 1 1 0 0 0 1
12 1 2 0 0 0 2
27 1 3 1 0 0 3
 . 1 4 2 1 0 4
 . 1 5 3 2 1 5
27 1 1 0 0 0 1
27 1 2 0 0 0 2
22 1 3 1 0 0 3
36 1 4 2 1 0 4
 . 1 5 3 2 1 5
19 1 1 0 0 0 1
32 1 2 0 0 0 2
27 1 3 1 0 0 3
26 1 4 2 1 0 4
 . 1 5 3 2 1 5
 0 0 1 0 0 0 1
 3 0 2 0 0 0 2
 8 0 3 1 0 0 3
 6 0 4 2 1 0 4
 4 0 5 3 2 1 5
 0 0 1 0 0 0 1
 1 0 2 0 0 0 2
 4 0 3 1 0 0 3
 . 0 4 2 1 0 4
 . 0 5 3 2 1 5
19 1 1 0 0 0 1
21 1 2 0 0 0 2
 . 1 3 1 0 0 3
 . 1 4 2 1 0 4
 . 1 5 3 2 1 5
 2 0 1 0 0 0 1
 4 0 2 0 0 0 2
11 0 3 1 0 0 3
20 0 4 2 1 0 4
13 0 5 3 2 1 5
18 0 1 0 0 0 1
15 0 2 0 0 0 2
45 0 3 1 0 0 3
32 0 4 2 1 0 4
 7 0 5 3 2 1 5
15 0 1 0 0 0 1
 9 0 2 0 0 0 2
 . 0 3 1 0 0 3
15 0 4 2 1 0 4
14 0 5 3 2 1 5
 7 1 1 0 0 0 1
11 1 2 0 0 0 2
 9 1 3 1 0 0 3
 . 1 4 2 1 0 4
 8 1 5 3 2 1 5
 9 1 1 0 0 0 1
 1 1 2 0 0 0 2
 2 1 3 1 0 0 3
 8 1 4 2 1 0 4
 8 1 5 3 2 1 5
 8 1 1 0 0 0 1
11 1 2 0 0 0 2
15 1 3 1 0 0 3
16 1 4 2 1 0 4
11 1 5 3 2 1 5
 5 0 1 0 0 0 1
 4 0 2 0 0 0 2
 3 0 3 1 0 0 3
 7 0 4 2 1 0 4
 3 0 5 3 2 1 5
 2 0 1 0 0 0 1
 2 0 2 0 0 0 2
 3 0 3 1 0 0 3
 . 0 4 2 1 0 4
 . 0 5 3 2 1 5
end

[Clyde Schechter]

I don' t understand why we need to consider that time1, time2...are continuous because they took only one value, resulting from the mkspline command?

That's not true! Just take a look at your data and you will see that they take on several values and if you treat them as discrete variables you create the colinearity that causes your problem.

The -dataex- you sent cannot possibly represent the data on which you ran the code, because you do not have an ID variable, without which the code cannot ran. Your actual data presumably contains that variable. The best way to interpret your outputs and graph them is to follow the -mixed- command with:

forvalues t = 1/5 {
    local at`t'
    forvalues i = 1/4 {
        summ time`i' if Time == `t', meanonly
        local at`t' `at`t'' time`i' = `r(mean)'
    }
    display "at(`at`t'')"
}

tempfile results
margins porteur, at(`at1') at(`at2') at(`at3') at(`at4') at(`at5') saving(`results')
use `results', clear
keep _m1 _at _margin _ci*
rename _margin dep
rename _at time
reshape wide dep _ci*, i(time) j(_m1)
label var dep0 "Non-porteur"
label var dep1 "Porteur"
graph twoway (connect dep* time) (rcap _ci_lb0 _ci_ub0 time, legend(off)) ///
    (rcap _ci_lb1 _ci_ub1 time, legend(off))

The code is complicated due to the use of splines, which -margins- and -marginsplot- do not "play with" nicely.

If i have well understood, the interpretation is that non carier at time 1 show a decrease (-2.49) of depression scores and there is a significant interaction as the estimate for carier (porteur in the table) at time 1 is significant. right?

Yes, this is correct. But I don't think it's the best way to think about it or present it. Interactions are not intuitive. Nobody but statisticians really understand them. What makes more sense, actually, is to show the graph created by the code above, and in the table output by -margins-. Anybody can look at that and see what is going on.

[carole fantini]

Hi Clyde
thanks again :-) Indeed the code for graph is very complicated. not sure i will be able to write that by myself !!! :-)
here is a correct set of data.

* Example generated by -dataex-. To install: ssc install dataex
clear
input str6 ID double(dep porteur) byte(time1 time2 time3 time4) double Time
"4050" 25 1 1 0 0 0 1
"4050" 32 1 2 0 0 0 2
"4050"  . 1 3 1 0 0 3
"4050"  . 1 4 2 1 0 4
"4050"  . 1 5 3 2 1 5
"4111"  0 0 1 0 0 0 1
"4111"  0 0 2 0 0 0 2
"4111"  4 0 3 1 0 0 3
"4111"  9 0 4 2 1 0 4
"4111"  . 0 5 3 2 1 5
"4118" 12 0 1 0 0 0 1
"4118"  . 0 2 0 0 0 2
"4118"  . 0 3 1 0 0 3
"4118"  . 0 4 2 1 0 4
"4118"  . 0 5 3 2 1 5
"4120" 12 0 1 0 0 0 1
"4120" 15 0 2 0 0 0 2
"4120" 23 0 3 1 0 0 3
"4120" 15 0 4 2 1 0 4
"4120" 11 0 5 3 2 1 5
"4121"  0 1 1 0 0 0 1
"4121"  0 1 2 0 0 0 2
"4121"  6 1 3 1 0 0 3
"4121"  . 1 4 2 1 0 4
"4121"  . 1 5 3 2 1 5
"4124"  . 1 1 0 0 0 1
"4124"  . 1 2 0 0 0 2
"4124"  . 1 3 1 0 0 3
"4124"  . 1 4 2 1 0 4
"4124"  . 1 5 3 2 1 5
"4128"  7 1 1 0 0 0 1
"4128"  . 1 2 0 0 0 2
"4128"  . 1 3 1 0 0 3
"4128"  . 1 4 2 1 0 4
"4128"  . 1 5 3 2 1 5
"4129"  . 0 1 0 0 0 1
"4129"  6 0 2 0 0 0 2
"4129"  . 0 3 1 0 0 3
"4129"  . 0 4 2 1 0 4
"4129"  . 0 5 3 2 1 5
"4130"  . 0 1 0 0 0 1
"4130" 16 0 2 0 0 0 2
"4130" 19 0 3 1 0 0 3
"4130"  . 0 4 2 1 0 4
"4130" 26 0 5 3 2 1 5
"4131"  9 0 1 0 0 0 1
"4131"  . 0 2 0 0 0 2
"4131"  . 0 3 1 0 0 3
"4131"  . 0 4 2 1 0 4
"4131"  . 0 5 3 2 1 5
"4134" 14 1 1 0 0 0 1
"4134"  . 1 2 0 0 0 2
"4134"  . 1 3 1 0 0 3
"4134"  . 1 4 2 1 0 4
"4134"  . 1 5 3 2 1 5
"4135"  4 1 1 0 0 0 1
"4135"  2 1 2 0 0 0 2
"4135"  2 1 3 1 0 0 3
"4135"  6 1 4 2 1 0 4
"4135"  5 1 5 3 2 1 5
"4136"  5 1 1 0 0 0 1
"4136"  4 1 2 0 0 0 2
"4136"  3 1 3 1 0 0 3
"4136" 11 1 4 2 1 0 4
"4136"  5 1 5 3 2 1 5
"4144" 14 1 1 0 0 0 1
"4144" 33 1 2 0 0 0 2
"4144"  9 1 3 1 0 0 3
"4144" 17 1 4 2 1 0 4
"4144"  9 1 5 3 2 1 5
"4145"  5 0 1 0 0 0 1
"4145"  0 0 2 0 0 0 2
"4145"  6 0 3 1 0 0 3
"4145"  6 0 4 2 1 0 4
"4145"  9 0 5 3 2 1 5
"4146" 13 0 1 0 0 0 1
"4146" 11 0 2 0 0 0 2
"4146"  . 0 3 1 0 0 3
"4146"  . 0 4 2 1 0 4
"4146"  . 0 5 3 2 1 5
"4148"  3 0 1 0 0 0 1
"4148" 14 0 2 0 0 0 2
"4148"  4 0 3 1 0 0 3
"4148" 10 0 4 2 1 0 4
"4148"  3 0 5 3 2 1 5
"4149" 28 0 1 0 0 0 1
"4149"  . 0 2 0 0 0 2
"4149"  . 0 3 1 0 0 3
"4149"  . 0 4 2 1 0 4
"4149"  . 0 5 3 2 1 5
"4150" 19 1 1 0 0 0 1
"4150" 22 1 2 0 0 0 2
"4150" 18 1 3 1 0 0 3
"4150" 22 1 4 2 1 0 4
"4150" 14 1 5 3 2 1 5
"4151" 11 0 1 0 0 0 1
"4151" 10 0 2 0 0 0 2
"4151" 30 0 3 1 0 0 3
"4151"  9 0 4 2 1 0 4
"4151"  7 0 5 3 2 1 5
end

------------------ copy up to and including the previous line ------------------

[carole fantini]

Hi, just another question. I haven't centering the time knots. so I suppose that the intercept refers to the overall mean? maybe it would have been interested to fix the intercept at the time 1 knot..but I'm not sure about how to do this..any advice? thanks a lot
carole

[Clyde Schechter]

If you are referring to the intercept in the model you show in #3, that is the expected value of dep at time 0 with porteur = 0.

I don't see any reason to center the time spline variables. In fact, I think doing that would make the interpretation of the results even more complicated than it already is.

[carole fantini]

thanks Clyde. I perfectly understand that. but time zero do not have any meaning here I suppose because each knots (time 1....time 4) represent the period of time between two measures : baseline to 2 month, 2 month to 6, 6 to ten and 12 to 24. I think I don't need to extrapolate about the meaning of the intercept and just stay focused on the slopes between time and carriers...right?

[Clyde Schechter]

Correct. If time = 0 has no meaning, then the intercept also has no meaning and you should not try to explain it. Focus on times that actually occur.

[carole fantini]

thanks a lot Clyde. that was very helpful for me...

[carole fantini]

Hi, always on the same dataset but i'm troubled actually regarding the missing data. Indeed we have 5 waves of evaluation. I have read that for growth curve with HLM, we need to have 3 data points minimum when collecting 5 waves. it's ok for me. but by reading further, it seems that mixed is not adequate with missing values. is that correct?? and if so, using sem seems better? but know how to draw a growth curve with sem... any help? thanks.

[Clyde Schechter]

Well, if your data are missing at random, then you could used -mixed- with multiple imputation, or you can use -sem- with full information maximum likelihood estimation. Either approach should effectively reduce bias attributable to data that are missing at random. If the data are missing not at random, then neither solution is really effective and there are no truly good approaches unless you have a good statistical model of the missingness generation process.

If you are familiar with -mixed- and not with -sem- then I would say go with -mixed- (and multiple imputation).

[carole fantini]

Hi Clyde, it is what i was thinking. so i don't know how to do multiple imputation with mixed. for Sem, i'm used to do that but not with piecewise growth curve like the model with have talked about (in the previous post). any suggestions??? thanks for your help...carole

[Oded Mcdossi]

You can read in the Stata manual how to do multiple imputation, the manual is clear and has many examples.
As Clyde already note switching to SEM with FIML to reduce bias from missing values, wouldn't be trivial with your model. You can find more information in this link. Richard Williams and his colleagues wrote program xtdpdml (ssc) which is a shell for SEM and make the syntax of this complicated model much easier. I have no experience with this model, let's hope Richard will chime in and will provide some information about how to estimate your model with his program.

[carole fantini]

Thanks for your feedback. multiple implutation with missing values on more than two variables in my models appears, to me, very difficult to do. So i have a look on SEM with FIML but is it Ok with non linear trends and further I have interaction terms? this is the reason why i have used piecewise model with knots....maybe using listwise deletion is easier even not optimal....