Construct multiple linear mixed effect model

[Hang Tran]

Hello Statalister,

I'm struggling with constructing mixed effect model for my thesis
I'm using Stata 18 BE software, and I first time run mixed effect model by applying the book " Multilevel modeling in plain language"
I don't know how to construct mixed effect model as table 4 in previous publication

I constructed model as:
[Model 1: mixed Physical_P|| ID:
Model 2: mixed Physical_P totaldate || ID:
Model 3: mixed Physical_P toddler yc child ado totaldate immu education
Model 4:mixed Physical_P toddler yc child ado totaldate immu education]
I have 4 aged groups with total assessment as following: toddler (aged 2-4): 104 ; yc- young child (aged 5-7): 120; child (aged 8-12): 134; ado: adolescent (aged 13-18): 54
my adviser told me construct model for toddler group separately with 3 other groups (young child, child, adolescent).
I would like to construct model for 3 aged group as table 4 but I couldn't figure out how to do ..


I greatly appreciated your time and your advance help...

[Hang Tran]

[Erik Ruzek]

Welcome to the forum! I think we need a little more information before we can help you. Here are a couple of questions I have.

1. What is the hierarchical nesting of the data. Is it that children are measured repeatedly (longitudinal)? Or is it that children are nested within larger groups such as schools, regions, or something else?
2. Can you please use the dataex command to provide us with a sample of your data - the outcome, key predictors and the grouping variable that leads you to believe you need a multilevel model?

[Hang Tran]

I greatly appreciate your quick respones,
the hierarchical nesting of the data is childrean are measured repeatedly but not as classic longitudinal, my study is cross-sectional design, some of them with repeated measurements in 2 years
I did the dataex command

dataex ID CAge AgeGroup Physical_P Physical education immu employment

The outcome is Quality of life (QoL) fuctionings, i.e, Physical functioning
Key predictors here is times after diagnosis (totaldate: continuous variable), proxy's education(education: binary variable), immuno type- acute lymphoblastic leukemia B/ ALL T, (immu: binary varialbe), proxy's employment (employment: binary variable)

my hypothesis is.., different age group have different in dynamic changes in QoL ( based on figure about dynamic change in QoL domains which make by applying kernel smoothing mean)

[Erik Ruzek]

dataex produces output syntax that you can copy and post to the forum. When you do so, paste it between two code blocks (click on the # in the post control options.

[Hang Tran]

	
	
* Example generated by -dataex-. For more info, type help dataex
clear
input float(ID totaldate) byte CAge double Physical_P float(education immu employment)
 1  402  3             93.75 0 1 1
 1  895  4            84.375 0 1 1
 1  920  4             81.25 0 1 1
 2  105  3             56.25 . . .
 2  598  4               100 1 1 0
 4  281  7                75 . 1 .
 4  774  8               100 0 1 1
 5  181  5                75 0 1 1
 5  708  6            65.625 0 1 0
 7 1309  5             43.75 . 1 .
 8 2400  7               100 0 1 0
 9  321  3             81.25 0 1 0
10    8 12                50 0 1 0
11   30  4             93.75 0 1 1
12  198 14             93.75 0 2 1
12  710 15            15.625 0 2 1
17  457  8            59.375 0 . 1
17  510  9             81.25 0 1 1
20  291  3             93.75 1 1 0
21  584  4             56.25 1 1 1
21 1128  5              12.5 1 1 .
23 1318  6                75 1 . 1
23 1350  7             81.25 1 . 1
25   26  8            71.875 0 . 1
26   86  4            71.875 0 1 1
27   23 13               100 0 1 0
29  441  6               100 1 1 1
30    5 13            59.375 0 . 1
30  482 14             68.75 0 1 1
33  614  3            71.875 0 1 1
33 1118  4            78.125 1 1 1
35  289  3             81.25 0 1 1
35  766  4            71.875 0 1 1
35  969  5            65.625 0 1 1
36   63  3             3.125 0 1 1
36  478  4             81.25 0 1 1
38  149  5              62.5 0 . 1
38  605  6               100 0 1 1
40  113  4             56.25 0 . 1
40  548  5            71.875 0 . 1
42  218 11             68.75 0 2 1
42  513 12               100 0 2 1
42  722 12               100 0 2 1
42  744 12              87.5 0 2 1
42  773 12             68.75 0 2 1
43   63  6              62.5 0 . 1
43  456  7                50 0 1 1
44  266  3              87.5 1 1 1
44  751  4                75 0 1 1
44  841  4             93.75 . 1 .
45  380  6            84.375 0 1 1
46  933  9            78.125 0 1 1
47  240 10             56.25 0 1 0
48  457 13              87.5 0 2 0
49  291 13             56.25 0 1 1
51  750  9             43.75 1 1 1
52  142  4               100 0 1 1
53  355  4               100 0 1 1
54  111  9             81.25 0 1 0
55  173  2                75 0 1 1
56   40  7              87.5 0 1 1
57  337  3 92.85714285714286 0 1 1
58   70  4               100 0 1 1
58  353  4             93.75 0 1 0
60   98  3             93.75 0 1 1
61  110  8 53.57142857142857 . 1 .
62   54  4            65.625 0 1 1
62  548  5            46.875 0 . 1
64  216  3                50 0 1 1
64  709  4             56.25 0 1 1
64  773  4            59.375 1 1 1
66   83 15             68.75 0 1 1
66  590  6            96.875 0 1 1
68  375  3             81.25 1 1 1
68  889  4               100 1 1 1
68  918  4               100 1 1 1
69  484 10            40.625 0 1 1
69  957 11            78.125 0 1 1
70  114  6             93.75 0 1 1
70  144  7              87.5 0 . 0
71  142  5             93.75 0 1 1
72   44  4            78.125 1 1 1
73  252 11            84.375 0 1 1
74   18 15             56.25 1 1 1
75   20 13              87.5 0 2 1
77  125  2              87.5 0 1 1
83   51  6            53.125 0 1 1
83  309  7              62.5 0 1 1
83  689  8             81.25 . 2 1
84  890  6               100 1 1 1
85  655 12             68.75 1 2 1
86  514  9            90.625 0 . 1
86  523 10             81.25 0 . 1
87  243  6             81.25 0 1 1
89 1175  7             81.25 0 1 1
90 1101 13               100 0 1 1
96  664  3            71.875 1 . 1
96  844  4            71.875 1 2 1
96  938  4             56.25 1 1 1
97  534  4            78.125 0 1 1
end

[Hang Tran]

Thank you, Erik.
Above is primary data from proxy data set (I also have child data set). From child age (CAge) variable, I seperated to four groups for proxy dataset as toddler (aged 2-4), young child (aged 5-7), child (aged 8-12), and adolescent (aged 13-18)
I constructed the mixed effect model again for 3 groups (young child, child, adolescent) that help me easy to compare with result from child dataset which also separate same 3 aged groups

 quietly mixed Physical_P|| ID:
estat icc
estadd scalar icc = r(icc2)
estimates store nullphy
quietly mixed Physical_P totaldate || ID:
estat icc
estadd scalar icc = r(icc2)
estimates store mod2
quietly mixed Physical_P youngc child ado totaldate immu education employment|| ID: 
estat icc
estadd scalar icc = r(icc2)
estimates store mod3
quietly mixed Physical_P youngc child ado totaldate immu education employment|| ID: totaldate, covariance(unstructured) nolog
estadd scalar icc = r(icc2)
estimates store mod4
esttab nullphy mod2 mod3 mod4, label se aic bic scalars(icc ll df_m) transform(ln*: exp(@)^2 exp(@)^2) eqlabels("" "var(Constant)" "var(Residual)" )

and here is the result

------------------------------------------------------------------------------------
(1) (2) (3) (4)
Physical_P Physical_P Physical_P Physical_P
------------------------------------------------------------------------------------

totaldate 0.0142*** 0.0142*** 0.0142***
(0.00277) (0.00277) (0.00277)

Constant 70.40*** 63.90*** 63.90*** 63.90***
(1.308) (1.775) (1.775) (1.775)
------------------------------------------------------------------------------------
var(Constant)
Constant 319.9*** 270.8*** 270.8*** 270.8***
(23.68) (22.34) (22.34) (22.34)
------------------------------------------------------------------------------------
var(Residual)
Constant 263.8*** 268.6*** 268.6*** 268.6***
(14.80) (15.21) (15.21) (15.21)
------------------------------------------------------------------------------------
Observations 460 460 460 460
AIC 4178.6 4155.5 4155.5 4155.5
BIC 4190.9 4172.0 4172.0 4172.0
icc 0.548 0.502 0.502 0.502
ll -2086.3 -2073.8 -2073.8 -2073.8
df_m 0 1 1 1
------------------------------------------------------------------------------------
Standard errors in parentheses
* p<0.05, ** p<0.01, *** p<0.001

but seems I'm doing wrong.

Thank you for your help

[Erik Ruzek]

Thank you for sharing that. Can you say more about why you think you are doing something wrong? What is that makes you think this?

The analysis should be start with a statement of your guiding research question and any hypotheses you have.

My guess is that you want to know whether the association between time from diagnosis (totaldate) and the outcome varies by child age category. If so, then this is a start:

* First, recode the CAge variable
recode CAge (2/4 = 0 "toddler") (5/7 = 1 "young child") (8/12 = 2 "child") ///
            (13/18 = 3 "adolescent"), gen(agecat)
tab agecat

* Next, look at the data, for example:
twoway (scatter Physical_P totaldate) (lowess Physical_P totaldate) if agecat==0
twoway (scatter Physical_P totaldate) (lowess Physical_P totaldate) if agecat==1
twoway (scatter Physical_P totaldate) (lowess Physical_P totaldate) if agecat==2
twoway (scatter Physical_P totaldate) (lowess Physical_P totaldate) if agecat==3

* Then run some candidate models that address your research question and test your hypotheses
*Note, you have missing data in the predictors
mixed Physical_P i.agecat totaldate immu education || ID:
eststo m2


mixed Physical_P i.agecat##c.totaldate immu education || ID:
eststo m3
lrtest m3 m2, stats     // no real improvement in model fit by including interactions

margins agecat, at(totaldate = (0(100)1000))
marginsplot

[Hang Tran]

Hi Erik,
Thank you very much, your question is so helpful for me to identify clearly my problem
At the beginning, I want to know whether the association between time after diagnosis and physical functioning ? and I also want to identify determinants related to QoL (i.e, Physical functiong)
I ran as your suggestions and obtained results

Likelihood-ratio test
Assumption: m2 nested within m3

LR chi2(3) = 7.80
Prob > chi2 = 0.0503

Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
Model | N ll(null) ll(model) df AIC BIC
-------------+---------------------------------------------------------------
m2 | 363 . -1641.232 9 3300.465 3335.515
m3 | 363 . -1637.331 12 3298.663 3345.395
-----------------------------------------------------------------------------

Should I pick a model with smaller AIC or smaller BIC?
The result of margin predictive as

margins AgeGroup, at(totaldate = (0(100)1000))

Predictive margins Number of obs = 363

Expression: Linear prediction, fixed portion, predict()
1._at: totaldate = 0
2._at: totaldate = 100
3._at: totaldate = 200
4._at: totaldate = 300
5._at: totaldate = 400
6._at: totaldate = 500
7._at: totaldate = 600
8._at: totaldate = 700
9._at: totaldate = 800
10._at: totaldate = 900
11._at: totaldate = 1000

------------------------------------------------------------------------------
| Delta-method
| Margin std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_at#AgeGroup |
1 0 | 68.78552 3.396808 20.25 0.000 62.1279 75.44315
1 1 | 67.39163 3.837148 17.56 0.000 59.87096 74.9123
1 2 | 54.03324 3.761237 14.37 0.000 46.66135 61.40513
1 3 | 63.82427 5.631117 11.33 0.000 52.78749 74.86106
2 0 | 70.07559 2.964575 23.64 0.000 64.26513 75.88605
2 1 | 68.3336 3.467213 19.71 0.000 61.53799 75.12922
2 2 | 56.63208 3.362134 16.84 0.000 50.04242 63.22175
2 3 | 63.41164 4.860423 13.05 0.000 53.88538 72.93789
3 0 | 71.36566 2.614599 27.30 0.000 66.24114 76.49018
3 1 | 69.27558 3.140113 22.06 0.000 63.12107 75.43008
3 2 | 59.23093 3.013341 19.66 0.000 53.32489 65.13697
3 3 | 62.999 4.260472 14.79 0.000 54.64863 71.34937
4 0 | 72.65573 2.383394 30.48 0.000 67.98436 77.3271
4 1 | 70.21755 2.870527 24.46 0.000 64.59142 75.84368
4 2 | 61.82978 2.73418 22.61 0.000 56.47088 67.18867
4 3 | 62.58636 3.910654 16.00 0.000 54.92162 70.2511
5 0 | 73.9458 2.306952 32.05 0.000 69.42426 78.46735
5 1 | 71.15953 2.675896 26.59 0.000 65.91486 76.40419
5 2 | 64.42862 2.547645 25.29 0.000 59.43533 69.42192
5 3 | 62.17372 3.879238 16.03 0.000 54.57056 69.77689
6 0 | 75.23587 2.400106 31.35 0.000 70.53175 79.93999
6 1 | 72.1015 2.573284 28.02 0.000 67.05796 77.14504
6 2 | 67.02747 2.474771 27.08 0.000 62.17701 71.87793
6 3 | 61.76109 4.173421 14.80 0.000 53.58133 69.94084
7 0 | 76.52594 2.644997 28.93 0.000 71.34184 81.71004
7 1 | 73.04347 2.573721 28.38 0.000 67.99907 78.08787
7 2 | 69.62632 2.525416 27.57 0.000 64.67659 74.57604
7 3 | 61.34845 4.732873 12.96 0.000 52.07219 70.62471
8 0 | 77.81601 3.004751 25.90 0.000 71.92681 83.70521
8 1 | 73.98545 2.677155 27.64 0.000 68.73832 79.23258
8 2 | 72.22516 2.692619 26.82 0.000 66.94773 77.5026
8 3 | 60.93581 5.476897 11.13 0.000 50.20129 71.67033
9 0 | 79.10608 3.443554 22.97 0.000 72.35684 85.85532
9 1 | 74.92742 2.872483 26.08 0.000 69.29746 80.55739
9 2 | 74.82401 2.956671 25.31 0.000 69.02904 80.61898
9 3 | 60.52317 6.340852 9.54 0.000 48.09533 72.95102
10 0 | 80.39615 3.935051 20.43 0.000 72.68359 88.10871
10 1 | 75.8694 3.142616 24.14 0.000 69.70998 82.02881
10 2 | 77.42286 3.294367 23.50 0.000 70.96602 83.8797
10 3 | 60.11054 7.282177 8.25 0.000 45.83773 74.38334
11 0 | 81.68622 4.461861 18.31 0.000 72.94113 90.4313
11 1 | 76.81137 3.470128 22.14 0.000 70.01005 83.6127
11 2 | 80.0217 3.685518 21.71 0.000 72.79822 87.24519
11 3 | 59.6979 8.274507 7.21 0.000 43.48016 75.91563
------------------------------------------------------------------------------

Could you please give me more suggestion in case my hypothesis is : is it better physical function (or emotion, social...) by different timeline in different age group
My previous model result couldn't answer for my hypothesis because it did not show the result in different timeduration as below line graph showed.
As you can see there are 3 line, 5-7 yrs( blue), 8-12 yrs(pink), 13-18 yrs(brown), and my hypothesis come from that line graph. the group aged 8-12 seem linear increase by time, not two other groups.
To Emotion, all 3 groups fluctuated, then I confused how to construct mixed effect model for this domain

Months after Dx	3	6	9	12	15	18	21	27	27	30	33	36	39	42	45	48	51	54	57	60
Total assessments
5-7 yrs	16	23	14	15	4	4	6	2	3	9	4	3	5	1	6	2	0	1	0	1
8-12 yrs	23	16	24	6	7	9	2	4	8	3	10	8	1	4	2	3	2	2	1	1
13-18 yrs	12	12	6	3	3	5	2	2	3	2	2	1	1	0	0	0	0	0	0	1

Thanks for your time and advance help

[Erik Ruzek]

The likelihood ratio test suggests that the model with the interactions provides a better fit (with the LR test, you divide the p-value by 2). When doing this kind of model comparison, you should ideally stick to using the same criteria for making judgements, whether that is the LR test, AIC, or BIC is really a matter of preference.

Could you please give me more suggestion in case my hypothesis is : is it better physical function (or emotion, social...) by different timeline in different age group

Are you trying to replicate the analyses shown in those graphs? If so, they are estimating something other than a linear relation to produce those graphs. Perhaps some kind of spline or maybe even a general additive model? At any rate, it seems that you need to run the same set of models as you did for physical on emotional and social. Then you can produce the graphs for those as you did with physical.

[Hang Tran]

Hi Erik,
Thanks for your answer... I'm so appreciate your patient with my difficutites

Are you trying to replicate the analyses shown in those graphs?

Those are all line graphs I did in my study desmontrated the dynamice change in Quality of life in children with acute lymphoblastic leukemia by using poly kernel smoothing.
At the first time of posting question, I don't know how to upload graph, and I'm so panic after construct and reconstruct model, my teacher said the result of model seems not fit with what the line graph showed, and the hypothesis underlying that graph
I think I should simplify problem by simple statement as your understand before

My guess is that you want to know whether the association between time from diagnosis (totaldate) and the outcome varies by child age category.

However, I want to seperately construct model for different age group
first I create dummy variables according to Age group, so I have 4 dummy variable for 4 age groups, then I run dummy1 (toddler)
I wonder should I seperate Physical into Age Group as

separate Physical_P, by(AgeGroup)

then I have four new variable sepearate from Physical_P by age group label as : Physical_P0 (aged 2-4), Physical_1 (age5-7)

Then I construct model for only toddler aged 2-4 as

1- construct with separated dependent variable according to age group
mixed Physical_P0 i.toddler totaldate immu education || ID:

Or should I construct model as

2- construct model for only toddler age 2-4 but no seperate dependent variable...
xtset ID totaldate
mixed Physical_P i.toddler totaldate immu education || ID:
eststo m2
mixed Physical_P i.toddler##c.totaldate immu education || ID:
eststo m3
lrtest m3 m2, stats
margins toddler, at(totaldate = (0(100)1000))
marginsplot

3- construct model for only 3 groups (5-7 vs 8-12) (8-12 vs 13-18) and (5-7 vs 13-18)
egen childado = group (child ado)
~ I have 2 dummy variable as child (aged 8-12) and ado (aged 13-18)
mixed Physical_P i.childado totaldate immu education || ID:
eststo m2
mixed Physical_P i.childado##c.totaldate immu education || ID:
eststo m3
lrtest m3 m2, stats

Am I on the right way fit to the aims as to see the association between time and QoL (i.e Physcial domain) in different pair of age groups, in case I choose group aged 8-12 for reference


Thanks for your previous suggestions, I leared more about regression coefficient in terms of interaction
However, I pretty confuses about that

mixed Physical_P i.agecat##c.totaldate immu education || ID: eststo m3 lrtest m3 m2, stats // no real improvement in model fit by including interactions

Could you please help me to understand why we need to compare model with/ without in term of interaction and choose one
How about comparison between model with control variables and model a random coefficient

mixed Physical_P youngc child ado totaldate immu education employment|| ID:
quietly mixed Physical_P youngc child ado totaldate immu education employment|| ID: totaldate, covariance(unstructured) nolog

Thank you very much for your help and supports

[Erik Ruzek]

The model I presented, based on your example data,

mixed Physical_P i.agecat##c.totaldate immu education || ID:

allows you to model all the children in the same model. The i.agecat##c.totaldate interaction gives you parameter estimates (coefficients) for each child group for the association between totaldate and Physical_P. This is the most parsimonious way to model the data and address your research question. The margins and marginsplot give you the model predicted Physical_P given different totaldates for each of the groups.

I compared this more complicated model with the interactions to a simpler model with no interactions. In this case, it is a joint test of whether the interactions are significantly different from 0. It is like the test of an interaction from an anova model. More generally we use this kind of testing to determine whether the inclusion of more parameters in the model improves the fit of the model to the data relative to a model with fewer parameters.

If you want to run this regression model separately for each age group, then you could do something like this,

foreach n if numlist 0/3 {
mixed Physical_P totaldate immu education || ID: if agecat==`n'
eststo mix_age_`n' }

which runs the model separately by agecat. They key thing to realize is that to restrict the sample for the regression, you use an if statement.

[Hang Tran]

Thanks for your concise explaination, I understood. I got another problem when I run mixed model for only toddler group aged 2-4
They said collinearity I don't know what I should do now @@. Could you please have a look and give me suggestions ..
Code:

xtset ID totaldate

Panel variable: ID (unbalanced)
Time variable: totaldate, 0 to 2582, but with gaps
Delta: 1 unit

mixed Physical_P totaldate immu education employment|| ID: if AgeGroup ==0


mixed Physical_P i.AgeGroup##c.totaldate immu education employment || ID: if AgeGroup ==0
note: 0.AgeGroup omitted because of collinearity.
note: 0.AgeGroup#c.totaldate omitted because of collinearity.
Performing EM optimization ...

Performing gradient-based optimization:
Iteration 0: Log likelihood = -530.31972
Iteration 1: Log likelihood = -530.29601
Iteration 2: Log likelihood = -530.296

Computing standard errors ...

Mixed-effects ML regression Number of obs = 119
Group variable: ID Number of groups = 88
Obs per group:
min = 1
avg = 1.4
max = 3
Wald chi2(4) = 8.94
Log likelihood = -530.296 Prob > chi2 = 0.0625


Physical_P Coefficient Std. err. z P>z [95% conf. interval]

0.AgeGroup 0 (omitted)
totaldate .0108845 .0060592 1.80 0.072 -.0009914 .0227604

AgeGroup#c.totaldate
0 0 (omitted)

immu 5.348036 6.711405 0.80 0.426 -7.806077 18.50215
education .4277383 4.633424 0.09 0.926 -8.653605 9.509082
employment -12.31359 5.654502 -2.18 0.029 -23.39621 -1.230967
_cons 73.28106 9.016 8.13 0.000 55.61003 90.9521



Random-effects parameters Estimate Std. err. [95% conf. interval]

ID: Identity
var(_cons) 241.9951 77.22085 129.4754 452.2991

var(Residual) 231.7099 56.72232 143.4073 374.3844

LR test vs. linear model: chibar2(01) = 10.93 Prob >= chibar2 = 0.0005

.
. eststo m3

.
. lrtest m3 m2, stats
df(unrestricted) = df(restricted) = 7
r(498);

Thank you very much for your time and support..

[Erik Ruzek]

Since you are restricting the regressions to just children from a single age category, there is no need for the interaction. You just include totaldate in the model. That is the totaldate slope for children in the age category you have restricted your sample to. See my code below.

recode CAge (2/4 = 0 "toddler") (5/7 = 1 "young child") (8/12 = 2 "child") ///
            (13/18 = 3 "adolescent"), gen(agecat)
*Mixed model for each agecat
mixed Physical_P totaldate immu education || ID: if agecat==0

*Or more parsimoniously using a foreach loop to loop through 0,1,2,3
foreach n of numlist 0/3 {
    mixed Physical_P totaldate immu education || ID: if agecat==`n'
    eststo mix_age_`n'
    }