use simulated_multilevel_longitudinal_data.dta5 Longitudinal Multilevel Models
5.1 The Data
The data employed in these examples are a longitudinal extension of the data described in Section 1.2.
5.2 The Equation
\[\text{outcome}_{itj} = \beta_0 + \beta_1 \text{parental warmth}_{itj} + \beta_2 \text{physical punishment}_{itj} + \beta_3 \text{time}_{itj} \ + \tag{5.1}\]
\[\beta_4 \text{identity}_{itj} + \beta_5 \text{intervention}_{itj} + \beta_6 \text{HDI}_{itj} +\]
\[u_{0j} + u_{1j} \times \text{parental warmth}_{itj} \ + \]
\[v_{0i} + v_{1i} \times \text{time}_{itj} + e_{itj}\]
5.3 Growth Trajectories
Remember, following Section 6.4, that in longitudinal multilevel models, the variable for time assumes an important role as we are often thinking of a growth trajectory over time.
As discussed in Section 6.4, think about a model where identity is a (1/0) variable for membership in one of two groups:
\[\text{outcome} = \beta_0 + \beta_t \text{time} + \beta_\text{identity} \text{identity} + \beta_\text{interaction} \text{identity} \times \text{time} + u_{0i} + e_{it}\]
Then, each identity group has its own intercept and time trajectory:
| Group | Intercept | Slope (Time Trajectory) |
|---|---|---|
| 0 | \(\beta_0\) | \(\beta_t\) |
| 1 | \(\beta_0 + \beta_\text{identity}\) | \(\beta_t + \beta_\text{interaction}\) |
Main Effects and Interactions
Thus, again following Section 6.4, in longitudinal multilevel models, main effects modify the intercept of the time trajectory, while interactions with time, modify the slope of the time trajectory. Below, we run models with main effects only, then models with main effects, and interactions with time.
5.4 Run Models
5.4.0.1 Get The Data
5.4.0.2 Run The Model
5.4.0.2.1 Main Effects Only
mixed outcome t warmth physical_punishment i.identity i.intervention HDI || country: warmthPerforming EM optimization ...
Performing gradient-based optimization:
Iteration 0: Log likelihood = -28739.506
Iteration 1: Log likelihood = -28739.506
Computing standard errors ...
Mixed-effects ML regression Number of obs = 9,000
Group variable: country Number of groups = 30
Obs per group:
min = 300
avg = 300.0
max = 300
Wald chi2(6) = 1119.81
Log likelihood = -28739.506 Prob > chi2 = 0.0000
-------------------------------------------------------------------------------------
outcome | Coefficient Std. err. z P>|z| [95% conf. interval]
--------------------+----------------------------------------------------------------
t | .9443446 .0756408 12.48 0.000 .7960914 1.092598
warmth | .9123903 .0430042 21.22 0.000 .8281035 .996677
physical_punishment | -.9881587 .0451732 -21.87 0.000 -1.076696 -.8996209
2.identity | -.1241465 .1242225 -1.00 0.318 -.367618 .1193251
1.intervention | .8575839 .1245179 6.89 0.000 .6135332 1.101635
HDI | -.0025173 .0191696 -0.13 0.896 -.0400891 .0350544
_cons | 50.54528 1.304146 38.76 0.000 47.9892 53.10136
-------------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects parameters | Estimate Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
country: Independent |
var(warmth) | .0229349 .0135353 .0072136 .0729194
var(_cons) | 3.0009 .8550708 1.716768 5.245553
-----------------------------+------------------------------------------------
var(Residual) | 34.31935 .5130963 33.3283 35.33988
------------------------------------------------------------------------------
LR test vs. linear model: chi2(2) = 767.22 Prob > chi2 = 0.0000
Note: LR test is conservative and provided only for reference.5.4.0.2.2 Interactions With Time
mixed outcome c.t##(c.warmth c.physical_punishment i.identity i.intervention c.HDI) || country: warmthPerforming EM optimization ...
Performing gradient-based optimization:
Iteration 0: Log likelihood = -28738.554
Iteration 1: Log likelihood = -28738.554
Computing standard errors ...
Mixed-effects ML regression Number of obs = 9,000
Group variable: country Number of groups = 30
Obs per group:
min = 300
avg = 300.0
max = 300
Wald chi2(11) = 1122.75
Log likelihood = -28738.554 Prob > chi2 = 0.0000
---------------------------------------------------------------------------------------
outcome | Coefficient Std. err. z P>|z| [95% conf. interval]
----------------------+----------------------------------------------------------------
t | .7537359 .3719996 2.03 0.043 .0246301 1.482842
warmth | .8198365 .0911059 9.00 0.000 .6412723 .9984008
physical_punishment | -1.000348 .1198049 -8.35 0.000 -1.235162 -.7655353
2.identity | -.2340191 .3271243 -0.72 0.474 -.875171 .4071327
1.intervention | .6597456 .3275877 2.01 0.044 .0176856 1.301806
HDI | -.0005531 .0210866 -0.03 0.979 -.041882 .0407757
|
c.t#c.warmth | .0463746 .0402459 1.15 0.249 -.0325059 .1252551
|
c.t#|
c.physical_punishment | .0061255 .0551491 0.11 0.912 -.1019647 .1142157
|
identity#c.t |
2 | .0548965 .1513015 0.36 0.717 -.241649 .3514421
|
intervention#c.t |
1 | .0990704 .151503 0.65 0.513 -.19787 .3960108
|
c.t#c.HDI | -.0009791 .0043888 -0.22 0.823 -.0095811 .0076229
|
_cons | 50.92503 1.494157 34.08 0.000 47.99654 53.85352
---------------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects parameters | Estimate Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
country: Independent |
var(warmth) | .0228292 .0135078 .0071588 .0728013
var(_cons) | 3.001849 .8552796 1.71738 5.247001
-----------------------------+------------------------------------------------
var(Residual) | 34.31227 .5129896 33.32141 35.33258
------------------------------------------------------------------------------
LR test vs. linear model: chi2(2) = 767.35 Prob > chi2 = 0.0000
Note: LR test is conservative and provided only for reference.5.4.0.3 Get The Data
library(haven)
dfL <- read_dta("simulated_multilevel_longitudinal_data.dta")5.4.0.4 Run The Model
Caution
lme4 does not directly provide p values in results, because of some disagreement over exactly how these p values should be calculated. Therefore, in this Appendix, I also call library lmerTest to provide p values for lme4 results.
Tip
R prefers to use scientific notation when possible. I find that the use of scientific notation can be confusing in reading results. I turn off scientific notation by setting a penalty for its use: options(scipen = 999).
5.4.0.4.1 Main Effects Only
library(lme4)
library(lmerTest)
options(scipen = 999)
fit2A <- lmer(outcome ~ t + warmth + physical_punishment +
identity + intervention + HDI +
(1 | country/id),
data = dfL)
summary(fit2A)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula:
outcome ~ t + warmth + physical_punishment + identity + intervention +
HDI + (1 | country/id)
Data: dfL
REML criterion at convergence: 57022.7
Scaled residuals:
Min 1Q Median 3Q Max
-3.6850 -0.6094 -0.0035 0.6133 3.6792
Random effects:
Groups Name Variance Std.Dev.
id:country (Intercept) 8.438 2.905
country (Intercept) 3.675 1.917
Residual 26.036 5.103
Number of obs: 9000, groups: id:country, 3000; country, 30
Fixed effects:
Estimate Std. Error df t value
(Intercept) 50.5161891 1.4296454 31.1737942 35.335
t 0.9433806 0.0658755 5998.3764616 14.321
warmth 0.9140307 0.0379336 4745.3496835 24.096
physical_punishment -1.0087537 0.0497972 6483.6770989 -20.257
identity -0.1319548 0.1517350 2968.7829265 -0.870
intervention 0.8591494 0.1520510 2971.8112001 5.650
HDI 0.0007909 0.0207656 28.0001836 0.038
Pr(>|t|)
(Intercept) < 0.0000000000000002 ***
t < 0.0000000000000002 ***
warmth < 0.0000000000000002 ***
physical_punishment < 0.0000000000000002 ***
identity 0.385
intervention 0.0000000175 ***
HDI 0.970
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) t warmth physc_ idntty intrvn
t -0.091
warmth -0.088 -0.002
physcl_pnsh -0.090 -0.007 -0.012
identity -0.156 0.000 -0.013 -0.003
interventin -0.055 0.000 0.039 0.019 -0.018
HDI -0.941 0.000 -0.004 0.005 0.000 0.0025.4.0.4.2 Interactions With Time
fit2B <- lmer(outcome ~ t *(warmth + physical_punishment +
identity + intervention + HDI) +
(1 | country/id),
data = dfL)
summary(fit2B)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula:
outcome ~ t * (warmth + physical_punishment + identity + intervention +
HDI) + (1 | country/id)
Data: dfL
REML criterion at convergence: 57042.8
Scaled residuals:
Min 1Q Median 3Q Max
-3.7118 -0.6092 -0.0024 0.6150 3.6779
Random effects:
Groups Name Variance Std.Dev.
id:country (Intercept) 8.436 2.905
country (Intercept) 3.675 1.917
Residual 26.046 5.104
Number of obs: 9000, groups: id:country, 3000; country, 30
Fixed effects:
Estimate Std. Error df t value
(Intercept) 51.0036725 1.6087742 49.9583024 31.703
t 0.6989769 0.3746882 6131.2125222 1.865
warmth 0.8170912 0.0805355 8274.9994610 10.146
physical_punishment -1.0097729 0.1113557 8084.6085126 -9.068
identity -0.2446453 0.3041604 8695.8966126 -0.804
intervention 0.6604671 0.3046286 8697.0843469 2.168
HDI 0.0026692 0.0221295 36.1037721 0.121
t:warmth 0.0486211 0.0356217 6404.8722416 1.365
t:physical_punishment 0.0004964 0.0494590 6753.0158846 0.010
t:identity 0.0563140 0.1318043 5993.4518199 0.427
t:intervention 0.0995037 0.1319917 5994.1433047 0.754
t:HDI -0.0009379 0.0038233 5993.9091197 -0.245
Pr(>|t|)
(Intercept) <0.0000000000000002 ***
t 0.0622 .
warmth <0.0000000000000002 ***
physical_punishment <0.0000000000000002 ***
identity 0.4212
intervention 0.0302 *
HDI 0.9047
t:warmth 0.1723
t:physical_punishment 0.9920
t:identity 0.6692
t:intervention 0.4510
t:HDI 0.8062
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) t warmth physc_ idntty intrvn HDI t:wrmt t:phy_
t -0.466
warmth -0.169 0.285
physcl_pnsh -0.183 0.313 -0.005
identity -0.278 0.450 -0.013 -0.002
interventin -0.100 0.162 0.039 0.019 -0.017
HDI -0.892 0.230 -0.007 0.012 -0.001 0.003
t:warmth 0.150 -0.324 -0.882 0.001 0.011 -0.035 0.006
t:physcl_pn 0.164 -0.351 0.004 -0.894 -0.001 -0.017 -0.010 -0.003
t:identity 0.242 -0.519 0.011 0.000 -0.867 0.014 0.001 -0.013 0.002
t:intervntn 0.087 -0.187 -0.035 -0.017 0.014 -0.867 -0.003 0.041 0.019
t:HDI 0.310 -0.666 0.015 -0.027 0.002 -0.007 -0.346 -0.016 0.029
t:dntt t:ntrv
t
warmth
physcl_pnsh
identity
interventin
HDI
t:warmth
t:physcl_pn
t:identity
t:intervntn -0.016
t:HDI -0.002 0.0085.4.0.5 Get The Data
using Tables, MixedModels, StatFiles, DataFrames, CategoricalArrays, DataFramesMeta
dfL = DataFrame(load("simulated_multilevel_longitudinal_data.dta"))5.4.0.6 Run The Model
5.4.0.6.1 Change Country To Categorical
@transform!(dfL, :country = categorical(:country))5.4.0.6.2 Main Effects Only
m2A = fit(MixedModel, @formula(outcome ~ t + warmth +
physical_punishment +
identity + intervention +
HDI +
(1 | country) +
(0 + warmth | country) +
(1 | id)), dfL)Linear mixed model fit by maximum likelihood
outcome ~ 1 + t + warmth + physical_punishment + identity + intervention + HDI + (1 | country) + (0 + warmth | country) + (1 | id)
logLik -2 logLik AIC AICc BIC
-28499.6031 56999.2063 57021.2063 57021.2356 57099.3610
Variance components:
Column Variance Std.Dev. Corr.
id (Intercept) 8.387351 2.896092
country (Intercept) 3.166939 1.779590
warmth 0.010760 0.103732 .
Residual 26.027290 5.101695
Number of obs: 9000; levels of grouping factors: 3000, 30
Fixed-effects parameters:
───────────────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
───────────────────────────────────────────────────────────────
(Intercept) 50.5949 1.35491 37.34 <1e-99
t 0.943864 0.0658716 14.33 <1e-45
warmth 0.913496 0.0423739 21.56 <1e-99
physical_punishment -1.0079 0.0497622 -20.25 <1e-90
identity -0.127692 0.151584 -0.84 0.3996
intervention 0.858997 0.15191 5.65 <1e-07
HDI -0.000565882 0.0196433 -0.03 0.9770
───────────────────────────────────────────────────────────────5.4.0.6.3 Interactions With Time
m2B = fit(MixedModel, @formula(outcome ~ t * (warmth +
physical_punishment +
identity + intervention +
HDI) +
(1 | country) +
(0 + warmth | country) +
(1 | id)), dfL)Linear mixed model fit by maximum likelihood
outcome ~ 1 + t + warmth + physical_punishment + identity + intervention + HDI + t & warmth + t & physical_punishment + t & identity + t & intervention + t & HDI + (1 | country) + (0 + warmth | country) + (1 | id)
logLik -2 logLik AIC AICc BIC
-28498.3091 56996.6182 57028.6182 57028.6788 57142.2979
Variance components:
Column Variance Std.Dev. Corr.
id (Intercept) 8.391746 2.896851
country (Intercept) 3.170031 1.780458
warmth 0.010609 0.102999 .
Residual 26.015906 5.100579
Number of obs: 9000; levels of grouping factors: 3000, 30
Fixed-effects parameters:
──────────────────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
──────────────────────────────────────────────────────────────────
(Intercept) 51.0751 1.54284 33.10 <1e-99
t 0.702771 0.374539 1.88 0.0606
warmth 0.817076 0.0826636 9.88 <1e-22
physical_punishment -1.00903 0.111293 -9.07 <1e-18
identity -0.238714 0.303996 -0.79 0.4323
intervention 0.660761 0.30445 2.17 0.0300
HDI 0.00136065 0.0210842 0.06 0.9485
t & warmth 0.0483635 0.0356074 1.36 0.1744
t & physical_punishment 0.000542203 0.0494355 0.01 0.9912
t & identity 0.0554385 0.131745 0.42 0.6739
t & intervention 0.0992809 0.131925 0.75 0.4517
t & HDI -0.000955067 0.00382162 -0.25 0.8027
──────────────────────────────────────────────────────────────────5.5 Interpretation
The main effects only model suggests that time is associated with increases in the outcome. In the main effects model, main effects other than time, indicate whether a particular variable is associated with higher or lower intercepts of the time trajectory, at the beginning of the study time. Warmth is again associated with increases in the outcome, while physical punishment is associated with decreases in the outcome. Identity is again not associated with the outcome, while the intervention is associated with higher levels of the outcome. The Human Development Index is again not associated with the outcome.
The second model adds interactions with time to the first model. Results are largly similar to the prior model. However, here we not only examine whether main effects other than time are associated with higher or lower time trajectories, but also whether particular variables are associated with differences in the slope of the time trajectory. In this case, we find that no independent variable is associated with changes in the slope of the time trajectory.
However, it may be illustrative to imagine how we would interpret the results had a particular interaction term been statistically significant. Let us consider one of the interaction terms with the largest coefficient, intervention#time. The interaction of the intervention with time is positive. Had this coefficient been statistically signifcant, it would have indicated that the intervention was associated with more rapid increases in the outcome over time in addition to the fact that the intervention is associated with higher initial levels of the outcome.