Type 3 Tests of Fixed Effects and Mixed Procedure

[Flo Bzh]

Dear all,

In a longitudinal study lasting 6 months, randomised with 2 groups (treated and control), and 3 measurement times point (baseline T0, at 3 month, T3 and at 6 month, T6), I want to know if the change in quality of life (QoL) score is significant in the treated vs control group.

The design of the planed analysis is as follows (I didnt planed it myself; I am not biostatistician but I try to learn with STATA):

Mixed procedure: mean changes from T0 (baseline) in QoL analyzed using a restricted maximum likelihood (REML)-based repeated measures approach including effects of Treatment group, Time (T3 and T6) and Treatment group x Time interaction as well as the covariates of mean centered T0 value of QoL, mean centered T0 value of QoL x Time interaction, Sex, Age, and years of education.
An unstructured covariance matrix will be used to model within-patient errors and the Kenward-Roger approximation will be used to estimate the degrees of freedom of the denominator.

To do so, I first created the variable "Change from_T0" (both for T3 and T6), then I created the new variable QoL_meancentered_T0 and then I reshape my wide data base in long database with the 3 time points measurement for QoL.

I honestly don’t know how to program this analysis in STATA.
I would like first the Type 3 Tests of Fixed Effects. And then the mixed procedure and estimates for treated vs control at T3 and at T6. And the estimates mean change from T0 to T3 and T0 to T6 for each group.
I tried this but I don’t see the “repeated” option…

mixed QoL_score group##time i.sex_t0 age_t0 education_t0 centered_QoL_T0 c. centered_QoL_T0#time, || PARTICIPANTS_ID:, noconstant covariance(unstructured) reml residuals(exchangeable) dfmethod(kroger)

-I feel like I'm forgetting steps
-I don't know how to find the right estimates

If any of you have an idea to perform this analysis I would be gratefull for this precious help

thank you

Flo

[Jeff Wooldridge]

How is QoL measured? Because of randomization, you can start with a linear analysis that includes T0 as a predictor (to possibly help with precision).

I can't say more unless, as discussed in the FAQ, you show Stata output after the commands you typed.

[Flo Bzh]

Thank you Jeff for your reply
the QoL score is a continuous variable between 0 to 30
How do you copy past the Stata output after the commands has been typed ?
I tried something below ... but it doesn't seem very comprehensible.


. mixed QoL_score training_group##time i.sex_t0 age_t0 education_t0 centered_QoL_T0 c.centered_QoL_T0#time if Clinical_group==1, || ID:, noconstant covariance(unstructured) reml residuals(exchangeable) dfmethod(kroger)

Obtaining starting values by EM:

Performing gradient-based optimization:

Iteration 0: log restricted-likelihood = -525.22478
Iteration 1: log restricted-likelihood = -523.85758
Iteration 2: log restricted-likelihood = -523.85564
Iteration 3: log restricted-likelihood = -523.85564

Computing standard errors:

Computing degrees of freedom:

Mixed-effects REML regression Number of obs = 299
Group variable: ID Number of groups = 114

Obs per group:
min = 1
avg = 2.6
max = 3
DF method: Kenward-Roger DF: min = 98.68
avg = 183.30
max = 280.04

F(11, 232.50) = 38.12
Log restricted-likelihood = -523.85564 Prob > F = 0.0000

----------------------------------------------------------------------------------------
QoL_score | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-----------------------+----------------------------------------------------------------
2.training_group | .0383067 .2562507 0.15 0.881 -.4661154 .5427287
|
time |
2 | .5212608 .2577052 2.02 0.044 .0130818 1.02944
3 | .4445715 .2672312 1.66 0.098 -.0822943 .9714373
|
training_group#time |
2 2 | .3563837 .3626092 0.98 0.327 -.3586678 1.071435
2 3 | .3396429 .3705453 0.92 0.360 -.3909721 1.070258
|
1.sex_t0 | -.3532893 .1930922 -1.83 0.070 -.7364414 .0298628
age_t0 | -.0388041 .0116656 -3.33 0.001 -.0619147 -.0156935
education_t0 | .065591 .026659 2.46 0.015 .012748 .118434
centered_QoL_T0 | .9280986 .0609839 15.22 0.000 .8080359 1.048161
|
time#c.centered_QoL_T0 |
2 | -.4038387 .0846852 -4.77 0.000 -.5708087 -.2368687
3 | -.5240611 .0853565 -6.14 0.000 -.6923478 -.3557744
|
_cons | 28.52157 .8331291 34.23 0.000 26.87088 30.17225
----------------------------------------------------------------------------------------

------------------------------------------------------------------------------
Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
-----------------------------+------------------------------------------------
ID: (empty) |
-----------------------------+------------------------------------------------
Residual: Exchangeable |
var(e) | 1.826028 .1540502 1.547736 2.154357
cov(e) | .1861708 .1209263 -.0508405 .4231821
------------------------------------------------------------------------------
LR test vs. linear model: chi2(1) = 2.74 Prob > chi2 = 0.0980

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.

[Leonardo Guizzetti]

In #1, you said you wanted to estimate a mixed model for repeated measures with unstructured covariance matrix. This is a common, short-handed way of meaning unstructured residual covariance structure over (discrete) time measurements.

mixed QoL_score training_group##time i.sex_t0 age_t0 education_t0 centered_QoL_T0 c.centered_QoL_T0#time if Clinical_group==1, || ID:, noconstant covariance(unstructured) reml residuals(exchangeable) dfmethod(kroger)

Note that the covariance() option here is ignored, as Stata tells you when it estimates the model, because you don't have any level-2 random effects (they were suppressed with no constant). This model instead estimates an exchangeable residual covariance structure over measurement times. There's nothing wrong with this set-up per se, but it's a decidedly different structure. Repeated measures over time for QoL tend to diminish in correlation with greater separation in time, so exchangeable structures may be less credible (at least to start with).

Instead, to get the unstructured residual covariance structure over time, change the model to

mixed QoL_score training_group##time i.sex_t0 age_t0 education_t0 centered_QoL_T0 c.centered_QoL_T0#time if Clinical_group==1, || ID:, noconstant covariance(unstructured) reml residuals(un, t(time)) dfmethod(kroger)

where time is as you have coded it, but will be interpreted as a categorical measure of time/visit number.

[Flo Bzh]

Thank you Leonardo
are you sure about this last ligne code ?
because when I type it, it goes on non-stop :

....
numerical derivatives are approximate
flat or discontinuous region encountered
Iteration 24: log restricted-likelihood = 2833.9083 (not concave)
numerical derivatives are approximate
flat or discontinuous region encountered
numerical derivatives are approximate
flat or discontinuous region encountered
numerical derivatives are approximate
flat or discontinuous region encountered
numerical derivatives are approximate
flat or discontinuous region encountered
numerical derivatives are approximate
flat or discontinuous region encountered
Iteration 25: log restricted-likelihood = 2833.9083 (not concave)

[Leonardo Guizzetti]

Thank you Leonardo
are you sure about this last ligne code ?
because when I type it, it goes on non-stop :

....
numerical derivatives are approximate
flat or discontinuous region encountered
Iteration 24: log restricted-likelihood = 2833.9083 (not concave)
numerical derivatives are approximate
flat or discontinuous region encountered
numerical derivatives are approximate
flat or discontinuous region encountered
numerical derivatives are approximate
flat or discontinuous region encountered
numerical derivatives are approximate
flat or discontinuous region encountered
numerical derivatives are approximate
flat or discontinuous region encountered
Iteration 25: log restricted-likelihood = 2833.9083 (not concave)

I am sure that the code modification I showed adheres to your question in #1

An unstructured covariance matrix will be used to model within-patient errors and the Kenward-Roger approximation will be used to estimate the degrees of freedom of the denominator.

However, seeing your results now suggests that this model is poorly identified and will be unable to be fit, so you will need to adapt the residual covariance structure. You may consider exchangeable as you have done, but may also want to consider other structures that allow for diminishing correlation or variance with measurements further apart in time, at least to see if exchangeable is reasonable. These options could be banded, toeplitz or auto-regressive (ar).

[Jeff Wooldridge]

Flo: I assume by baseline you mean the randomization occurred after observing the baseline value, correct?

How come you’re not just using standard regression analysis and computing robust standard errors?

[Flo Bzh]

Hello Jeff,
thank you for helping the poor STATA beginner that I am )
yes baseline values were performed before randomization (we randomized the participants after that they completed all the baseline test procedure).
I am agree with your point 2. There is other statistical method but I would like to stick with the initial statistical plan done by a biostatistician with SAS program. To me, it is an interesting exercice to try how to use STATA to reproduce the analysis that he did with SAS.

[Flo Bzh]

Hello Leonardo Guizzetti
you help me months ago with my question above and your line command works perfectly now. thank you again
After the mixed command I wanna have the effect size but I dont know how to get it from this model
mixed QoL_score training_group##time i.sex_t0 age_t0 education_t0 centered_QoL_T0 c.centered_QoL_T0#time if Clinical_group==1, || ID:, noconstant covariance(unstructured) reml residuals(un, t(time)) dfmethod(kroger) I used first the margins command to have the estimate mean values for each group but I also want the effect size for each group. any idea ? thank you Florent

[Leonardo Guizzetti]

Sorry, I don't use Cohen's d much, and if so, never from the mixed model framework. You also don't specify for which group you want the effect size, and if it's part of the interaction. In any case, following from this paper suggests that you can probably compute an analogue of Cohen's d as follows:

d = (difference in group means) / sqrt(total variance of this difference including random effects)

If you are able to use margins to give you either (a) the group means or (b) the difference in those group means, then this will also include your standard error(s) estimates.

[Flo Bzh]

Yes I have margins results : (2 groups in the design with 2 time point measurement)
other question, do you know how is it possible to have the number of participants for each line (n) ?
is the effect size would be : margin mean / Standard deviation ? (where SD= SE x sqrtn ?


. margins Training_group#time

Predictive margins Number of obs = 190

Expression : Linear prediction, fixed portion, predict()

-------------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
--------------------+----------------------------------------------------------------
Training_group#time |
combined#4 | .2053785 .1053978 1.95 0.051 -.0011974 .4119544
combined#5 | .2110065 .1143886 1.84 0.065 -.013191 .4352041
physical_only#4 | .1413386 .1090014 1.30 0.195 -.0723002 .3549775
physical_only#5 | .1584896 .1190921 1.33 0.183 -.0749265 .3919058
-------------------------------------------------------------------------------------

[Dave Airey]

I think the standard "MMRM" (R side only repeated measures model done with SAS Proc Mixed or Open Pharma MMRM) approach in Pharma in Stata is like:

mixed outcome i.group##i.time /*...other fixed effects like baseline sex etc. */, || personid:, noconstant residuals(unstructured, t(time)) reml dfmethod(kroger, oim)

contrast r.group@time, small