[carole fantini]

Hi Clyde, I hope this will be my last question about this topic. I have simplified my model by using a group variable with four level leading to a "simple" interaction between one categorical IV and one continuous (time splines, marginal). I'm always struggling with slopes comparisons.after carefully reading the book of Mitchell (excellent !!!!!!), I have used the contrast command - contrast r2b4.group#c.(time1 time2 time3 time4 time5), effects- in order to compare the slopes of group 2 (Non carrier HBOC) against group4 (carrier HBOC). I have a significant effect for time 2 which is attended) I just don't have any trust in my interpretation which is : the slope of carrier vs non carrier is significant for the knot 2. BUT, as knots have been defined with the marginal command, I m not sure about this interpretation at all!!! does it means that each group differ on the time 2 knot respective of their own previous knot (time1)? thanks

[Clyde Schechter]

I think it would be best if you showed the code and output so I can be sure I understand what you've described. Off hand, it sounds like you're on the right track, but I'd prefer the increased confidence that comes with knowing that what's in my imagination is the same as what you're working with.

[carole fantini]

hi Clyde, here is the information I obtained first from the mixed and then to the contrast...

 mixed dep i.group##c.(time1 time2 time3 time4 time5) CRc ESc ||ID: time1 time2 time3 time4 time5 CRc ESc, 

Performing EM optimization: 

Performing gradient-based optimization: 

Iteration 0:   log likelihood = -2622.0521  
Iteration 1:   log likelihood = -2600.5559  (not concave)
Iteration 2:   log likelihood = -2600.1394  (not concave)
Iteration 3:   log likelihood = -2599.7892  
Iteration 4:   log likelihood = -2598.1474  
Iteration 5:   log likelihood =  -2596.834  
Iteration 6:   log likelihood = -2596.7866  
Iteration 7:   log likelihood = -2596.7822  
Iteration 8:   log likelihood =  -2596.782  

Computing standard errors:

Mixed-effects ML regression                     Number of obs     =        737
Group variable: ID                              Number of groups  =        156

                                                Obs per group:
                                                              min =          3
                                                              avg =        4.7
                                                              max =          6

                                                Wald chi2(25)     =      84.56
Log likelihood =  -2596.782                     Prob > chi2       =     0.0000

-------------------------------------------------------------------------------
          dep |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------+----------------------------------------------------------------
        group |
      0 HBOC  |   5.835396    3.76847     1.55   0.122     -1.55067    13.22146
     1 HNPCC  |   3.324368   4.084699     0.81   0.416    -4.681495    11.33023
      1 HBOC  |   6.291214   4.299187     1.46   0.143    -2.135038    14.71747
              |
        time1 |   .8381888   1.720544     0.49   0.626    -2.534016    4.210394
        time2 |   -1.34031   3.025853    -0.44   0.658    -7.270872    4.590252
        time3 |   4.281338   3.128444     1.37   0.171    -1.850298    10.41297
        time4 |  -.9496321   3.450392    -0.28   0.783    -7.712277    5.813012
        time5 |  -6.941403   3.826365    -1.81   0.070    -14.44094    .5581354
              |
group#c.time1 |
      0 HBOC  |  -1.576868   2.106796    -0.75   0.454    -5.706112    2.552376
     1 HNPCC  |  -1.456218   2.292118    -0.64   0.525    -5.948687     3.03625
      1 HBOC  |  -2.568753   2.448139    -1.05   0.294    -7.367018    2.229511
              |
group#c.time2 |
      0 HBOC  |  -.7178708   3.691509    -0.19   0.846    -7.953095    6.517354
     1 HNPCC  |   4.337357   4.005089     1.08   0.279    -3.512474    12.18719
      1 HBOC  |   7.305894   4.314151     1.69   0.090    -1.149686    15.76147
              |
group#c.time3 |
      0 HBOC  |   2.790232   3.840908     0.73   0.468    -4.737811    10.31827
     1 HNPCC  |  -8.165621   4.145545    -1.97   0.049    -16.29074   -.0405032
      1 HBOC  |  -16.35531   4.482354    -3.65   0.000    -25.14056   -7.570054
              |
group#c.time4 |
      0 HBOC  |  -1.178429   4.225278    -0.28   0.780    -9.459822    7.102963
     1 HNPCC  |   7.565985   4.661469     1.62   0.105    -1.570326     16.7023
      1 HBOC  |    17.0437   5.035661     3.38   0.001     7.173988    26.91342
              |
group#c.time5 |
      0 HBOC  |   .4423934   4.588062     0.10   0.923    -8.550044     9.43483
     1 HNPCC  |  -2.770583   5.305468    -0.52   0.602    -13.16911    7.627944
      1 HBOC  |   -3.96626   5.315317    -0.75   0.456    -14.38409    6.451569
              |
          CRc |    .128873   .0961463     1.34   0.180    -.0595702    .3173162
          ESc |   .2775783   .1383555     2.01   0.045     .0064066    .5487501
        _cons |   8.538832   3.060521     2.79   0.005     2.540321    14.53734
-------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
ID: Independent              |
                  var(time1) |   1.97e-11          .             .           .
                  var(time2) |   8.54e-06          .             .           .
                  var(time3) |   1.18e-08          .             .           .
                  var(time4) |   7.327946          .             .           .
                  var(time5) |   1.16e-08          .             .           .
                    var(CRc) |   1.70e-14          .             .           .
                    var(ESc) |    .249485          .             .           .
                  var(_cons) |    53.5015          .             .           .
-----------------------------+------------------------------------------------
               var(Residual) |   41.01446          .             .           .
------------------------------------------------------------------------------
LR test vs. linear model: chi2(8) = 327.51                Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.

. contrast r2b4.group#c.(time1 time2 time3 time4 time5), effects

Contrasts of marginal linear predictions

Margins      : asbalanced

-------------------------------------------------
              |         df        chi2     P>chi2
--------------+----------------------------------
dep           |
group#c.time1 |          1        0.22     0.6405
              |
group#c.time2 |          1        4.62     0.0316
              |
group#c.time3 |          1       24.01     0.0000
              |
group#c.time4 |          1       17.12     0.0000
              |
group#c.time5 |          1        0.97     0.3245
-------------------------------------------------

-------------------------------------------------------------------------------------
                    |   Contrast   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------------+----------------------------------------------------------------
dep                 |
      group#c.time1 |
(0 HBOC vs 1 HBOC)  |   .9918852   2.123992     0.47   0.641    -3.171064    5.154834
                    |
      group#c.time2 |
(0 HBOC vs 1 HBOC)  |  -8.023765   3.731973    -2.15   0.032     -15.3383   -.7092329
                    |
      group#c.time3 |
(0 HBOC vs 1 HBOC)  |   19.14554   3.907408     4.90   0.000     11.48716    26.80392
                    |
      group#c.time4 |
(0 HBOC vs 1 HBOC)  |  -18.22213   4.404216    -4.14   0.000    -26.85424   -9.590026
                    |
      group#c.time5 |
(0 HBOC vs 1 HBOC)  |   4.408653   4.474398     0.99   0.324    -4.361005    13.17831

[Clyde Schechter]

The time* variables being a marginal spline, I would interpret this as: in the time periods from time 2 to time 3, time 3 to time 4 and time 4 to time 5, the rate of change per unit of time changes by a statistically significantly different amount in the two HBOC groups. That's a mouthful. At the risk of some unclarity, one could say this a bit more succinctly as the acceleration of dep is different in the two groups from time 2 through time 5. (Yes, that's right, it's the difference in second derivatives!)

To be honest, it's hard to imagine how you can explain this in words to any but the most non-technical audiences. I think that in your (oral and written) presentations of this work you will need to rely heavily on graphs.

[carole fantini]

yes it's very hard to translate such a thing with a meaningful way for psychologist ! i'm manly searching after a difference between carrier and non carriers on depression scores, on some specific splines (which is difficult with marginal mksplines). say that the rate of change between the two group is different is not helpful except if we find a solution to show that depressions scores are higher (with piecewise). i'm wondering if it s not more relevant to use mkspline to define the knots without "marginal", just to be more simple and test really if one group increase (carriers HBOC) while the other one decrease (non carrier) after the announcement of the genetic test result. Further, i'm troubled by the fact that the genetic test result which is a non time varying variable, is also considered for the knot one (where people doesn't knew their result). But i remember that you suggest to include all main effect and interactions terms.

I show you the graph of the depression scores for my populations (the period of time are not equally distant between each point : it was : base, 15 days, 2 montth, 6 12 and 24).

[Clyde Schechter]

Yes, it might be easier to explain the results using the non-marginal spline.

As for the genetic test being considered at the beginning, this is not a problem. Whatever difference you find between them represents some difference that is either just noise or perhaps some effect of the genotype itself that acts in the absence of knowing you have it.

[carole fantini]

thanks Clyde. very helpful for me.

[Vicky Bajpai]

Hi Clyde,
I am also running the similar analysis. I am trying to do piecewise multilevel modelling for a continuous outcome (say pain score) and trying to find association with BMI over time (baseline, 12, and 24 months) by adjusting other confounders. I exactly followed as you suggested in this conversation. My codes are given below [confounders not shown here];

mkspline time = time, cubic nknots(3)
mixed a_1 i.bmi_grp time1 time2 i.bmi_grp#c.(time1 time2) || surveyid: time1 time2, reml
forvalues t = 1/3 {
    local at`t'
    forvalues i = 1/2 {
        summ time`i' if time == `t', meanonly
        local at`t' `at`t'' time`i' = `r(mean)'
    }
    display "at(`at`t'')"
}
tempfile results
margins bmi_grp, at(`at1') at(`at2') at(`at3') saving(`results')

But the last line [margins.........(results)] is not running and prompting the error 'option at() not allowed'. Please help me to sort out this issue.

[Clyde Schechter]

I don't think it's possible to troubleshoot this without an example of your data (use -dataex-), and all of the output you got from Stata from the code you show in #38.

Added: You might find it helpful to look at a different approach to using -margins- with splines at https://www.statalist.org/forums/for...ne-of-best-fit.

[Vicky Bajpai]

Hi Clyde,
Thank you for your reply. Please see below

example of data, and all of the output got from Stata;

mkspline tc = time, cubic nknots(3)

input int surveyid byte a_1 float(bmi_grp tc1 tc2)
1  7 30  0  0
1  8  30 12  3
1  7  25 24 18
2  9 30  0  0
2  4 30 12  3
2  4 30 24 18
3 10 25  0  0
3  6 25 12  3
3  2 30 24 18
4  . 25  0  0
end

mixed a_1 i.bmi_grp tc1 tc2 i.bmi_grp#c.(tc1 tc2) || surveyid: tc1 tc2, reml

. mixed a_1 i.bmi_grp tc1 tc2 i.bmi_grp#c.(tc1 tc2) || surveyid: tc1 tc2, reml

Performing EM optimization:

Performing gradient-based optimization:

Iteration 0:   log restricted-likelihood = -3260.1463 
Iteration 1:   log restricted-likelihood = -3243.8207 
Iteration 2:   log restricted-likelihood = -3237.7075 
Iteration 3:   log restricted-likelihood =  -3237.187 
Iteration 4:   log restricted-likelihood = -3237.1855 
Iteration 5:   log restricted-likelihood = -3237.1855 

Computing standard errors:

Mixed-effects REML regression                   Number of obs     =      1,472
Group variable: surveyid                        Number of groups  =        621

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

                                                Wald chi2(8)      =    2001.61
Log restricted-likelihood = -3237.1855          Prob > chi2       =     0.0000

------------------------------------------------------------------------------
         a_1 |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     bmi_grp |
    25-29.9  |   .2504804     .18312     1.37   0.171    -.1084282     .609389
        >30  |   .4815063   .2059685     2.34   0.019     .0778155    .8851971
             |
         tc1 |   -.540444   .0234114   -23.08   0.000    -.5863294   -.4945585
         tc2 |   .4386729   .0300224    14.61   0.000     .3798301    .4975156
             |
     bmi_grp#|
       c.tc1 |
    25-29.9  |   .0009325   .0316658     0.03   0.977    -.0611313    .0629963
        >30  |   .0681842   .0345966     1.97   0.049     .0003761    .1359924
             |
     bmi_grp#|
       c.tc2 |
    25-29.9  |   .0073628    .040896     0.18   0.857    -.0727919    .0875175
        >30  |  -.0765332   .0442284    -1.73   0.084    -.1632192    .0101528
             |
       _cons |   7.578771   .1355195    55.92   0.000     7.313158    7.844385
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
surveyid: Independent        |
                    var(tc1) |   .0069362   .0009787      .0052604    .0091458
                    var(tc2) |   5.28e-22          .             .           .
                  var(_cons) |   .9599879   .1815316      .6626823    1.390677
-----------------------------+------------------------------------------------
               var(Residual) |   2.915068   .1692614      2.601502    3.266428
------------------------------------------------------------------------------
LR test vs. linear model: chi2(3) = 174.34                Prob > chi2 = 0.0000

mixed a_1 i.bmi_grp tc1 tc2 i.bmi_grp#c.(tc1 tc2) || surveyid: tc1 tc2, reml
forvalues t = 1/3 {
    local at`t'
    forvalues i = 1/2 {
        summ tc`i' if time == `t', meanonly
        local at`t' `at`t'' time`i' = `r(mean)'
    }
    display "at(`at`t'')"
}

. margins bmi_grp, at(`at1') at(`at2') at(`at3')
option at() not allowed

I also tried your other suggestion but

margins, `ats'

is not producing results at different time points. please see below;

. margins, `ats' nose saving(`margins_results')

Predictive margins                              Number of obs     =      1,472

Expression   : Linear prediction, fixed portion, predict()

------------------------------------------------------------------------------
             |     Margin
-------------+----------------------------------------------------------------
       _cons |   4.901581
------------------------------------------------------------------------------

[Clyde Schechter]

The example data you show cannot be the data you are running the code you show on, because it does not contain the variable time. So the code you show would have stopped long before ever getting to -margins- with "time not found" errors.

Looking into it farther, the first spline variable in Stata's cubic splines is always equal to the original variable. So by running -gen time = tc1- I can presumably recreate the data you actually started with. Doing that, the problem becomes pretty clear. The outermost loop is -forvalues t = 1/3 {-, and within the loop we see -summ tc`i' if time == `t', meanonly-. But within the data example you sent, the variable time never takes on the values 1, 2, or 3. So the -summ- command will return nothing for the mean. The resulting locals at1, at2, and at3 are ill-formed. For example, local at1 comes out as at(time1 = timie2 =), which is clearly invalid syntax and is the cause of the error message you are getting.

You need to revise that -forvalues- loop to iterate over the actual values of your time variable. That said, do look at #2 in https://www.statalist.org/forums/for...ne-of-best-fit. That approach does not require you to know the values of your time variable in advance and enumerate them in the loop: it finds them for you and incorporates them as needed. I think you will find that approach more convenient.

[Vicky Bajpai]

Hi Clyde,
Thank you for your wonderful suggestion. I modified the time variable accordingly and rerun the analysis. Now loop seems working but 'margin' have some problem as it did not generate means at each time point. It's just generating overall means. Please see all results below and suggest how should I get that.

I'd also like to compare BMI categories at 12 and 24 months from baseline. How could generate p-value?

. mixed a_1 i.bmi_grp#c.(tc1 tc2) age i.gender i.b_1 i.g7a i.g8a || surveyid: tc1 tc2, reml

Performing EM optimization:

Performing gradient-based optimization:

Iteration 0:   log restricted-likelihood = -3177.1571 
Iteration 1:   log restricted-likelihood = -3156.3513 
Iteration 2:   log restricted-likelihood = -3155.0736 
Iteration 3:   log restricted-likelihood =  -3155.004 
Iteration 4:   log restricted-likelihood = -3155.0039 

Computing standard errors:

Mixed-effects REML regression                   Number of obs     =      1,445
Group variable: surveyid                        Number of groups  =        613

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

                                                Wald chi2(14)     =    2113.95
Log restricted-likelihood = -3155.0039          Prob > chi2       =     0.0000

------------------------------------------------------------------------------
         a_1 |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     bmi_grp#|
       c.tc1 |
        <25  |  -6.372076    .188325   -33.84   0.000    -6.741187   -6.002966
    25-29.9  |  -6.156112   .1789593   -34.40   0.000    -6.506866   -5.805359
        >30  |  -5.833497   .1865246   -31.27   0.000    -6.199079   -5.467916
             |
     bmi_grp#|
       c.tc2 |
        <25  |   5.255249   .2802262    18.75   0.000     4.706016    5.804483
    25-29.9  |   5.085107   .2606108    19.51   0.000     4.574319    5.595894
        >30  |   4.613157   .2916723    15.82   0.000      4.04149    5.184825
             |
         age |   .0034696   .0074851     0.46   0.643    -.0112009    .0181401
             |
      gender |
       Male  |  -.2851943   .1484266    -1.92   0.055    -.5761051    .0057166
             |
         b_1 |
        Yes  |   .4140827   .1792336     2.31   0.021     .0627913     .765374
             |
         g7a |
 Previously  |  -.0679927   .1421341    -0.48   0.632    -.3465704     .210585
  Currently  |   .4466982   .3029811     1.47   0.140    -.1471339     1.04053
             |
         g8a |
1-2 times..  |   .1232898   .1820733     0.68   0.498    -.2335673     .480147
Special o..  |   .8667378   .1937989     4.47   0.000      .486899    1.246577
      Never  |   .5257494   .2229372     2.36   0.018     .0888007    .9626982
             |
       _cons |   13.05482     .60982    21.41   0.000      11.8596    14.25005
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
surveyid: Independent        |
                    var(tc1) |   .4448684          .             .           .
                    var(tc2) |   3.13e-13          .             .           .
                  var(_cons) |   9.58e-12          .             .           .
-----------------------------+------------------------------------------------
               var(Residual) |   3.278111          .             .           .
------------------------------------------------------------------------------
LR test vs. linear model: chi2(3) = 151.97                Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.

. forvalues t = 1/3 {
  2.     local at`t'
  3.     forvalues i = 1/2 {
  4.         summ tc`i' if time1 == `t', meanonly
  5.         local at`t' `at`t'' time1`i' = `r(mean)'
  6.     }
  7.     display "at(`at`t'')"
  8. }
at(time11 = 1 time12 = 0)
at(time11 = 2 time12 = .25)
at(time11 = 3 time12 = 1.5)

It still does not allow to run the following command;

. margins bmi_grp, at(`at1') at(`at2') at(`at3') saving(`results')
option at() not allowed

If I try at ('at1') and ('at2') separately then it produced the same results.


.

 margins bmi_grp, at(`at1')

Predictive margins                              Number of obs     =      1,445

Expression   : Linear prediction, fixed portion, predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     bmi_grp |
        <25  |   4.587498    .121318    37.81   0.000     4.349719    4.825276
    25-29.9  |   4.905767    .106845    45.91   0.000     4.696355    5.115179
        >30  |    5.27571   .1327491    39.74   0.000     5.015526    5.535893
------------------------------------------------------------------------------

. margins bmi_grp, at(`at2')

Predictive margins                              Number of obs     =      1,445

Expression   : Linear prediction, fixed portion, predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     bmi_grp |
        <25  |   4.587498    .121318    37.81   0.000     4.349719    4.825276
    25-29.9  |   4.905767    .106845    45.91   0.000     4.696355    5.115179
        >30  |    5.27571   .1327491    39.74   0.000     5.015526    5.535893
------------------------------------------------------------------------------

[Clyde Schechter]

. forvalues t = 1/3 { 2. local at`t' 3. forvalues i = 1/2 { 4. summ tc`i' if time1 == `t', meanonly 5. local at`t' `at`t'' time1`i' = `r(mean)' 6. } 7. display "at(`at`t'')" 8. } at(time11 = 1 time12 = 0) at(time11 = 2 time12 = .25) at(time11 = 3 time12 = 1.5)

makes no sense here, as your regression model has no variables time11 and time12. So it is not surprising that -margins- will not run it. What does surprise me is that it does run -at(`at1')- separately. In any case, you cannot simply copy/paste the code I have shown at the link. You have to adapt it to match your model.

[Vicky Bajpai]

Thanks Clyde for your guidance. I managed to customise codes to my scenario.