Problem with mixed-effects regression model

[Giulia Nalesso]

Hi everyone,
I have a database of 120 patients with repeated measure (the ild variable identify the patient) of forced vital capacity (fvc) at different timepoints for each patient (pft_date), I also calculated the time from enrolment for each measure in years (fu_time_yrs). I would like to estimate the effect of the variable gerd on fvc decline over time adjusted for the other independent variables (age, bmi, pkyrs)

I'm using statanow 18.5.
This is an example of my dataset

* Example generated by -dataex-. For more info, type help dataex
clear
input byte ild double fvc byte(gerd age pkyrs bmi) int(en_date pft_date) double fu_time_yrs
1  2.1 1 55  5 40 22742 22742                    0
1  1.9 1 55  5 40 22742 22836   .25735797399041754
1  1.7 1 55  5 40 22742 23013    .7419575633127995
1  2.2 1 55  5 40 22742 23108    1.002053388090349
1  2.6 1 55  5 40 22742 23276    1.462012320328542
2  2.7 0 74  0 35 22405 22407 .0054757015742642025
2  2.6 0 74  0 35 22405 22561    .4271047227926078
2  2.6 0 74  0 35 22405 22710    .8350444900752909
2  2.4 0 74  0 35 22405 22997    1.620807665982204
2  2.3 0 74  0 35 22405 23061   1.7960301163586585
3    2 1 68 15 25 22833 22834 .0027378507871321013
3  2.2 1 68 15 25 22833 22927   .25735797399041754
3  2.3 1 68 15 25 22833 23083    .6844626967830253
3  1.2 1 68 15 25 22833 23176    .9390828199863107
4  1.1 1 71 12 28 21967 21977   .02737850787132101
4  1.1 1 71 12 28 21967 22173    .5639972621492129
4  1.1 1 71 12 28 21967 22376   1.1197809719370295
4  1.1 1 71 12 28 21967 22560    1.623545516769336
5  1.1 0 83 30 32 23013 23013                    0
5  1.1 0 83 30 32 23013 23183    .4654346338124572
5  2.8 0 83 30 32 23013 23382   1.0102669404517455
6 3.16 1 78 50 21 22935 22935                    0
6 3.17 1 78 50 21 22935 23045   .30116358658453113
6 3.43 1 78 50 21 22935 23330     1.08145106091718
end
format %tddd-Mon-YY en_date
format %tddd-Mon-YY pft_date

The command that I used is

mixed fvc i.gerd##c.fu_time_yrs bmi age pkyrs || ild:

My result with the example dataset look like this:

Performing EM optimization ...

Performing gradient-based optimization: 
Iteration 0:  Log likelihood = -16.801092  
Iteration 1:  Log likelihood = -16.801092  

Computing standard errors ...

Mixed-effects ML regression                          Number of obs    =     24
Group variable: ild                                  Number of groups =      6
                                                     Obs per group:
                                                                  min =      3
                                                                  avg =    4.0
                                                                  max =      5
                                                     Wald chi2(6)     =  13.36
Log likelihood = -16.801092                          Prob > chi2      = 0.0376

------------------------------------------------------------------------------------
               fvc | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------------+----------------------------------------------------------------
            1.gerd |  -3.487185   1.313645    -2.65   0.008    -6.061882   -.9124885
       fu_time_yrs |   .2345998   .2491368     0.94   0.346    -.2536993     .722899
                   |
gerd#c.fu_time_yrs |
                1  |  -.1533919   .3290906    -0.47   0.641    -.7983977    .4916139
                   |
               bmi |  -.1777524   .0775982    -2.29   0.022    -.3298422   -.0256627
               age |  -.2244968   .0779432    -2.88   0.004    -.3772627   -.0717309
             pkyrs |   .0559209   .0173716     3.22   0.001     .0218732    .0899685
             _cons |   24.70601   8.371693     2.95   0.003     8.297794    41.11423
------------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
ild: Identity                |
                  var(_cons) |   .0889041   .0850662      .0136291    .5799323
-----------------------------+------------------------------------------------
               var(Residual) |   .1814541   .0616237      .0932589     .353056
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 2.25          Prob >= chibar2 = 0.0668

My questions are:
1 Does the command make sense?
2 Am I interpreting the results correctly if I say that gerd impact the fvc change over time with a coefficient of -0.15 and that the impact of gerd on fvc is -3.49?

I apologize if my question is unclear but I'm new on Stata.
Thanks for your response

[Clyde Schechter]

Does the command make sense?

Well, it would with a different data set. But I don't think this kind of model is appropriate with only 6 individuals and an average of four observations for each. For one thing,-mixed- is estimated by maximum likelihood, which is asymptotically correct, but cannot be relied on for small samples. On top of that, with only 6 individuals, you are really not getting a good sample of individual-space, so I think that the ild-level variance component estimation is not really useful.

Furthermore, your covariates are all constant within ild and what you are trying to estimate, namely the rates of change in fvc with fu_time in the presence and absence of gerd, are within-person effects. So you would have a better model using fixed-effects regression instead. Because the covariates are invariant within ild, you can just omit them from the model, and get a consistent estimate for the marginal effects of fu_time on fvc in those two conditions. (In a fixed-effects model, any observed or unobserved attributes of the ild's that do not vary over time are automatically adjusted for.) So I would do this analysis as:

xtset ild
xtreg fvc i.gerd##c.fu_time_yrs, fe
lincom fu_time_yrs // RATE OF FVC CHANGE PER YEAR WITHOUT GERD
lincom fu_time_yrs + 1.gerd#fu_time_yrs // RATE OF FVC CHANGE PER YEAR WITH GERD

Am I interpreting the results correctly if I say that gerd impact the fvc change over time with a coefficient of -0.15...

Since I don't really endorse using that model in this data, I wouldn't put much stock in any interpretation of its results. But just to review how one would interpret this model had the data been more appropriate to it, yes, the coefficient of -0.15 represents the estimated difference between fvc change over time with gerd and without it.

...and that the impact of gerd on fvc is -3.49?

No, that's not correct. The -3.49 is the estimated difference between fvc with and without gerd at follow-up time zero years. It does not apply to fvc measured at other times.

[Joseph Coveney]

I have a database of 120 patients . . . I would like to estimate the effect of the variable gerd on fvc decline over time adjusted for the other independent variables (age, bmi, pkyrs) . . .

The command that I used is
mixed fvc i.gerd##c.fu_time_yrs bmi age pkyrs || ild:

. . . Does the command make sense?

Not if your objective is as you stated, namely, to adjust the differential rate of change for the covariates. As Clyde mentions, that command adjusts only the time-zero intercepts for the two GERD-diagnosis groups.

If you wish to adjust the differential rate of change for those covariates, then you'll need to include each in a three-way interaction involving GERD diagnosis and the follow-up interval for a set of three three-way interaction terms in the regression model. If that's too messy for you, then you'll probably be better off taking Clyde's advice to go with a fixed-effects regression model.

. . . if I say that gerd impact . . . and that the impact of gerd . . . is . . .?

You have an observational study, and so it's a little bold to claim "GERD impacts" and "the impact of GERD is". Adjusting for confounders (age, smoking history, obesity) might help with your particular outcome measure, but endogeneity will still be a concern.

Why didn't you show the model fitted to all 120 patients' data?

[Erik Ruzek]

Just to add on to Clyde's and Joseph's excellent points, if you wanted to go beyond estimating the differential rate of change at the time-zero intercepts for the diagnosis groups, you would need to alter your original model to allow each individual to have a unique rate of change. In the altered model, the gerd#fu_time_yrs interaction fully tests whether the time-zero intercept and rate of change (slope) differs between gerd groups:

mixed fvc i.gerd##c.fu_time_yrs bmi age pkyrs || ild: fu_time_yrs, cov(un) reml

margins gerd, at(fu_time_yrs=(.15(.25)1.10))
marginsplot 

* Compare the above marginsplot with the one from your original model
mixed fvc i.gerd##c.fu_time_yrs bmi age pkyrs || ild: ,  reml

margins gerd, at(fu_time_yrs=(.15(.25)1.10))
marginsplot

Further, as Joseph said, to adjust for the covariates on the differential rates of change by gerd, you need to add interactions between each of the covariates and c.fu_time_yrs#i.gerd. This model would be overkill for the data example you provided, but will probably run with your larger dataset.

[Giulia Nalesso]

Thank you for all your answers, they were very helpful and they clarified most of my doubts. I know the -mixed- command could't fit for a dataset like the one I presented in the example, I have a larger dataset but I was not sure if I could publish the data.