A Two-Step Robust Estimation Approachfor Inferring Within-person Relations inLongitudinal Design: Tutorial and
Simulations
Satoshi Usami
Department of Education, University of Tokyo
usami s@p.u-tokyo.ac.jp
January 26, 2023

1

A Two-Step Robust Estimation Approach forInferring Within-Person Relations in Longitudinal
Design: Tutorial and Simulations

Abstract
Psychological researchers have shown an interest in disaggregating within-person
variability from between-person diﬀerences. Especially for inferring reciprocal relations among variables at the within-person level, applications of the random-intercept
cross-lagged panel model (RI-CLPM) with stable trait factors has increased rapidly.
This paper provides a tutorial, simulation, and illustrative example of another recent
approach proposed by Usami (2022a). This approach consists of a two-step procedure: within-person variability scores (WPVS) for each person, which are disaggregated from the stable traits of that person, are calculated using structural equation
modeling, and causal parameters are then estimated via a potential outcome approach, such as by using structural nested mean models (SNMMs). This method
assumes a data-generating process similar to that in RI-CLPM, and has several advantages: (i) the flexible inclusion of curvilinear and interaction eﬀects for WPVS
as latent variables, (ii) more accurate estimates of causal parameters for reciprocal
relations can be obtained under certain conditions owing to them being doubly robust, even if unobserved time-varying confounders and model misspecifications exist,
and (iii) the risk of obtaining improper solutions is minimized. After explaining thedata-generating process and the analysis procedure using the R package DTRreg for
SNMMs, estimation performances are compared with RI-CLPM through large-scale
simulations. We show that the proposed approach works well in many conditions
if longitudinal data with T ≥ 4 are available, and that the accuracy increases as T
becomes larger. An analytic example using data regarding sleep habits and mental
health statuses from the Tokyo Teen Cohort (TTC) study is also provided.

Keywords: within-person relation, longitudinal data, structural equation modeling, structural nested mean model, causal inference1

1

INTRODUCTION

When analyzing relations among variables in data, psychological researchers diﬀerentiate
between within-group (unit-level) relations and between-group (group-level) relations, andbetween within-person relations and between-person relations. Particularly in longitudinal
design, researchers have shown an interest in inferring within-person relations: how changes
in one variable influence another for the same person. Within-person relations may exhibit
statistically diﬀerent (or even opposite) tendencies from between-person relations, and thisis one reason that statistical inference for disaggregating within- and between-person relations has long been a concern in psychology. On the other hand, estimands that are definedat the within-person level have been less common in the causal inference literature (Lüdtke& Robitzsch, 2021).
Multilevel modeling (Hoﬀman, 2014; Wang & Maxwell, 2015) and structural equation
modeling (SEM) (Hamaker, Kuiper, & Grasman, 2015) are two statistical methods thathave been popular for investigating within-person relations. The SEM approach might beadvantageous if researchers wish to (i) include common factors in measurement models
for multiple indicators, (ii) assume measurement errors for measurements with imperfectreliability, (iii) treat multiple outcomes to evaluate reciprocal (or mutual) relations, or (iv)use model fit indices to evaluate how the model implied mean and covariance structure canreproduce the observed mean vector and covariance matrix. After the critique of applying
a cross-lagged panel model (CLPM) to infer reciprocal relations at the within-person level(Hamaker et al., 2015), applications of random-intercept CLPM (RI-CLPM), which includes common factors called stable trait factors to control for between-person diﬀerences,have been growing rapidly. It has been empirically shown that model choice of CLPM or
RI-CLPM can be critical in terms of the signs, statistical significance, and magnitudes ofkey parameter estimates (i.e., cross-lagged parameters) for reciprocal relations (e.g., Orth,Clark, Donnellan, & Robins, 2021; Usami, Murayama, & Hamaker, 2019). RI-CLPM is a2

useful analytic option, but various kinds of (SEM-based) statistical models are available
for examining reciprocal relations, and model choice is still an ongoing issue (see Andersen,
2022; Lucas, 2022; Lüdtke & Robitzsch, 2021; Usami, 2021,2022ab; Usami, Murayama, etal., 2019 and later discussion).
Along with these increasing applications and theoretical interests, Usami (2022a) aimed
to synthesize the SEM-based approach traditionally used in psychology and the potential
outcome approaches used in epidemiology, in order to enable flexible and robust inferenceof within-person relations. This method consists of a two-step procedure that assumes thesimilar data-generating process (DGP) as assumed in the RI-CLPM. In this approach, thewithin-person variability scores (WPVS) for all participants, which are disaggregated fromtheir stable trait factor scores, are first calculated in each variable using SEM. Then, causal
parameters modeled at the within-person level are estimated using a potential outcomesapproach with the calculated WPVS, such as marginal structural models (MSMs; Robins,
1999; Robins, Hernán, & Brumback, 2000) or structural nested (mean) models (SNMMs;
Robins, 1989, 1992). Usami (2022a) originally assumed a situation where researchers sought
to investigate the unidirectional relation between variables, but it is straightforward toextend the approach to reciprocal relations. For brevity, we will refer to this two-stepestimation approach as the TS method.
The TS method has some advantages over RI-CLPM. If researchers assume a DGP that
includes stable trait factors, then they need to model the relations among WPVS as latentvariables. This corresponds to the fact that RI-CLPM is sometimes classified as a deviation
model rather than an observation model (e.g., Andersen, 2022). However, in the standard
SEM as covariance structure analysis, estimating curvilinear (e.g., quadratic) eﬀects andinteraction eﬀects for latent variables can be challenging in terms of implementation, especially if multiple indicators are not available for each variable at every time point (e.g.,
Kenny & Judd, 1984; Klein & Moosbrugger, 2000; Kline, 2016 pp.437-444; Marsh, Wen,
& Hau, 2006). Even if multiple indicators are available, identification may require strong3

parameter constraints, and implementation could be very diﬃcult, especially for cases inwhich researchers wish to include interaction eﬀects between variables from diﬀerent time
points (e.g., X3∗ Y1∗ for variables X ∗ at t = 3 and Y ∗ at t = 1 at the within-person level).
Bayesian approach is a useful alternative (e.g., Lee, 2007; Ozkok et al., 2022), but plausible
specifications of priors and computation time can be challenging. For these reasons, usualapplications of RI-CLPM do not assume curvilinear and interaction eﬀects of variables at
the within-person level (WPVS). On the other hand, in TS method, such eﬀects can be
introduced in a straightforward manner because the calculated WPVS are directly used forcausal parameter estimation.
Another potential advantage of the TS method is the manner in which it accounts
for time-varying confounders from the view of causal inference. RI-CLPM (and otherSEM-based models) was originally proposed to uncover reciprocal relations between focal variables at the within-person level (X ∗ and Y ∗ ), and it does not directly account for
time-varying confounders. One reasonable procedure to account for observed time-varying
confounders in RI-CLPM is to simply include them in the lagged regression (i.e., structuralmodels) for reciprocal relations like ANCOVA. However, along with the potential aforementioned restrictions in modeling curvilinear and interaction eﬀects for observed time-varying
confounders, this approach requires researchers to correctly specify all structural models
at each time point, and parameters might be seriously biased if researchers miss some
unobserved time-varying confounders. On the other hand, as we will illustrate, a more accurate estimate of causal parameters can be obtained under certain conditions using the TS
method, especially with (linear) SNMMs (we call these TS-SNMMs later for brevity), even
if unobserved time-varying confounders and model misspecifications relating to observed
confounders exist.
The TS method is also potentially advantageous for handling improper solutions. Though
this risk of these occurring in RI-CLPM is empirically known to be smaller than in some
SEM-based approaches, it can still frequently produce improper solutions (e.g., the non-4

Hessian matrix of trait factor scores and negative variances), especially when either the
sample size or the number of time points is small (Orth et al., 2021; Usami, Todo, & Murayama., 2019). The TS method takes a two-step estimation approach, so not all parameters
are simultaneously estimated as in SEM-based approaches. This can minimize the risk ofobtaining improper solutions caused by sample fluctuations and model misspecifications.
For these reasons, the TS method can be considered as a viable alternative of RICLPM (and its variants) for inferring within-person (reciprocal) relations, especially when
researchers wish to account for time-varying confounders. Comparing estimation results
between these approaches could also be useful as a kind of sensitivity analysis. However,
MSMs and SNMMs have been developed in epidemiology and we have seen few applications of these methods in psychology, while useful introductions and applications of MSMs
are increasing (e.g., VanderWeele, Hawkley, Thisted & Cacioppo, 2011). Actual applications of SNMMs have also been infrequent in other disciplines because of limited software
availability and tutorials (Vansteelandt & Joﬀe, 2014), as well as its typically strongly theoretical presentation and challenging implementation (Wallace, Moodie & Stephens, 2017a).
SNMs, however, are suitable and robust for handling violation of the usual assumptions of
no unobserved confounders and sequential ignorability (Robins, 1999; Robins & Hernan,
2009; Vansteelandt & Joﬀe, 2014). Only a brief introduction of MSMs and SNMMs wasprovided in Usami (2022a), and there are not many applied researchers in psychology whoare capable of applying such novel methods.
On a related note, another important limitation of Usami (2022a) is that their estimation performance was assessed under restricted scenarios of model (mis)specifications withlimited conditions for the number of time points and the magnitudes of WPVS and stable
trait factor (co)variances. WPVS need to be calculated in the first step of the TS method,which may result in less accurate causal estimates, especially when the number of time
points is limited under a misspecified measurement model. One of the primary motivations
for performing the simulation in Usami (2022a) was to compare the performance of the TS5method (which aims to conduct true score centering [e.g., Asparouhov & Muthén, 2018]
based on stable trait factors) over diﬀerent centering methods (i.e., observed person-meancentering). Therefore, how the estimation performances (e.g., the bias of causal parameters
and the frequency of improper solutions) compare between the TS method and RI-CLPM
is still an open question.
In this paper, we provide a tutorial, simulation, and illustrative example of how to use
the TS method to infer within-person (reciprocal) relations, focusing especially on using TSSNMMs to estimate causal parameters by utilizing the R package DTRreg (Wallace et al.,2017ab), which was not used in Usami (2022a). We demonstrate in large-scale simulations
that TS-SNMMs can flexibly and accurately estimate causal parameters due to them being
doubly robust, even if unobserved time-varying confounders and model misspecifications
exist. We also show that TS-SNMMs work well in many conditions if longitudinal datawith T ≥ 4 are available, while accuracy increases as T becomes larger. The R code for
TS-SNMMs is available in a supplementary document to make the method more accessible
for applied researchers.
The remainder of this paper is organized as follows. In Section 2 we start our discussion
by introducing assumed DGPs and the definition of causal eﬀects, while appropriately
referring to Usami (2022a). Readers familiar with these may skip this (sub)section. TSSNMMs are introduced in Section 3, and simulations are provided in Section 4. Section 5
describes an empirical application of TS-SNMMs using data from the Tokyo Teen Cohort(TTC) study. A summary is provided in the final section along with a discussion of our
future research agenda.

6

2

ASSUMED DATA-GENERATING PROCESS AND
DEFINITIONS OF CAUSAL EFFECT

2.1

Data-generating process

We suppose that data are generated at fixed time points and let Xit and Yit denote focal(continuous) variables at t (t = 1, . . . , T ) for person i: researchers wish to infer their
within-person relation. Also, let Lit be the observed time-varying confounders. We assume
a single confounder here for explanation purposes. Suppose that a time-varying confounder
has three characteristics: it is independently associated with future focal variables as well
as future confounders, and it is aﬀected by earlier focal variables and confounders.
Figure 1 is a directed acyclic graph (DAG) that expresses the causal relations among
variables in the assumed DGP of T = 4. Because we suppose the presence of stable trait
factors in DGP, as in RI-CLPM, this DAG includes stable trait factors I as time-invariant
factors. This DAG is similar to the one presented in Usami (2022a, Figure 1b), but it
now assumes that Xit , Yit , and Lit are measured at the same time, and thus, no directcausal relation is assumed among them within each time point. This setting is analogous
to the path diagrams of the SEM-based statistical models (e.g., RI-CLPM). Each solid
single-headed arrow represents a direct causal relation, and a dashed double-headed arrow
indicates the covarying relations due to unobserved confounding. Dashed circles are used
to express latent variables. To keep the illustration simple, here we temporarily assume
first-order linear lagged eﬀects of variables.
By the properties of stable trait factors (specifically that the diﬀerence between the
expected value of measurement of each person and the temporal group mean is invariant
over time; see Usami, 2022a), these factors have zero means (E(I (Y ) ) = E(I (A) ) = E(I (L) ) =0) and additively influence respective measurements. For the same reason, the values of
the coeﬃcients corresponding to the paths from these factors to the measurements are7

restricted to one.
In this DAG, directed edges from WPVS, which are expressed by variable names withan asterisk (e.g., Y3∗ ), are drawn to the corresponding measurements. Directed edges are
assumed between WPVS as (time-varying) latent variables rather than between measurements. Without loss of generality, we assume that WPVS for all variables have zero means.
It is important to note that WPVS are also assumed to be uncorrelated with stable trait
factors (we will revisit this issue in the Discussion). As a result, like RI-CLPM, stable
trait factors I have only direct eﬀects on measurements, and each measurement can be
decomposed into a (linear) sum of time-invariant (stable traits) and time-varying factors(i.e., WPVS) that are mutually uncorrelated.
Under the linear causal DAG model, the DGP can be represented by the following
equations, which correspond to RI-CLPM that includes L:(Y )Yit = µt

(Y )

+ Ii

+ Yit∗ , Xit = µt

(X)

(X)

+ Ii

+ Xit∗ , Lit = µt

(L)

(L)

+ Ii

+ L∗it(1)

for t ≥ 1, and∗∗Yit∗ = αt Yi(t−1)+ βt Xi(t−1)+ γt L∗i(t−1) + dit
(Y )

(Y )

∗
Xit∗ = αt Yi(t−1)
+ βt
(X)

(Y )

(X)

∗
Xi(t−1)
+ γt

(X)

(Y )

L∗i(t−1) + dit

(X)

(2)

∗
∗
L∗it = αt Yi(t−1)+ βt Xi(t−1)+ γt L∗i(t−1) + dit
(L)

(L)

(L)

(L)

for t ≥ 2. µ and d denote the temporal group means and residual terms, respectively, andthey are omitted in the DAG representation. As noted, stable trait factors are assumed to
be uncorrelated with WPVS. For example,
Cov(Ii , Yit∗ ) = 0,(Y )Cov(Ii , Xit∗ ) = 0,(Y )Cov(Ii , L∗it ) = 0
(Y )

(3)

(Y )

for a stable trait factor of a variable Y (Ii ).
As implied in Equation (1), under these specifications, Yit∗ , Xit∗ , and L∗it represent tem(Y )poral deviations from the expected score for person i at t (i.e., µt

8

(Y )

(X)

+ Ii , µt

(X)

+ Ii

,

(L)

(L)

and µt +Ii ), whereas the stable trait factors represent stable between-person diﬀerences
over time. To put it another way, WPVS can be characterized as the diﬀerence between ameasurement and its expected value for each person at each time point ( see Usami (2022a)
for further explanation). Additionally, (co)variances of measurements at t can be expressedas the sum of those of stable trait factor scores and WPVS. For a variable Y ,Cov(Yit , Yit′ ) = Cov(Yit∗ , Yit∗′ ) + V ar(Ii ).
(Y )

(4)

In RI-CLPM, the initial deviations (WPVS) are modeled as exogeneous variables, and theirvariances and covariances are estimated.
RI-CLPM can identify the parameters from data with T ≥ 3. If one is interested in∗inferring first-order lagged eﬀects of X ∗ on Y ∗ (Xt−1→ Yt∗ ), βt(Y )in Equation (2), whichis called a cross-lagged parameter, is key.

2.2

Definition of causal eﬀects at the within-person level

Below we assume a similar causal DAG model to that in Figure 1. However, unlike theprevious subsection, we relax some assumptions about WPVS, specifically that higherorder, curvilinear (e.g., quadratic), and interaction eﬀects for WPVS can exist. The current
focus is on evaluating the within-person relation between variables, that is, how one variable∗) influences another variable (e.g., Yt∗ ) and vice versa at the within-person level.
(e.g., Xt−1

Stable trait factors and WPVS are mutually uncorrelated, and each measurement canbe characterized as a collider in causal DAG (i.e., it is causally influenced by stable trait
factors and WPVS; Usami, 2022a), so causal eﬀects of WPVS can be defined independently
from stable trait factors. However, modeling stable trait factors through measurements isstill required for estimating causal parameters because WPVS are latent variables (Usami,
2022a). Below we first explain the definition of causal eﬀects of X ∗ (as predictor) on Y ∗(as outcome).
We use overbars Ȳt∗ = {Y1∗ , . . . , Yt∗ } to denote the history of Y ∗ through t and underbars

9∗Y ∗t = {Yt∗ , . . . , YT∗ } to denote the future of this variable. Let Yit∗ X̄i(t−1) (t = 1, . . . , T ) denotethe WPVS for the outcome that would take at time t for person i if this person had a∗∗∗history of predictors X ∗ at the within-person level X̄i(t−1)= {Xi1, . . . , Xi(t−1)} through∗t − 1. Yit∗ X̄i(t−1) is a potential outcome, which we connect to WPVS by the consistency
assumption∗

Yit∗ = Yit∗ x̄i(t−1)

(5)

∗

∗
if X̄i(t−1)= x̄∗i(t−1) ; otherwise, Yit∗ x̄i(t−1) is counterfactual. Note that Yit∗ is unobservable,

and thus can be predicted in the first step of the TS method, while potential outcomes formeasurements (i.e., Yit ) are considered in the standard potential outcome approach.
In the potential outcome approach, causal eﬀect refers to the contrast between potential
outcomes under diﬀerent levels of treatments/predictors (i.e. X ∗ ). The average causal eﬀect∗
on Yit∗ when Xi(t−1)
increases one unit from the reference value x∗ri(t−1) at time t − 1 can be

expressed as
∗

∗r

∗

∗r

∗

∗r

∗

∗r

E(Yit∗ x̄i(t−2) ,xi(t−1) +1 − Yit∗ x̄i(t−2) ,xi(t−1) ) = E(Yit∗ x̄i(t−2) ,xi(t−1) +1 ) − E(Yit∗ x̄i(t−2) ,xi(t−1) ).

(6)

The standard assumption of no unobserved confounders (or sequential ignorability) indicates that∗∗∗∗Y ∗it x̄i(t−2) ,0 ⊥⊥ Xi(t−1)|Ȳi(t−1), L̄∗i(t−1) , X̄i(t−2)= x̄∗i(t−2) .

(7)Here, (x̄∗i(t−2) , 0) is the counterfactual history, that is, the history that agrees with x̄∗i(t−2)
through time t−2 and is zero thereafter. Equation (7) indicates that the potential outcome∗∗, given the ob= 0) is independent of Xi(t−1)
at time t (were the person receives Xi(t−1)∗∗
served confounders in past time (X̄i(t−2), Ȳi(t−1), L̄∗i(t−1) )1 . Along with the assumed causal

DAG above, as well as the consistency and sequential ignorability, we impose the stable
1∗∗and X̄i(t−2)
are the outcome and predictors, they can be characterized as confounders
Though Ȳi(t−1)∗
when it causes Yit∗ , and if it either causes or associates with Xi(t−1).

10unit treatment value assumption (SUTVA; no unmodeled spillovers, e.g., Hong, 2015) andassumptions of positivity (i.e., the probability of taking each level of predictor (X ∗ ) conditional on past variables is greater than zero) and modularity (i.e., mechanisms that are
not directly targeted by treatment/predictors are not altered). Under these assumptions,the average causal eﬀect in Equation (6) can be expressed by the diﬀerence in conditional
means given by information on observed confounders as
∗

∗r

∗

∗r

E(Yit∗ x̄i(t−2) ,xi(t−1) +1 ) − E(Yit∗ x̄i(t−2) ,xi(t−1) )∗∗∗=E(Yit∗ |Ȳi(t−1), L̄∗i(t−1) , X̄i(t−2)= x̄∗i(t−2) , Xi(t−1)= x∗ri(t−1) + 1)∗∗∗
− E(Yit∗ |Ȳi(t−1), L̄∗i(t−1) , X̄i(t−2)= x̄∗i(t−2) , Xi(t−1)= x∗ri(t−1) ).

(8)

In other words, the average causal eﬀect can be evaluated from the diﬀerence in conditional
∗r
∗
= x∗r
means of Yit∗ between persons who receive Xi(t−1)
i(t−1) +1 (i.e., levels that are xi(t−1) +1
(X)

(X)

larger than their expected scores µ(t−1) + Ii∗
) and Xi(t−1)= x∗ri(t−1) , given information on

the observed confounders’ history.
∗
on Yit∗ when they
Similarly, the average joint eﬀects for a sequence of predictors X̄i(t−1)

commonly increase one unit from the reference values x̄∗ri(t−1) can be generally expressed as∗r∗r

E(Yit∗ x̄i(t−1) +1 ) − E(Yit∗ x̄i(t−1) )
∗
∗
∗
∗
∗
∗
∗r
=E(Yit∗ |Ȳi(t−1), L̄∗i(t−1) , X̄i(t−1)= x̄∗ri(t−1) + 1) − E(Yit |Ȳi(t−1) , L̄i(t−1) , X̄i(t−1) = x̄i(t−1) ). (9)Suppose T = 3, and that the DGP can be represented by linear and first-order regression
models, as in Equations (1) and (2) (assuming no interaction eﬀects for WPVS). Then, a∗= x̄∗i2 ) at t = 3 can be expressed as a linear (weighted)
conditional mean E(Yi3∗ |Ȳi2∗ , L̄∗i2 , X̄i2

11

sum of the terms x∗i1 and x∗i2 :E(α3 Yi2∗ + β3 x∗i2 + γ3 L∗i2 + di3 )
(Y )

(Y )

(Y )

(Y )

=α3 (α2 E(Yi1∗ ) + β2 x∗i1 + γ2 E(L∗i1 )) + β3 x∗i2 + γ3 (α2 E(Yi1∗ ) + β2 x∗i1 + γ2 E(L∗i1 ))
(Y )

(Y )

(Y )

(Y )

(Y )

(Y )

(L)

(L)

(L)

= [α3 (α2 E(Yi1∗ ) + γ2 E(L∗i1 )) + γ3 (α2 E(Yi1∗ ) + γ2 E(L∗i1 ))]
|
{z
}
(Y )

(Y )

(Y )

(Y )

(L)

(L)

0

+ [α3 β2 + γ3 β2 ] x∗i1 + β3 x∗i2
|{z}
|
{z
}
(Y ) (Y )(Y ) (L)

(Y )

β32

β31

=β31 x∗i1 + β32 x∗i2 .

(10)

The first term in the third line becomes zero because WPVS have zero means. From∗∗
this result, the joint eﬀects of Xi1
and Xi2when commonly increasing one unit from the∗rreference values x∗r
i1 and xi2 become∗r∗r∗r(β31 (x∗ri1 + 1) + β32 (xi2 + 1)) − (β31 xi1 + β32 xi2 ) = β31 + β32 .
(Y ) (Y )Note that β31 (= α3 β2(11)+ γ3 β2 ; the eﬀect of X1∗ on Y3∗ ) can also be evaluated by(Y ) (L)tracing the two paths X1∗ → Y2∗ → Y3∗ (= α3 β2 ) and X1∗ → L∗2 → Y3∗ (= γ3 β2 ) that(Y ) (Y )(Y ) (L)start at X1∗ and end at Y3∗ , shown in Figure 1.
∗= x∗i1 ) at t = 2 canLikewise, when the eﬀect of X1∗ on Y2∗ is a focus, E(Yi2∗ |Yi1∗ , L∗i1 , Xi1

be expressed asE(α2 Yi1∗ + β2 x∗i1 + γ2 L∗i1 ) = [α2 E(Yi1∗ ) + γ2 E(L∗i1 )] + β2 x∗i1 = β21 x∗i1 ,|{z} |{z}
(Y )

(Y )

(Y )

(Y )

(Y )

0

(Y )

(12)

β21

∗
when increasing one unit from the reference value x∗rand thus the causal eﬀect of Xi1i1 at the(Y )within-person level becomes β21 = β2 , which is equivalent to the cross-lagged parameterin Equation (2). Note that because stable trait factors are not associated with the WPVS
in the assumed DGP, the eﬀects of X ∗ on Y ∗ can be numerically equivalent to those on Y
(i.e., measurements) (Usami, 2022a).
Compared with other disciplines, such as epidemiology, investigations of joint eﬀects in
psychology have been limited in applications of SEM-based approaches (e.g., RI-CLPM)

12

where only the cross-lagged parameters gather attention from researchers. However, useful
information can be gained by interpreting the eﬀects of time-varying treatments/predictors∗∗(e.g., how a change in X(t−2)
at two previous time points influences Yt∗ , fixing X(t−1)).

In investigating reciprocal relations, the causal eﬀect of Y ∗ on X ∗ is also a focus. In a∗∗similar manner, the average causal eﬀect of Yi(t−1)
on Xit∗ when Yi(t−1)
increases one unit∗rfrom the reference value yi(t−1)
at t − 1 can be expressed as
∗

∗r

∗

∗r

E(Xit∗ ȳi(t−2) ,yi(t−1) +1 − Xit∗ ȳi(t−2) ,yi(t−1) ).

(13)∗
The average joint eﬀects of a sequence of Ȳi(t−1)
on Xit∗ when they commonly increase one∗runit from the reference values ȳi(t−1)
can be expressed as∗r∗r

E(Xit∗ ȳi(t−1) +1 ) − E(Xit∗ ȳi(t−1) )
∗
∗
∗r
∗
∗
∗r
=E(Xit∗ |X̄i(t−1), L̄∗i(t−1) , Ȳi(t−1)= ȳi(t−1)+ 1) − E(Xit∗ |X̄i(t−1), L̄∗i(t−1) , Ȳi(t−1)= ȳi(t−1)). (14)

We can now summarize the assumptions for identifying causal parameters in the TS
method: (i) measurements (T ≥ 3) are expressed by a linear sum of stable trait factors
and WPVS that are mutually uncorrelated, (ii) WPVS are expressed by functions of thosein past time, (iii) consistency, (iv) no unobserved confounders, (v) SUTVA, (vi) positivity,and (vii) modularity. The first assumption is unique to the TS method and RI-CLPM, andwe will revisit this issue in the Discussion.
Regarding the second assumption, if the DGP can be represented by linear equations,such as in Equations (1) and (2) (assuming no interaction eﬀects for WPVS), then RICLPM (that includes L) can identify causal parameters. However, modeling curvilinear(e.g., quadratic) and interaction eﬀects for WPVS as latent variables can be challenging
in RI-CLPM. Furthermore, all structural models (as in Equation (2)) should be correctly
specified, which might be very restrictive. The risk of obtaining improper solutions canalso be an obstacle when applying RI-CLPM. TS-SNMMs can overcome these potential
problems: they can provide a more flexible and robust inferential framework toward model13

misspecifications, resulting in more accurate estimates of causal parameters with a lower
risk of improper solutions.

3

INTRODUCTION OF THE TS METHODIn TS-SNMMs, WPVS are first calculated for each variable by SEM that only models
measurement parts with stable trait factors. Causal parameters in the structural model
are then estimated by SNMMs using calculated WPVS. Usami (2022a) argues that thiskind of two-step approach has some strengths compared with simultaneous estimation, suchas (i) minimizing the risk of potential confounding in interpreting the estimation results
when either the measurement or the structural model is misspecified, (ii) having a greater
feasibility for analyses (because common factors or unit eﬀects are not explicitly modeledin SNMMs), and (iii) having less risk of obtaining improper solutions. Regarding the firstpoint, when a simultaneous estimation procedure like RI-CLPM is used, misspecificationin the structural models at the within-person level may greatly aﬀect parameter estimatesin the measurement model ((co)variances of stable factors and WPVS), and vice versa.

3.1

Step 1: Estimation of measurement models and predictionof WPVSThe first step is subdivided into two sub-steps: (i) specification and estimation of themeasurement models and (ii) prediction of WPVS.
3.1.1

Specification and estimations of measurement models

The RI-CLPM and TS approaches commonly assume measurement models like Equation (1). This equation can be viewed as a model similar to the factor analysis model,
which includes a single common factor I (whose factor loadings are all one) and a unique14

factor in the form of WPVS. In vector notation, this equation for outcome Y becomes
Yi = µ(Y ) + Ii 1T + Yi∗ ,
(Y )

(Y )

(15)
(Y )

where µ(Y ) is a T × 1 mean vector, E(Ii ) = 0, V ar(Ii ) = ϕ2(Y ) , and 1T denotes the
T × 1 vector whose elements are all one. Yi∗ is a T × 1 vector of WPVS, and E(Yi∗ ) = 0
and Cov(Ii , Yit∗ ) = 0. We denote as Ψ(Y ) a T × T variance-covariance matrix of WPVS.
(Y )This implies that the variance-covariance matrix of Y (denoted as Σ(Y ) ) is of the form
Σ(Y ) = ϕ2(Y ) 1T 1tT + Ψ(Y ) .
Unlike the standard factor analysis model, Ψ(Y ) has a dependence structure and is not
diagonal. Therefore, in using SEM to estimate the parameters in Equation (15), somestructure, such as the autoregressive (AR) structure, must be specified in Ψ(Y ) to enablemodel identification. The specified model can then be diagnosed via model fit indices.
Similarly, we also set measurement models for the other variables X and L separately
in this sub-step, then estimate parameters for the mean vectors (µ(X) and µ(L) ), stable
trait factor variances (ϕ2(X) and ϕ2(L) ), and variance-covariance matrices of WPVS (Ψ(X)
and Ψ(L) ).
3.1.2

Predicting WPVSt
Let Zi = (Yit , Xit , Lti )t and Zi∗ = (Yi∗t , Xi∗t , L∗ti ) be vectors of measurements and WPVS,

respectively, and let µ = (µ(Y )t , µ(X)t , µ(L)t )t be a mean vector. Also, let Σ and Ψ becovariance matrices for measurements Zi and WPVS Zi∗ .
We consider linear prediction of WPVS Ẑi∗ under the condition that Σ and Ψ are known.
Consider a 3T × 3T weight matrix W that provides WPVS from measurementsẐi∗ = W t (Zi − µ),(16)E(Ẑi∗ Ẑi∗t ) = W t E[(Zi − µ)(Zi − µ)t ]W = W t ΣW = Ψ.

(17)satisfying the relation

15

Unlike standard applications of factor analysis, we are interested in predicting WPVS
(unique factor scores) rather than stable trait factor (common factor) scores. With thispoint in mind, the W that can minimize the risk function defined as the trace of a residual
covariance matrix (i.e., the mean squared error MSE(Ẑi∗ )=E[(Ẑi∗ − Zi∗ )t (Ẑi∗ − Zi∗ )]) andalso satisfy the relation in Equation (17) can be directly obtained by the so-called (linear)correlation preserving predictor (e.g., ten Berge, Krijinen, Wansbeek, & Shapiro, 1999
pp.317)−1/2

W t = Ψ1/2 (Ψ3/2 Σ−1 Ψ3/2 )Ψ3/2 Σ−1 .

(18)Here, for a positive (semi)definite matrix C, we denote as C 1/2 the positive (semi)definite
matrix whose square equals C. Matrices C −1/2 and C 3/2 are the inverse (if it exists) andthird power of C 1/2 , respectively.
We use the sample means Z̄ and the unbiased variance-covariance matrix S as estimatorsof µ and Σ. As implied from the relation in Equation (4), we use the estimated stable trait
factor variances to estimate Ψ as
Ψ̂ = S − Φ̂ ⊗ 1T 1tT ,where Φ̂ is an estimator of a 3 × 3 stable trait factor variance-covariance matrix Φ:2ϕ̂(Y,X) ϕ̂(Y,L)ϕ̂
 (Y )Φ̂ =  ϕ̂(Y,X) ϕ̂2(X) ϕ̂(X,L) 2ϕ̂(Y,L) ϕ̂(X,L) ϕ̂(L) .

(19)

(20)

Since stable trait factor covariances are not estimated in the previous sub-step, covariances
between calculated linear correlation preserving predictors from each variable are used (seeUsami (2022a) for further details). From Equations (16) and (18)–(20), we can calculate
WPVS Ẑi∗ without specifying the structural models among variables at the within-personlevel, successfully maintaining independence from the next step.

16

3.2

Applying SNMMs

SNMMs were developed in epidemiology to estimate the causal eﬀects of a sequence of timevarying treatments/predictors in longitudinal studies, with eﬀectively handling violation ofno unobserved confounders (U ) under the assumed causal DAG where U influences boththe outcomes Y and the observed confounders L (Robins & Hernán, 2009; Vansteelandt &Joﬀe, 2014). For notational simplicity, we use the symbols Y ∗ , X ∗ , and L∗ below to denote
WPVS that are calculated in the first step.
For explanation purposes, we assume that one is interested in evaluating the joint
eﬀects of X̄T∗ −1 on the outcome of the last time point (YT∗ ), while allowing its curvilinear(e.g., quadratic) and interaction eﬀects. SNMMs simulate the sequential removal of theeﬀect (called blip) that X̄T∗ −1 has on YT∗ , after having removed the eﬀects of all subsequent
treatments/predictors. More specifically, SNMMs models a blip in X̄t∗ on YT∗ while holding
all future treatments/predictors at t′ ≥ t fixed at a reference level 0 (i.e., the level that isequal to the expected scores of a person in the current context).
∗

Linear SNMMs parameterize contrasts of YiT∗ x̄i(t−1) and YiT∗ x̄i(t−2) ,0 conditionally on confounder histories through t − 1 (t = 2, . . . , T ) as
∗
∗
∗
∗
∗
∗
E(YiT∗ x̄i(t−1) − YiT∗ x̄i(t−2) ,0 |Ȳi(t−1)= ȳi(t−1), L̄∗i(t−1) = ¯li(t−1), X̄i(t−1)= x̄∗i(t−1) )∗∗= ht (ȳi(t−1), ¯li(t−1), x̄∗i(t−1) ; τ ),

(21)

∗
∗
where ht−1 (ȳi(t−1), ¯li(t−1), x̄∗i(t−1) ; τ ) is a known function with parameter vector τ (Vanstee-landt & Joﬀe, 2014).
Suppose we assume a DGP similar to the one in Figure 1 that has not only the firstorder (linear) eﬀect of X ∗ but also its interaction eﬀect with L∗ for blip. In the empirical17

example using the data for T = 4, a linear SNMM may be given by
∗

∗

∗

∗

∗

∗
∗
∗
∗
E(Yi4∗ xi1 ,xi2 ,xi3 − Yi4∗ xi1 ,xi2 ,0 |Ȳi3∗ = ȳi3
, L̄∗i3 = ¯li3
, X̄i3= x̄∗i3 ) = (β3 + γ3 li3)x∗i3 ,
∗
∗
∗
∗
∗
∗
∗
E(Yi4∗ xi1 ,xi2 ,0 − Yi4∗ xi1 ,0,0 |Ȳi2∗ = ȳi2
, L̄∗i2 = ¯li2
, X̄i2= x̄∗i2 ) = (β2 + γ2 li2)x∗i2 ,∗∗∗∗∗)x∗i1 .
E(Yi4∗ xi1 ,0,0 − Yi4∗ 0,0,0 |Yi1∗ = yi1, L∗i1 = li1, Xi1= x∗i1 ) = (β1 + γ1 li1(22)∗Here, the first equation models the eﬀect of Xi3
on Yi4∗ , the second models the eﬀect of∗∗Xi2
on Yi4∗ , and the third models the eﬀect of Xi1
on Yi4∗ . This is just one simple example,∗∗∗
but another specification, say (β3 + γ31 li1+ γ32 li2
+ γ33 li3)x∗i3 , could be possible to express∗
the interaction eﬀects between Xi3
and the past confounders. From Equation (22), the∗(conditional) average joint eﬀects of X̄i3
on Yi4∗ when increasing one unit from the reference

values in each predictor become∗∗∗
β3 + γ3 li3+ β2 + γ2 li2+ β1 + γ1 li1.

(23)

If no interaction eﬀects are assumed, this becomes
β41 + β42 + β43 .

(24)∗(τ ) (t = 2, . . . , T ), the mean value of which is
SNMMs consider a transformation Ui(t−1)

equal to the mean that would be observed if treatments/predictors were stopped from timet − 1 onward, in the sense that∗∗∗∗
E(Ui(t−1)(τ )|Ȳi(t−1), L̄∗i(t−1) , X̄i(t−2)= x̄∗i(t−2) , Xi(t−1)
)
∗

∗
∗
∗
=E(YiT∗ x̄i(t−2) ,0 |Ȳi(t−1), L̄∗i(t−1) , X̄i(t−2)= x̄∗i(t−2) , Xi(t−1)).
∗
Here, Ui(t−1)(τ ) is a vector with components YiT∗ −

(25)

∑T −1

∗
∗
∗
k=t−1 hk (Ȳik , L̄ik , Āik ; τ ). For instance,in the current example for T = 4,∗
Ui3∗ (τ ) = Yi4∗ − (β3 + γ3 L∗i3 )Xi3,∗∗
Ui2∗ (τ ) = Yi4∗ − (β3 + γ3 L∗i3 )Xi3
− (β2 + γ2 L∗i2 )Xi2
,
∗
∗
∗
.
− (β1 + γ1 L∗i1 )Xi1
− (β2 + γ2 L∗i2 )Xi2
Ui1∗ (τ ) = Yi4∗ − (β3 + γ3 L∗i3 )Xi3

18

(26)

The assumption of no unobserved confounders (Equation (7)) together with the identity(Equation (25)) implies that
∗
∗
∗
∗
∗
∗
E(Ui(t−1)(τ )|Ȳi(t−1), L̄∗i(t−1) , X̄i(t−1)) = E(Ui(t−1)(τ ∗ )|Ȳi(t−1), L̄∗i(t−1) , X̄i(t−2)
),

(27)

∗
∗
indicating the conditional independence between Xi(t−1)
and Ui(t−1). In the current example

for T = 4, this indicates∗Ui1∗ ⊥⊥ Xi1|Yi1∗ , L∗i1 ,∗∗,
Ui2∗ ⊥⊥ Xi2|Ȳi2∗ , L̄∗i2 , Xi1∗∗.
Ui3∗ ⊥⊥ Xi3|Ȳi3∗ , L̄∗i3 , X̄i2

(28)

The parameters τ can be estimated from these moment conditions by solving the estimating
equation EN [f (τ ; η̂; κ̂)] = 0, where EN denotes the empirical average function:f (τ ; η̂; κ̂) =

T
∑

∗
∗
∗
∗
∗
∗
[dt−1 (Ȳi(t−1), L̄∗i(t−1) , X̄i(t−1)) − E(dt−1 (Ȳi(t−1), L̄∗i(t−1) , X̄i(t−1))|Ȳi(t−1), L̄∗i(t−1) , X̄i(t−2)
)]◦

t=2
∗
∗
∗
∗
[Ui(t−1)(τ ) − E(Ui(t−1)(τ )|Ȳi(t−1), L̄∗i(t−1) , X̄i(t−2)
)],

(29)

where
[
∗
∗
dt−1 (Ȳi(t−1), L̄∗i(t−1) , X̄i(t−1)
)=E

∗
(τ )
∂Ui(t−1)

∂τt−1

]
∗
∗
Ȳi(t−1), L̄∗i(t−1) , X̄i(t−1)

.

(30)

Here, τt−1 denotes the elements in τ that are relevant to the assumption of conditional
independence (Equation (28)) at t − 1. For instance, in the current example for T = 4, we
can apply Equation (26) to obtain
] ( ∗
)t
[
∂Ui3 (τ ) ∂Ui3∗ (τ )∂Ui3∗ (τ ) ∗ ∗∗∗∗∗Ȳi3 , L̄i3 , X̄i3 =,d3 (Ȳi3 , L̄i3 , X̄i3 ) = E
∂τ3
∂β3
∂γ3
(
)t
∂
∂
∗
∗
∗
∗
∗
∗
=
[Y − (β3 + γ3 Li3 )Xi3 ] ,[Y − (β3 + γ3 Li3 )Xi3 ]∂β3 i4∂γ3 i4∗∗ t= (−Xi3, −L∗i3 Xi3))t( ∗∂Ui2 (τ ) ∂Ui2∗ (τ )
∗
∗ t
∗
∗
∗
,
= (−Xi2, −L∗i2 Xi2)d2 (Ȳi2 , L̄i2 , X̄i2 ) =
∂β2
∂γ2
( ∗
)t
∂Ui1 (τ ) ∂Ui1∗ (τ )
∗
∗
∗
∗
∗ t
d1 (Yi1 , Li1 , Xi1 ) =,= (−Xi1, −L∗i1 Xi1
).
∂β1
∂γ1

19

(31)

This estimating equation essentially sets the sum across the time points of the condi∗∗∗
tional covariances between Ui(t−1)(τ ) and the function dt−1 (Ȳi(t−1), L̄∗i(t−1) , X̄i(t−1)), given

confounders, to zero. More specifically, in this example, the elements in Equation (29)essentially become
∗
∗
∗
∗
)]
)] × [Ui3∗ (τ ) − E(Ui3∗ (τ )|Ȳi3∗ , L̄∗i3 , X̄i2
[Xi3
− E(Xi3|Ȳi3∗ , L̄∗i3 , X̄i2
∗
∗
∗
∗
∗
[L∗i3 Xi3
− Li3
E(Xi3
|Ȳi3∗ , L̄∗i3 , X̄i2)] × [Ui3∗ (τ ) − E(Ui3∗ (τ )|Ȳi3∗ , L̄∗i3 , X̄i2
)]
∗
∗
∗
∗
[Xi2
− E(Xi2|Ȳi2∗ , L̄∗i2 , Xi1)] × [Ui2∗ (τ ) − E(Ui2∗ (τ )|Ȳi2∗ , L̄∗i2 , Xi1
)]
∗
∗
∗
∗
)]
)] × [Ui2∗ (τ ) − E(Ui2∗ (τ )|Ȳi2∗ , L̄∗i2 , Xi1[L∗i2 Xi2
− L∗i2 E(Xi2|Ȳi2∗ , L̄∗i2 , Xi1∗∗− E(Xi1|L∗i1 , L∗i1 )] × [Ui1∗ (τ ) − E(Ui1∗ (τ )|Yi1∗ , L∗i1 )][Xi1∗∗[L∗i1 Xi1
− L∗i1 E(Xi1|Yi1∗ , L∗i1 )] × [Ui1∗ (τ ) − E(Ui1∗ (τ )|Yi1∗ , L∗i1 )].

(32)∗
If homoscedasticity of the conditional variance Ui(t−1)(τ ∗ ) given confounders is satisfied, the∗∗),, L̄∗i(t−1) , X̄i(t−1)
local semiparametric eﬃciency under the SNMM is attained by choosing dt−1 (Ȳi(t−1)

as in Equation (30) (Vansteelandt & Joﬀe, 2014).
As seen in Equation (32), solving the estimating equation (29) requires a treatment∗∗∗∗; η). It also
, L̄∗i(t−1) , X̄i(t−2)|Ȳi(t−1): f (Xi(t−1)
model A for the treatment/predictor Xi(t−1)

requires a treatment-free model B (Wallace et al., 2017b) for the conditional mean of∗∗∗∗; κ). Notably, when the parameters, L̄∗i(t−1) , X̄i(t−2)(τ )|Ȳi(t−1)(τ ), namely, f (Ui(t−1)
Ui(t−1)η and κ are variation-independent, the estimator that solves Equation (29) (called Gestimator), obtained by substituting η and κ with consistent estimators, are doubly robust(Robins & Rotnitzky, 2001, cited from Vansteelandt & Joﬀe, 2014), meaning that estimatesof causal parameters are consistent when either model A or model B is correctly specified.
∗Though we have supposed situations with joint eﬀects of X̄3∗ on Y4∗ , the eﬀects of X̄t−1

on Yt∗ for t < 4 can also be modeled and estimated in a similar manner. Likewise, the∗opposite relation for the eﬀects of Ȳt−1
on Xt∗ can also be investigated with a separate

procedure.
To summarize, in the second step of using (linear) SNMMs, we first specify the blip

20

model (e.g., Equation (22)) for time-varying treatments/predictors that may have curvilin∗2∗
ear and interaction eﬀects (e.g., Xt−1
and Xt−1
L∗t−1 ). We then consider a transformation∗
Ui(t−1)(τ ) (e.g., Equation (26)), and the parameters for blip (τ ) are estimated based on

moment conditions implied by conditional independence (e.g., Equation (28)). When estimating the parameters, the treatment model A and the treatment-free model B need to
be specified based on the assumed DGP, but SNMMs have the property of being doubly
robust. We further explain below the specifications of these models using the R package
DTRreg (Wallace et al., 2017b).

3.3

SNMMs via the R package DTRreg

The framework of SNMMs is usually presented using estimation equations like Equation (29). However, Wallace et al (2017a) explained that the same calculation may infact be conducted using a relatively straightforward series of matrix equations (see Wallace & Moodie, 2015, and the Web Appendix2 provided by Wallace et al., 2017a for moredetails). More specifically, the final stage (i.e., t = T ) blip parameters are estimated first,
before working backwards until every stage at t has been analyzed. By working recursively,
we are able to calculate each potential outcome by plugging in all future blip parameters(the Web Appendix in Wallace et al., 2017a).
Suppose we assume a DGP similar to that in Figure 1, and first-order (linear) eﬀects
of X ∗ and its interaction eﬀect with L∗ are assumed for the blip model in T = 4, as in∗Equation (22). β3 and γ3 (i.e., the eﬀects associated with Xi3
on Yi4∗ ) are first estimated∗and Ui3∗
using the moment condition induced by the conditional independence between Xi3(Equation (28)). After plugging in these estimates, β2 and γ2 are estimated next based on∗
the conditional independence between Xi2and Ui2∗ . Finally, using the estimates obtained∗and Ui1∗ .
so far, β1 and γ1 are estimated based on the conditional independence between Xi1

2

http://links.lww.com/EDE/B134

21

Even if we take this matrix equations approach, specifications of the treatment model A
and the treatment-free model B are required. These models are expressed in the R package
DTRreg as treat.mod and tf.mod, respectively.
In the current example, from the assumed DGP A can be specified as an (linear) AR(1)model. If a researcher is interested in the joint eﬀect of X̄3∗ on Y4∗ , treat.mod may be
specified in the R package DTRreg as
treat.mod <- list(X1~1, X2~X1+L1+Y1, X3~X2+L2+Y2).
Note that X1~1 indicates that only an intercept is included because X1∗ is now treated as
an exogenous variable. In B, because the outcomes Y ∗ are modeled by the similar (linear)AR(1) model (with interaction eﬀects between X ∗ and L∗ ) in the assumed DGP, if X3∗ = 0,then the conditional means of Y4∗ used for U3 (τ ) (given Ȳ3∗ , L̄∗3 , and X̄2∗ ) can be expressedby a linear (weighted) sum of Y3∗ and L∗3 . Next, if X2∗ = X3∗ = 0, then the conditional
means of Y4∗ used for U2 (τ ) (given Ȳ2∗ , L̄∗2 , and X1∗ ) becomes a linear (weighted) sum of Y2∗
and L∗2 . Likewise, if X̄3∗ = 0, then the conditional means of Y4∗ used for U1 (τ ) (given Y1∗
and L∗1 ) becomes a linear (weighted) sum of Y1∗ and L∗1 . Thus, tf.mod can be specified in
the R package DTRreg as
tf.mod <- list(~Y1+L1, ~Y2+L2, ~Y3+L3).
Note that one can allow curvilinear and interaction eﬀects for WPVS in A and B, but thiscan be challenging when using RI-CLPM as covariance structure analysis.
In addition to these, the blip model needs to be specified in the R package DTRreg.
Since the interaction eﬀect between X ∗ and L∗ is now assumed, the blip model (expressedas blip.mod) can be expressed as
blip.mod <- list(~L1, ~L2, ~L3).
If no interaction eﬀect is assumed, then the specification becomes
blip.mod <- list(~1, ~1, ~1).

22

The causal parameters τ can now be estimated in the R package DTRreg by three models
specified so far:mod<-DTRreg(Y4, blip.mod, treat.mod, tf.mod, var.est="sandwich", type="alt"),
where the option var.est="sandwich" specifies the robust (or sandwich) variance estimator obtained using standard estimating equation theory (Web Appendix in Wallace et al.,2017a; see Hardin & Hilbe, 2013 pp.30-34 for further details). Though nuisance parametersη and κ need to be accounted for to obtain better estimates of standard errors, Wallace
et al. (2017a) explained that sandwich estimators that ignore the nuisance parameter estimation typically perform as well as the bootstrap or the nuisance parameter corrected
standard error. The package DTRreg also oﬀers a variety of more complex options, including bootstrap, and it can be specified by a combination of var.est="boot" and B=n for
the number of bootstrap replications (Wallace et al., 2017b).
Using the command type="alt", summary(mod) and coef(mod) return summaries andblip parameter estimates in a fashion similar to more familiar commands like lm and glm(Web Appendix in Wallace et al., 2017a). Additionally, the package DTRreg will automatically ignore any persons with missing data (thereby carrying out a complete-cases analysis),but if the option missing = "ipcw" is specified, then the inverse probability of censoredweights (e.g., Hernan & Robins, 2021) is used. The probability of censoring is estimated via
logistic regression on the full covariate history up to that point (Web Appendix in Wallace
et al., 2017a).

4

SIMULATIONSThere are two main goals in the simulations. The first is to investigate the performanceof TS-SNMMs under various data-generating conditions in which data are generated by alinear sum of stable trait factors and WPVS, as in RI-CLPM. In the first simulation, we

23

assume there are no misspecifications in the structural models. The second is to demonstrate the robustness of TS-SNMMs and to compare its performance to RI-CLPM under
the presence of unobserved time-varying confounders U ∗ (that influence outcomes Y ∗ andobserved time-varying confounders L∗ ) and model misspecifications caused by observed
time-varying confounders L∗ (whose direct second-order eﬀect on the outcome is missing)in the structural model. Second and third simulations are performed for this purpose.
Throughout the simulation, we suppose a situation where researchers want to evaluate∗∗∗
the joint eﬀects of Xi(T−2) and Xi(T −1) (we do not focus on the eﬀects from X̄i(T −3) ) on

YiT∗ , and vice versa for the reciprocal relation. For simplicity, we also assume that theinteraction eﬀects of predictors with observed confounders are not present.

4.1

Scenario 1: No misspecifications in the structural models

4.1.1

MethodThree stable trait factors are first generated by the multivariate normal:


(Y )



 


0

Ii

ϕ2

ϕ2 rB ϕ2 rB     2
 (X) 2
2
 Ii  = M V N  0  ,  ϕ rB ϕ
ϕ rB   ,  


(L)
2
2
2
ϕ rB ϕ rB ϕ
0
Ii

(33)

where ϕ2 is a stable trait factor variance and rB is a correlation between factors. The initial
WPVS are generated independently from the stable trait factors as


Yi1∗





 



0
1−ϕ
(1 − ϕ )rW (1 − ϕ )rW     ∗ 
 Xi1  = M V N  0  ,  (1 − ϕ2 )rW
1 − ϕ2(1 − ϕ2 )rW  ,  222∗(1 − ϕ )rW (1 − ϕ )rW
1−ϕ
0
Li1
2

2

2

(34)

where rW is a correlation among the initial WPVS. Then, the WPVS at successive times
are sequentially generated via a first-order linear regression model with time-invariant co-24

eﬃcients:∗∗
Yit∗ = 0.30Yi(t−1)+ 0.20Xi(t−1)+ 0.20L∗i(t−1) + dit

(Y )

∗
∗
Xit∗ = 0.20Yi(t−1)+ 0.30Xi(t−1)+ 0.20L∗i(t−1) + dit

(X)

(35)

∗
∗
L∗it = 0.20Yi(t−1)+ 0.20Xi(t−1)+ 0.30L∗i(t−1) + dit .
(L)

When T = 3, this parameter setting indicates (see Equation (10) for the calculation)∗

E(Yi3∗ x̄i2 ) = E(0.30Yi2∗ + 0.20L∗i2 ) + 0.20x∗i2 = (0.06x∗i1 + 0.04x∗i1 ) + 0.20x∗i2 = 0.10x∗i1 + 0.20x∗i2∗∗ ȳi2∗∗E(Xi3) = 0.10yi1+ 0.20yi2
.

(36)

∗
∗
∗
More generally, the joint eﬀects of Xi(T−2) and Xi(T −1) (with X̄i(T −3) not manipulated)∗∗on YiT∗ can be evaluated as E(YiT∗ xi(T −2) ,xi(T −1) ) = 0.10x∗i(T −2) + 0.20x∗i(T −1) . Similarly,∗∗∗∗ yi(T −2) ,yi(T −1)∗
E(XiT) = 0.10yi(T−2) + 0.20yi(T −1) . Therefore, accurately estimating these fourcoeﬃcients (0.10, 0.20, 0.10, 0.20)t for within-person reciprocal relations is a shared goal
between TS-SNMMs and RI-CLPM.
The variance of the normal residual d for each variable was set so that ratio of thevariance of WPVS at time t (e.g., var(Yit∗ ) for Y ) to that at time t = 1 (e.g., var(Yi1∗ ))becomes var(Yit∗ )/var(Yi1∗ ) = T −t+(t−1)ω
. More specifically, when t = T , this ratio becomesT −1
ω (i.e., ω = var(YiT∗ )/var(Yi1∗ )). When ω = 1, the variance of WPVS is constant over time.
Note that residual variances were controlled over time in the simulations of Usami (2022a).
Like RI-CLPM, measurements are then generated as(Y )Yit = Ii

+ Yit∗ , Xit = Ii

(X)

+ Xit∗ , Lit = Ii

(L)

+ L∗it ,(37)

where temporal group means are set to zero at each time point for each variable.
In this simulation, we systematically changed the total number of persons to N =200, 600, 1000, the number of time points to T = 3, 4, 6, 9, the variances of stable trait
factors to ϕ2 = 0.1, 0.4, 0.7, the correlations between stable trait factors to rB = 0.1, 0.3, 0.5,
the correlations between initial WPVS to rw = 0.1, 0.3, 0.5, and the ratio of the variance25

of WPVS at time t = T to that at time t = 1 to ω = 1, 3, 5. Some major diﬀerencesfrom the simulations in Usami (2022a) are that more varied specifications of T and themanipulations of rB , rW , and ω. By crossing these factors, we generated 200 simulation
data for each combination of factors. Under each simulation condition, we calculated thebias and root mean squared error (RMSE) of estimates from TS-SNMMs and RI-CLPM(which assumes first-order lagged regressions that include L∗ ).
In the first step of TS-SNMMs, a model that assumes a linear AR(1) structure with
time-varying autoregressive coeﬃcients and residual variances was specified for the WPVS
in each variable when T = 3. The models with similar AR(2) and AR(3) structures werealso specified in the T = 4 and T = 6, 9 conditions, respectively. Note that in the current
DGP, a true measurement model for each variable has an AR(T − 1) structure in WPVS,
but such a measurement model cannot be identified. Therefore, the measurement model
is more or less misspecified in the current setting. In the second step, both the treatment
model A and the treatment-free models B were correctly specified. For example, in T = 9(the joint eﬀects of X7∗ and X8∗ on Y9∗ ), three models are specified in the package DTRreg
astreat.mod <- list(X7~X6+L6+Y6, X8~X7+L7+Y7)
tf.mod <- list(~Y7+L7, ~Y8+L8)
blip.mod <- list(~1, ~1).
Estimation results are discarded when improper solutions because of out-of-range parameter estimates (e.g., negative variance or the singular Hessian matrix for trait factors)
occur when applying either TS-SNMMs or RI-CLPM. The simulation was conducted in R
using the package lavaan (Rosseel, 2012) to estimate parameters with MLE in RI-CLPM
and in measurement models for the first step of TS-SNMMs. The package DTRreg wasused in the second step of TS-SNMMs. The simulation code is available in the Online
Supplemental Material.

26

4.1.2

Results

Figure 2 shows biases of causal parameter estimates from the TS-SNMMs (with an AR(1)measurement model in the first step) and RI-CLPM. Because similar tendencies were observed, we provide results for the eﬀects of X ∗ on Y ∗ (the sizes of the first-order eﬀect
from XT∗ −1 and the second-order eﬀect from XT∗ −2 are 0.20 and 0.10, respectively), andthose of Y ∗ on X ∗ are provided in the Online Supplemental Material (Figures S10-S11).
Also, the impacts of N and rW were not relatively large in terms of biases, especially inthe second-order eﬀect estimates, which exhibited a larger variability than the first-ordereﬀect estimates (see Table S1 of the ANOVA results, where calculated biases are set as
outcomes). Here, we provide results from N = 1000 and rW = 0.3, and we confirm thatwhole conclusions are not influenced by diﬀerences in N and rW (Figures S2-S5).
Since the data were generated by the process assumed in RI-CLPM, RI-CLPM canaccurately estimate the parameters. However, since improper solutions often arise, especially when N is small or ϕ2 = 0.1 (i.e., small true stable trait variances cause negative
estimates), small (negative) biases occur in RI-CLPM as a result of discarding estimateswith improper solutions. The magnitudes of RMSEs were similar between TS-SNMMs and
RI-CLPM (Figures S1,S6-S9), so we will focus particularly on the results of biases below.
TS-SNMMs show biases in some conditions because of the misspecified measurement
model in the first step. However, these biases become smaller as T becomes larger andϕ2 becomes smaller. When T is as small as T = 3, TS-SNMMs show large biases inmany conditions (e.g., the relative biases over true eﬀects exceed 10%), so even though TSSNMMs can be used in T = 3 for some conditions, such as ϕ2 = 0.1, it is generally advisable
not to do so. On the other hand, when T ≥ 4, biases were within 10% of the relative biases
over true eﬀects in most conditions, the exceptions being when trait factor correlations
were small and the proportions of variances explained by within-person fluctuations (i.e.,
WPVS) were kept small over time (i.e., large ϕ2 = 0.4, 0.7 and small rB = 0.1 and ω = 1).

27

The overall results become more accurate when T ≥ 6.
When an AR(2) measurement model was used in TS-SNMMs for T = 4, the positive
biases observed in some conditions were mitigated, but the overall diﬀerences were almost
ignorable (see Figures S12-S13 in the Online Supplemental Material). The same tendencies
were observed when either AR(2) or AR(3) measurement models were used in T ≥ 6.
Therefore, the choice of AR order in the measurement model does not have a large impact
on biases, at least in the current parameters setting.
Table 1 provides the proportions of improper solutions calculated when a total of 200
proper solutions were obtained in each condition. It shows that improper solutions very
often arise in RI-CLPM, especially when ϕ2 is small, but also that the sizes of both N and
T have large impacts. RI-CLPM produces improper solutions more frequently than TSSNMMs in all conditions. This diﬀerence results from the fact that TS-SNMMs calculate
WPVS from the measurement model of each variable, while RI-CLPM estimates structuralmodels and measurement models for all variables simultaneously. With the exception of
ϕ2 = 0.1, improper solutions arise in less than 1% of the simulations in TS-SNMMs whenT ≥ 4. When T = 4, 6 and ϕ2 =0.1, improper solutions arise in more than 15% in allconditions, even if N becomes larger in RI-CLPM. In cases of small N , the risk of improper
solutions in RI-CLPM becomes an issue if either T is small (T ≤ 4) or ϕ2 is small (ϕ2 =0.1).

4.2

Scenario 2: Unobserved confounders that influence observed
time-varying confounders and outcomes

4.2.1

MethodWe consider the DGP in which an unobserved time-varying confounder U ∗ influences boththe observed time-varying confounders L∗ and the outcomes Y ∗ . In this scenario, controllingfor L∗ (as a collider) to infer the causal eﬀects of X ∗ on Y ∗ introduces bias through U ∗ on
the path from X ∗ to Y ∗ (i.e., XT∗ −2 → L∗T −1 ← UT∗ −2 → UT∗ −1 → YT∗ ). This is known as

28

collider bias, and this DGP is in line with the one typically assumed in the epidemiologicalliterature on MSMs and SNMMs (e.g., Robins & Hernán, 2009; Vansteelandt & Joﬀe,2014). Such a DGP might be assumed when researchers perform a sequential randomized
trial and when the treatment model can be correctly specified (i.e., one is sure how thetreatments/predictors X ∗ can be explained by other observed variables, but the outcomesY ∗ and the observed confounders L∗ may be influenced by unobserved confounders U ∗ ).
∗For this DGP, the AR(1) process is assumed in the model for U ∗ as Uit∗ = 0.7Ui(t−1)+dit ,(U )and the variance of the residual dit is manipulated according to ω like other variables. Y ∗(U )and L∗ are then generated as∗∗∗Yit∗ = 0.30Yi(t−1)+ 0.20Xi(t−1)+ 0.20L∗i(t−1) + 0.30Ui(t−1)
+ dit

(Y )

∗
∗
∗
L∗it = 0.20Yi(t−1)+ 0.20Xi(t−1)+ 0.30L∗i(t−1) + 0.30Ui(t−1)+ dit ,
(L)

(38)

while X ∗ is generated as in the previous simulation. Both the treatment model A and thetreatment-free model B were specified as in the previous simulation, and we compared theestimation performances of TS-SNMMs and RI-CLPM that lacked U ∗ . Namely, in TSSNMMs B is misspecified in inferring eﬀects from X ∗ on Y ∗ , whereas A is still correctly
specified. Thus, TS-SNMMs can take doubly robust property.
4.2.2

Results

Figure 3 shows the biases of estimates for the eﬀects of X ∗ on Y ∗ in TS-SNMMs (thatassume an AR(1) measurement model in the first step) and RI-CLPM. Because the impactsof N and rW on the results were again smaller in terms of biases, we provide the resultsfor N = 1000 and rW = 0.3. We also confirm again that the choice of AR order inthe measurement model does not have a large impact on biases in TS-SNMMs (FigureS14). From Figure 3, we can see that TS-SNMMs show a similar amount of biases in eachcondition to those in the previous simulation. Namely, with the exception of the specificconditions shown above (large ϕ2 = 0.4, 0.7 and small rB = 0.1 and ω = 1), the biases were29

small in T = 4, and the overall results become more accurate as T became larger.
On the other hand, notable positive biases were observed in RI-CLPM for the secondorder eﬀect of X ∗ . With the exception of the conditions specified above, the biases became
larger in RI-CLPM than those in TS-SNMMs in T ≥ 6, and RI-CLPM exhibited morethan a 10% level of relative biases over the true eﬀects in all conditions. This demonstrates
that estimates from RI-CLPM can be seriously biased by unobserved confounders (thatinfluences both observed confounders and outcomes), while TS-SNMMs are robust to thisinfluence. Note that because the first-order eﬀect is not influenced by collider bias in
the current DGP, except for the conditions specified above, estimation performances arecompetitive between methods.
These methods are also competitive in estimates for the second-order eﬀects of Y ∗ onX ∗ , but they commonly show small amounts of biases (Figure S15) because in TS-SNMMs
both A and B are misspecified in the current DGP. More specifically, the model for YT∗−2 isinfluenced by the previous U ∗ , and so too is XT∗ via the causal path UT∗ −2 → L∗T −1 → XT∗ .
Therefore, in the current DGP, TS-SNMMs can more accurately estimate the eﬀects of X ∗on Y ∗ than RI-CLPM, but not the eﬀects of Y ∗ on X ∗ .
As shown in Table S2, similar tendencies were observed for the frequency of improper
solutions as in the previous simulation. RI-CLPM very often produces improper solutions:at least 30% of all cases when ϕ2 is as small as 0.1. Even for larger ϕ2 , when T and N
are as small as T = 4, 6 and N = 200, the observed proportions were around 10%. On theother hand, the risk of obtaining improper solutions in TS-SNMMs was much smaller than
RI-CLPM.
One might consider that when such unobserved confounders exist, one can avoid referring to such biased estimation results from RI-CLPM by diagnosing some model fit indices.
However, in the current scenario, 98.06% of the estimation results from RI-CLPM showed
an SRMR lower than .05. Likewise, 99.97% and 97.23% of the results from RI-CLPM
showed a CFI larger than .95 and an RMSEA lower than .05, respectively (see Table S3).

30

These proportions exceeded 70–80% on average, even when T and N were small. Therefore,
it would be very diﬃcult to detect collider biases caused by U ∗ through these model fit
indices alone.

4.3

Scenario 3: Ignored direct higher-order eﬀects of observed
time-varying confounders on outcomes

4.3.1

MethodIn this scenario, we assume that there are ignored direct higher-order eﬀects of L∗ on Y ∗ .
Although a first-order lagged eﬀect is often assumed when applying the (RI-)CLPM, anignored (direct) higher-order eﬀect leads to biased causal estimates unless it is zero, andassessing the number of order is a major task in inferring within-person relations. However,as the numbers of time points and observed time-varying confounders increase, this task
becomes more complex, and the risk of misspecifications can increase.
In this scenario, Y ∗ is generated as∗∗Yit∗ = 0.30Yi(t−1)+ 0.20Xi(t−1)+ 0.20L∗i(t−1) + 0.20L∗i(t−2) + dit ,
(Y )

(39)

while X ∗ and L∗ are generated similarly as in the first simulation. Because L∗i(T −2) influ∗∗
ences both Xi(T−1) and YiT , it is obviously a confounder in evaluating the relation between
∗
∗
∗
∗
X(T
−1) and YT , and we expect ignoring Li(T −2) in structural models for YT leads to biased

causal estimates in RI-CLPM. Both the A and B models were specified as in the firstsimulation, and we compared the estimation results between TS-SNMMs and RI-CLPM
that did not include L∗i(t−2) . Namely, B was misspecified, but A was still correctly specified
when estimating the eﬀects of X ∗ on Y ∗ in TS-SNMMs.

31

4.3.2

Results

Figure 4 shows the biases of causal eﬀect estimates of X ∗ on Y ∗ in TS-SNMMs (thatassume an AR(1) measurement model in the first step) and RI-CLPM. Since the impactsof N and rW on the results were again smaller in terms of biases, we provide here theresults when N = 1000 and rW = 0.3. We also confirmed again that the choice of AR
order in the measurement model does not have a large impact on biases in TS-SNMMs(Figure S16). From Figure 4, we can see that TS-SNMMs show a similar amount of biases
in each condition as the first simulation (i.e., Figure 2). Therefore, with the exception of
the conditions mentioned above, the biases were small in T ≥ 4, but the overall results
became more accurate when T ≥ 6.
Notably, both the first- and second-order eﬀects have large positive biases (that exceed a
10% level of relative biases over true eﬀects) in RI-CLPM, and the amounts of biases in RICLPM exceed those in TS-SNMMs in most conditions when T ≥ 4. This clearly indicatesthat causal estimates from RI-CLPM can be seriously biased by model misspecification in
the form of omitted direct higher-order lagged eﬀects of L∗ in the structural models, even if
there are no other variables as unobserved confounders U ∗ in the DGP. On the other hand,estimates for the opposite relation (Y ∗ on X ∗ ) were similar to those in Figure 2 for bothmethods due to there being no ignored direct higher-order eﬀects in the structural modelof X ∗ (Figure S17).
As shown in Table S4, similar tendencies were observed for the frequency of improper
solutions as in the previous simulations. RI-CLPM often produces improper solutions when
ϕ2 is as small as 0.1. Even for larger ϕ2 , when T and N are small as T = 4 and N = 200,
the observed proportions were more than 10%.
Like the previous simulation, we evaluated how model fit indices work to avoid referring
to estimation results from RI-CLPM. However, 94.63% of the estimation results from RICLPM (that ignore direct second-order lagged eﬀects) showed an SRMR lower than .05.

32

Likewise, 99.09% and 87.06% showed a CFI larger than .95 and an RMSEA lower than .05,
respectively (see Table S5). These proportions exceeded 30–50% on average, even when
T and N were small. Therefore, it seems very diﬃcult to detect estimation biases causedby misspecification in the form of omitted direct higher-order lagged eﬀects of L∗ in the
structural model by using these model fit indices.
We have demonstrated that TS-SNMMs can be eﬀectively used in many conditions
where longitudinal data with T ≥ 4 are available, and that more accurate causal estimates
can be obtained under some situations even if unobserved time-varying confounders andmodel misspecifications exist, with a lower risk of obtaining improper solutions comparedwith RI-CLPM.

5EMPIRICAL APPLICATIONIn this section, we describe an empirical application of TS-SNMMs using data from theTokyo Teen Cohort (TTC) study (Ando et al., 2019), which was also used in Usami (2022a)
with T = 3 waves. We assume a similar causal DAG model to that in Figure 1. Namely,
we assume that measurements are expressed by the linear sum of stable trait factors and
WPVS. TTC is a longitudinal cohort study for investigating the psychological and physicaldevelopment of (N =3,171) adolescents in the Tokyo metropolitan area, and data for whichhave been gathered in T = 4 waves: from 2012 to 2015 (age 10), from 2014 to 2017 (age12), from 2017 to 2019 (age 14), and from 2020 to 2022 (age 16).
The focus of this analysis is on the reciprocal relations between sleep duration (X ∗ )and depressive symptoms (Y ∗ ), and we estimate the eﬀects that past sleep duration (X ∗ )at ages 10, 12, and 14 has on later depressive symptoms (Y ∗ ) at age 16 (measured by the
Short Mood and Feelings Questionnaire; SMFQ, Angold et al., 1995), and vice versa, basedon the T = 4 waves. Several epidemiological studies have suggested a relationship between
sleep habits (sleep duration, bedtime, and bedtime regularity) and mental health status

33

(depression and anxiety) in adolescents, and their reciprocal relations were investigatedby Matamura et al. (2014). However, this study did not account for unit eﬀects (i.e.
stable traits) in sleep duration and depressive symptoms, and the relation at the withinperson level was also not investigated. In inferring the eﬀects of sleep duration on depressive
symptoms, the statistical control of other sleep habits like bedtime may be key. However, as
illustrated in the previous simulation, controlling for observed confounders L∗ can introduce
collider bias via unobserved confounders U ∗ that influence both L∗ and the outcome Y ∗ .
In this example, some life habits and the home environment (e.g., discipline from parents,
engagement in clubs/extracurricular activities in school) might aﬀect the level of depressive
symptoms and bedtime. Therefore, the use of TS-SNMMs can be considered a reasonable
way to robustly estimate the causal eﬀects of sleep duration over time.
The SMFQ consists of 13 items that assess depressive symptoms (0: not true, 1: sometimes true, 2: true) related to feelings and actions over the preceding two weeks. Higher
SMFQ scores suggest more severe symptoms. Sleep duration in hours was measured by
the question “How long do you usually sleep on weekdays?” Bedtime was used as observed
confounders (L∗ ), which was measured by the question “When do you usually go to bed on
weekdays?” In the present example, we focus on N = 1, 294 adolescents who consistently
responded to these items during four waves. Descriptive statistics of these variables areavailable in Table S6.
In the first step, we estimated the measurement model that assumes an AR(1) structurewith time-varying autoregressive coeﬃcients and residual variances in WPVS by MLE, and
WPVS were calculated for each variable. In the second step, we used the calculated WPVS
to estimate the joint eﬀect of sleep duration at ages 10, 12, and 14 on depressive symptoms
at age 16 (i.e., the joint eﬀect of X1∗ , X2∗ , X3∗ on Y4∗ ), and vice versa. Models A and B
were both specified based on first-order linear regression models (see the specifications in
the package DTRreg shown in Subsection 3.3). For the blip model, we specified one thatassumes interaction eﬀects between predictors and observed confounders (e.g., X3∗ L∗3 and34

X3∗ Y3∗ on Y4∗ ). For comparison, we also applied RI-CLPM assuming AR(1) regressions that
include L∗ (i.e., Equation (2)).
We confirmed that the first step did not produce improper solutions for each variable,
and that the specified AR(1) measurement model shows good fit in terms of the model fit
indices. Table S7 summarizes the model fit indices for each variable. Proportions of thevariances in measurements attributable to estimated stable trait factors at an initial time
(like the ϕ2 manipulated in the previous simulations) were calculated as 24.3%, 18.0%, and27.2% for depressive symptoms, sleep duration, and bedtime, respectively. The ratio of
variance in WPVS at time t = 4 to that at time t = 1 (ω) was 2.30, 1.45, and 2.22 fordepressive symptoms, sleep duration, and bedtime, respectively. The amounts of biases in
TS-SNMMs were pragmatically small in the previous simulations under the conditions forthe calculated ϕ2 and ω.
Table 2 gives the estimated joint eﬀect for reciprocal relations in TS-SNMMs (thatassume interaction eﬀects) and RI-CLPM, respectively. Bootstrap was used in RI-CLPM
to estimate standard errors.
As seen in Table 2, TS-SNMMs reveal that sleep duration at age 14 (X3∗ ) shows astatistically significant eﬀect for decreasing later depressive symptoms at age 16 (Y4∗ ) at
the within-person level (β̂y3 = −0.305, 95%CI [-0.595, -0.016], p <.05: 1 hour longer sleep
at age 14 decreases the SMFQ score by 0.305 points at age 16), while the main eﬀects
for age 12 (X2∗ ) and 10 (X1∗ ) do not. A statistically significant interaction eﬀect between
sleep duration and depressive symptoms at age 12 (X2∗ Y2∗ ) was observed (γ̂y2 = −0.173,
95%CI [-0.331, -0.016], p <.05: adolescents with higher depressive symptoms at age 12
take larger eﬀects on later depressive symptoms at age 14 from having more sleep), while
such interaction eﬀects were not observed at age 14 (X3∗ Y3∗ ). On the other hand, depressive
symptoms at age 12 (Y2∗ ) show a small, but statistically significant, main eﬀect at increasing
sleep durations at age 16 (X4∗ ) (β̂x2 = 0.018 [60× 0.018=1.08min increase], 95%CI [0.003,
0.032], p <.05), while the main eﬀects from ages 14 (Y3∗ ) and 10 (Y1∗ ) do not. On the35whole, the results were unchanged when using AR(2) measurement models for each variable(Table S8). Similar positive eﬀects of sleep duration on later depressive symptoms were
found in a previous study (Matamura et al., 2014). Usami (2022a), which mainly focusedon a clinical group comprising N = 416 adolescents with SMFQ scores of 6 or higher
during the study, also demonstrated that increased sleep duration has positive eﬀects onlater depressive symptoms.
RI-CLPM produced similar results to those from TS-SNMMs in that the eﬀects of sleep
duration at ages 12 (X2∗ ) and 10 (X1∗ ) were not significant, but sleep duration at 14 years old(X3∗ ) showed a significant eﬀect for decreasing later depressive symptoms (β̂4(Y )= −0.901,

95%CI [-1.533, -0.282], p <.05). Compared with TS-SNMMs, this point estimate is smaller
and indicates that a longer sleep duration has a more positive eﬀect on later depressive
symptoms. On the other hand, depressive symptoms did not exhibit any statisticallysignificant eﬀects on later sleep duration. The TTC study is now gathering new data for a
fifth wave, and further investigations with larger T are desired.
As illustrated in this example, the statistical significance, sign, and magnitude of estimates of causal eﬀects might change depending on how the eﬀects are modeled and how theparameters are estimated. TS-SNMMs that can flexibly and robustly estimate parameters
can be considered as potential alternatives of RI-CLPM (and its variants), and comparing
estimation results between these approaches could also be useful as a sensitivity analysis.

6

DISCUSSION

In this paper, we provided a tutorial, simulation, and illustrative example of TS-SNMMs
to infer within-person (reciprocal) relations. This method assumes a DGP similar to RICLPM, but we have shown through our simulations and the illustrative example that TSSNMMs have several advantages over RI-CLPM: (i) the flexible inclusion of curvilinear (e.g.,quadratic) and interaction eﬀects of WPVS as latent variables, (ii) more accurate estimates36

of causal parameters can be obtained under certain conditions due to them being doubly
robust, even if unobserved time-varying confounders and model misspecifications exist, and(iii) the risk of obtaining improper solutions is minimized. We showed in simulations that
TS-SNMMs work well in many conditions if longitudinal data with T ≥ 4 are available,
and more accurate estimates can be obtained if T becomes larger.
In psychology and related disciplines, the use of RI-CLPM has rapidly increased in
the last few years. Although researchers often wish to look at within-person (reciprocal)
relations from the perspective of causal inference, there are potential limitations in traditional SEM-based approaches like RI-CLPM when it comes to accounting for time-varying
confounders (which relate to the ability to flexibly and robustly estimate causal parameters), and the risk of obtaining improper solutions can present a challenge, especially when
the sample size is small. TS-SNMMs can be used as an alternative over RI-CLPM (andits variants), especially when data with T ≥ 4 are available and when improper solutions
occur in RI-CLPM. However, the purpose of this paper is not to completely denounce theuse of RI-CLPM. Comparing estimation results between these methods can be useful as asensitivity analysis, especially when either T is small or researchers are concerned with thepresence of unobserved time-varying confounders and model misspecifications.
There are other scenarios in which biased estimates can arise in RI-CLPM while TSSNMMs perform well. One is the case where some model misspecifications exist in the
structural model for the observed confounders L∗ . Even if researchers are sure that observed confounders were appropriately controlled for, the model for L∗ might be misspecified for some reasons (e.g., a misspecified functional form and/or omitted variables). In thestandard SEM approach as covariance structure analysis, one needs to specify the (linear)model for L∗ , even if the main focus is the inference of relations between the focal variables
X ∗ and Y ∗ . For this reason, model misspecification in L∗ can cause biases in the causal
parameters when applying RI-CLPM. Future investigations should study how a misspecification in L∗ can cause biased causal estimates in RI-CLPM depending on the DGP and37

model specifications, and how eﬀectively TS-SNMMs can mitigate this problem.
However, TS-SNMMs are also prone to possible pitfalls that practitioners should beaware of. The number of time points T is an especially critical aspect. Most research aimed
at inferring reciprocal (within-person) relations has used longitudinal data with T = 2, 3
(e.g., Usami, Todo et al., 2019), and imprecise predictions of WPVS can seriously degrade
the estimation performance in TS-SNMMs. However, in this case, the risk of producing
seriously biased estimates can increase even in RI-CLPM if model misspecifications exist.
Correct specification of the measurement model can be challenging. Fortunately, bothcurrent simulations and Usami (2022) have shown that even (misspecified) models with atime-varying AR(1) structure for WPVS produce accurate estimates, and that the choiceof order has little influence on estimation performance. However, further investigations ofestimation performance that account for various misspecifications in measurement models
are still required, and comparisons with other recent estimation approaches (e.g., Du &Bentler, 2021; Du, Bentler & Rosseel, 2022) for parameters in measurement models, especially for data with non-normal or small N , will also be an important topic for futurestudies.
Another important but still unresolved issue is how to establish the correct (or evena plausible) DAG model, or how one can validate the incorporation of time-invariant factors, such as stable trait factors, to infer within-person relations. To identify the causal
parameters, we assumed that measurements are expressed by the linear sum of stable trait
factor scores and WPVS, as in RI-CLPM. Usami, Murayama et al. (2019) explained that
there are two primary ways to statistically control for time-invariant factors as individual
diﬀerences: using stable trait factors included in RI-CLPM that have only direct eﬀects on
measurements and are uncorrelated with within-person processes (i.e., WPVS), and using
the accumulating factors included in several other statistical models (e.g., Bollen & Brand,2010), which have both direct and indirect eﬀects on measurements. In these models, accumulating factors are modeled with lagged regressions rather than being modeled separately

38

from them (see Usami 2022a for further details and the corresponding DGP).
The discussion of issues surrounding the appropriate model choice and the potential
diﬀerence of inferential results among statistical models when inferring within-person relations is ongoing, and it continues to gather attention from quantitative researchers in
psychology (e.g., Andersen, 2022; Lucas, 2022; Lüdtke & Robitzsch, 2021; Usami, 2022ab).
Though our attention in this paper was on TS-SNMMs that assume stable traits factors in
DGP, if an uncorrelated assumption of between- and within-person processes is violated,then both RI-CLPM and TS-SNMMs can produce seriously biased estimates. One canresort to substantial theory to choose the model (e.g., Shehata et al., 2021; who argued
that some statistical models that include stable trait factors, such as a constrained version of RI-CLPM, are well suited to model between- and within-person components when
capturing maintenance eﬀects in communication research). However, in many cases, researchers do not exactly know the true DGP and how time-invariant factors (if they exist)
influence measurements (e.g., linearly or nonlinearly, directly or indirectly, or both), andunambiguous specification of the theoretically derived expected relations for variables isquite challenging in practical applications (e.g., Curran & Bauer, 2011; Usami, 2022a).
Further discussion accounting for empirical investigations and the sensitivity of results will
be required in the future (Usami, 2022a). The extensions of TS-SNMMs that accountfor accumulating factors (i.e., correlated within-person and between-person processes) andmeasurement errors, along with developing a package to further increase the feasibility of
TS-SNMMs, are also important future research goals.

39

REFERENCES
Andersen, H.K. (2022). Equivalent approaches to dealing with unobserved heterogeneityin cross-lagged panel models? Investigating the benefits and drawbacks of the latent
curve model with structured residuals and the random intercept cross-lagged panel
model. Psychological Methods, https://doi.org/10.1037/met0000285
Ando, S., Nishida, A., Yamasaki, S., Koike, S., Morimoto, Y., Hoshino, A., Kanata, S.,
Fujikawa, S., Endo, K., Usami, S., Furukawa, T.A., Hiraiwa-Hasegawa, M., & Kasai,
K. (2019). Cohort profile: Tokyo Teen Cohort study (TTC). International Journalof Epidemiology, 48, 1414-1414g.
Angold, A., Costello, E.J., Messer, S.C., Pickles, A., Winder, F., & Silver, D. (1995).
Development of a short questionnaire for use in epidemiological studies of depressionin children and adolescents. International Journal of methods in Psychiatric Research,
5, 237-249.
Asparouhov, T., & Muthén, B. (2018). Latent variable centering of predictors andmediators in multilevel and time-series models. Structural Equation Modeling: A
Multidisciplinary Journal, 26, 1-24.
Bollen, K.A., & Brand, J.E. (2010). A general panel model with random and fixedeﬀects: A structural equations approach. Social Forces, 89, 1-34.
Curran, P.J., & Bauer, D.J. (2011). The disaggregation of within-person and betweenperson eﬀects in longitudinal models of change. Annual Review of Psychology, 62,
583-619.
Du, H., & Bentler, P.M. (2022). Distributionally weighted least squares in structural
equation modeling. Psychological Methods, 27, 519-540.
Du, H., Bentler, P.M. & Rosseel, Y. (2022). Distributionally-weighted least squares
in growth curve modeling. Structural Equation Modeling: A Multidisciplinary
Journal, 29, 1-22. DOI: 10.1080/10705511.2021.1931870
Hamaker, E.L., Kuiper, R.M., & Grasman, R.P.P.P. (2015). A critique of the cross-lagged40

panel model. Psychological Methods, 20, 102-116.
Hardin, J.W. & Hilbe, J.M. (2013) Generalized estimating equations (2nd ed). (pp.30-34)
Chapman and Hall/CRC Press, Boca Ratan, FL.
Hernán, M.A., & Robins, J.M. (2021). Causal Inference: What If. Boca Raton: Chapman& Hall/CRC.
Hoﬀman, L. (2014). Longitudinal analysis: Modeling within-person fluctuation and change.
New York, NY: Routledge/Taylor & Francis.
Hong, G. (2015). Causality in a social world: Moderation, mediation and spill-over.
West Sussex, UK: John Wiley & Sons, Ltd.
Kenny, D.A. & Judd, C.M. (1984). Estimating the nonlinear and interactive eﬀects oflatent variables. Psychological Bulletin, 96, 201-210.
Klein, A., & Moosbrugger, H. (2000). Maximum likelihood estimation of latent interactioneﬀects with the LMS method. Psychometrika, 65, 457-474.
Kline, R.B. (2016). Principles and practice of structural equation modeling (4th ed.).
Guilford Press.
Lee, S.Y. (2007). Structural Equation Modeling: A Bayesian Approach. John Wiley and
Sons.
Lucas, R.E. (2022). It’s time to abandon the cross-lagged panel model.
https://doi.org/10.31234/osf.io/pkec7
Lüdtke, O., & Robitzsch, A. (2021). A critique of the random intercept cross-lagged
panel model. https://doi.org/10.31234/osf.io/6f85c
Marsh H.W., Wen Z., & Hau, K.T. (2006). Structural equation models of latent interaction
and quadratic eﬀects. In Hancock G.R., Mueller R.O. (Eds.), Structural equation
modeling: A second course (pp. 225-265). Greenwich, CT: Information Age.
Matamura, M., Tochigi, M, Usami, S., Yonehara, H., Fukushima, M., Nishida, A.,
Togo, F., & Sasaki, T. (2014). Associations between sleep habits and mental healthstatus and suicidality in the longitudinal survey of monozygotic-twin adolescents.

41

Journal of Sleep Research, 23, 290-294.
Orth, U., Clark, D.A., Donnellan, M.B., & Robins, R.W. (2021). Testing prospective
eﬀects in longitudinal research: Comparing seven competing cross-lagged models.
Journal of Personality and Social Psychology, 120, 1013-1034.
Ozkok, O., Vaulont, M.J., Zyphur, M.J., Zhang, Z., Preacher, K.J., Koval, P., &Zheng, Y. (2022). Interaction eﬀects in cross-lagged panel models: SEM with latentinteractions applied to work-family conflict, job satisfaction, and gender.
Organizational Research Methods, 25, 673-715.
Robins, J.M. (1989). The analysis of randomized and non-randomized AIDS treatment
trials using a new approach to causal inference in longitudinal studies. In HealthService Research Methodology: A Focus on AIDS (L. Sechrest, H. Freeman and A.
Mulley, eds.) 113-159. U.S. Public Health Service, National Center for Health Services
Research, Washington, DC.
Robins, J.M. (1992). Estimation of the time-dependent accelerated failure time model inthe presence of confounding factors. Biometrika, 79, 321-334.
Robins, J.M. (1994). Correcting for non-compliance in randomized trials using structural
nested mean models. Communications in Statistics–Theory and Methods, 23,2379-2412.
Robins, J.M. (1999). Marginal structural models versus structural nested models as tools
for causal inference. Epidemiology, 116, 95-134.
Robins, J.M., & Hernán, M.A. (2009). Estimation of the causal eﬀects of time-varying
exposures. In G. Fitzmaurice et al. (Eds.), Handbooks of modern statistical methods:
Longitudinal data analysis (pp. 553-599). Boca Raton: CRC Press.
Robins, J.M., Hernán, M.A. & Brumback, B. (2000). Marginal structural models andcausal inference in epidemiology. Epidemiology, 11, 550-560.
Robins, J.M. & Rotnitzky, A. (2001). Comment on “Inference for semiparametric models:Some questions and an answer”, by P.J. Bickel and J. Kwon, Statistica Sinica, 11,42920-936.
Rosseel, Y. (2012). Lavaan: an R package for structural equation modeling. Journal of
Statistical software, 48, 1-36.
Shehata, A., Andersson, D., Glogger, I., Hopmann, D.N., Andersen, K., Kruikemeier, S.,& Johansson, J. (2021). Conceptualizing long-term media eﬀects on societal beliefs.
Annals of the International Communication Association, 45, 75-93.
ten Berge, T.M.F., Krijinen, W.P., Wansbeek, T.J., & Shapiro, A. (1999). Some new
results on correlation-preserving factor scores prediction methods. Linear Algebra andits Applications, 289, 311-318.
Usami, S. (2021). On the diﬀerences between general cross-lagged panel model and randomintercept cross-lagged panel model: Interpretation of cross-lagged parameters and model
choice. Structural Equation Modeling: A Multidisciplinary Journal, 28, 331-344.
Usami, S. (2022a). Within-person variability score-based causal inference: A two-stepestimation for joint eﬀects of time-varying treatments. Psychometrika. in press.
Usami, S. (2022b). Statistical models for the inference of within-person relations:
A random intercept cross-lagged panel model and its interpretation. The Japanese
Journal of Developmental Psychology. in press.
Usami, S., Murayama, K., & Hamaker, E.L. (2019). A unified framework of longitudinalmodels to examine reciprocal relations. Psychological Methods, 24, 637-657.
Usami, S., Todo, N., & Murayama, K. (2019). Modeling reciprocal eﬀects in medical
research: Critical discussion on the current practices and potential alternative models.
PLOS ONE, 14(9): e0209133.
VanderWeele, T.J., Hawkley, L.C., Thisted, R.A., Cacioppo, J.T. (2011). A marginal
structural model analysis for loneliness: implications for intervention trials and clinicalpractice. Journal of Consulting and Clinical Psychology, 79, 225-235.
Vansteelandt, S., & Joﬀe, M. (2014). Structural nested models and g-estimation: thepartially realized promise. Statistical Science, 29, 707-731.

43

Wallace, M.P., & Moodie, E.E.M. (2015). Doubly-robust dynamic treatment regimen
estimation via weighted least squares. Biometrics. 71, 636-644.
Wallace, M.P., Moodie, E.E.M., & Stephens, D.A. (2017a). An R package for G-estimationof structural nested mean models. Epidemiology, 28, e18-e20.
Wallace, M.P., Moodie, E.E.M., & Stephens, D.A. (2017b) Dynamic treatment regimen
estimation via regression-based techniques: Introducing R Package DTRreg.
Journal of Statistical Software, 80, 1-20.
Wang, L. & Maxwell, S.E. (2015). On disaggregating between-person and within-person
eﬀects with longitudinal data using multilevel models. Psychological Methods, 20,
63-83.

44

Table 1. Proportion of improper solutions in each condition.
T=3
2

𝜑𝜑 =0.12

𝜑𝜑 =0.4
𝜑𝜑 2 =0.7

T=4

T=6

T=9

RI-CLPMTS-SNMMs

RI-CLPMTS-SNMMs

RI-CLPMTS-SNMMs

RI-CLPMTS-SNMMs

N=200

86.53%

19.24%

76.82%

5.63%

58.00%

0.62%

37.15%

0.03%

N=600

74.45%

2.92%

55.45%

0.18%

28.27%

0.00%

10.11%

0.00%

N=100065.32%

0.66%

42.34%

0.00%

16.87%

0.00%

3.78%

0.00%

N=200

39.52%

1.39%

15.31%

0.01%

2.50%

0.00%

0.24%

0.00%

N=600

11.92%

0.00%

1.61%

0.00%

0.06%

0.00%

0.00%

0.00%

N=1000

5.67%

0.00%

0.45%

0.00%

0.00%

0.00%

0.00%

0.00%

N=200

26.01%

5.44%

7.46%

0.69%

0.63%

0.00%

0.04%

0.00%

N=600

2.23%

0.09%

0.20%

0.00%

0.00%

0.00%

0.00%

0.00%

N=1000
0.51%
0.00%
0.02%
0.00%
0.00%
0.00%
0.00%
0.00%
2
𝜑𝜑 , (proportion of) stable trait factor variances at t=1; T, the number of time points; N, sample size; RI-CLPM, random
intercept cross-lagged panel model; TS-SNMMs, two step estimation approach with structural nested mean model.

45

Table 2. Joint effect estimates for reciprocal relations in TS-SNMMs and the RI-CLPM (N=1,294)
TS-SNMMs
Estimate
SE-0.305
0.148
0.057
0.092
-0.0300.070-0.2890.249-0.601
0.629-0.173
0.080
0.300
0.514
0.178
1.324
0.104
0.095
Estimate
SE
0.000
0.007
0.017
0.017
0.010
0.018
0.018
0.007
-0.007
0.018
0.007
0.018
0.009
0.007
0.055
0.021
0.048
0.020

Sleep→SMFQ
Sleep 3 → SMFQ 4(Sleep 3 × Bedtime 3) → SMFQ 4(Sleep 3 × SMFQ 3) → SMFQ 4
Sleep 2 → SMFQ 4(Sleep 2 × Bedtime 2) → SMFQ 4(Sleep 2 × SMFQ 2) → SMFQ 4
Sleep 1 → SMFQ 4(Sleep 1 × Bedtime 1) → SMFQ 4(Sleep 1 × SMFQ 1) → SMFQ 4
SMFQ→Sleep
SMFQ 3 → Sleep 4(SMFQ 3 × Bedtime 3) → Sleep 4(SMFQ 3 × Sleep 3) → Sleep 4
SMFQ 2 → Sleep 4(SMFQ 2 × Bedtime 2) → Sleep 4(SMFQ 2 × Sleep 2) → Sleep 4
SMFQ 1 → Sleep 4(SMFQ 1 × Bedtime 1) → Sleep 4(SMFQ 1 × Sleep 1) → Sleep 4

RI-CLPM
Estimate
SE-0.9010.321-0.2550.157-0.071

0.059

Estimate
0.001

SE
0.008

0.000

0.003

0.002

0.001

Bold font indicates statistical significance. Sleep, sleep duration; SMFQ, Short Mood and Feelings Questionnaire; SE,standard error.

46

Figure 1. The linear causal diagrams (DAGs) for data-generating process in which stable trait
factors are included. Solid single-headed arrows (directed edges) are labeled with pathcoefficients that quantify direct causal effects. A dashed double-headed arrow (bidirected edge)
represents a correlation due to an unobserved common cause. Stable trait factors are represented
in dashed circles, indicating that these are latent variables.

47

𝜑𝜑 2 , (proportion of) stable trait factor variances at t=1; T, the number of time points; 𝑟𝑟𝐵𝐵 , correlation between stable

trait factors; 𝜔𝜔, ratio of variance of within-person variability scores at t=T to that at t=1; RI-CLPM, random intercept
cross-lagged panel model; TS-SNMMs, two-step estimation approach with structural nested mean model.

Figure 2. Averages of estimates and biases for the joint effects of X* on Y* in TS-SNMMs andthe RI-CLPM.

48

𝜑𝜑 2 , (proportion of) stable trait factor variances at t=1; T, the number of time points; 𝑟𝑟𝐵𝐵 , correlation between stable

trait factors; 𝜔𝜔, ratio of variance of within-person variability scores at t=T to that at t=1; RI-CLPM, random intercept
cross-lagged panel model; TS-SNMMs, two-step estimation approach with structural nested mean model.

Figure 3. Averages of estimates and biases for the joint effects of X* on Y* in TS-SNMMs andthe RI-CLPM when there exist unobserved confounders U* that influence the observed timevarying confounders L* and the outcomes Y*.

49

𝜑𝜑 2 , (proportion of) stable trait factor variances at t=1; T, the number of time points; 𝑟𝑟𝐵𝐵 , correlation between stable

trait factors; 𝜔𝜔, ratio of variance of within-person variability scores at t=T to that at t=1; RI-CLPM, random intercept
cross-lagged panel model; TS-SNMMs, two-step estimation approach with structural nested mean model.

Figure 4. Averages of estimates and biases for the joint effects of X* on Y* in TS-SNMMs andthe RI-CLPM when there exist ignored direct higher-order effects of observed time-varying
confounders L* on outcomes Y*.

50

Online Supplemental Material
Tables ············································································ p. 2
Figures ········································································ p. 11*R code for the simulation is provided in a separate file.

1

Table S1. Results of ANOVA where calculated biases are set as outcomes.
First-order effect of X* on Y*

Second-order effect of X* on Y*

T

3

0.01967

0.00656

𝜂𝜂 2

0.046

0.20790.06931

𝜂𝜂 2

N

2

0.000960.000480.002

0.01970.00984

0.009

2

2

0.01969

0.00984

0.046

0.12880.06440

0.061

𝑟𝑟𝐵𝐵

2

0.107410.053700.248

0.42730.21363

0.202

𝑟𝑟𝑊𝑊

2

0.070470.035230.163

0.09040.04519

0.043

ω

2

0.037810.018900.087

0.16940.08468

0.080

T×N

6

0.001130.000190.003

0.01930.003220.009

T×𝜑𝜑 2

df

𝜑𝜑

SS

MS

SS

MS

0.098

6

0.000390.00007

0.001

0.09690.01614

0.046

2

4

0.000640.000160.001

0.00140.00036

0.001

T×𝑟𝑟𝐵𝐵

6

0.019160.003190.044

0.23510.03919

0.111

N×𝑟𝑟𝐵𝐵

40.000540.000140.001

0.00030.00008

0.000

𝜑𝜑 2 ×𝑟𝑟𝐵𝐵

40.02665

0.006660.062

0.12970.03242

0.061

T×𝑟𝑟𝑊𝑊

6

0.01808

0.00301

0.042

0.01760.00294

0.008

N×𝑟𝑟𝑊𝑊

40.000690.000170.002

0.00010.00002

0.000

2

𝜑𝜑 ×𝑟𝑟𝑊𝑊

40.005640.001410.013

0.06470.01619

0.031

𝑟𝑟𝐵𝐵 ×𝑟𝑟𝑊𝑊

40.000250.00006

0.001

0.00360.00090

0.002

T×ω

6

0.00008

0.00001

0.000

0.11140.01856

0.053

N×ω

4

0.00004

0.00001

0.000

0.00300.00076

0.001

𝜑𝜑 ×ω

4

0.01232

0.00308

0.029

0.00340.00084

0.002

𝑟𝑟𝐵𝐵 ×ω

4

0.02303

0.005760.053

0.00030.00006

0.000

𝑟𝑟𝑊𝑊 ×ω

4

0.008190.00205

0.019

0.00470.00119

0.002

T×N×𝜑𝜑 2

120.002020.000170.005

0.01030.00086

0.005

T×N×𝑟𝑟𝐵𝐵

12

0.000810.00007

0.002

0.00060.00005

0.000

T×𝜑𝜑 2 ×𝑟𝑟𝐵𝐵

120.007300.00061

0.017

0.09660.00805

0.046

2

8

0.000240.00003

0.001

0.00060.00008

0.000

T×N×𝑟𝑟𝑊𝑊

N×𝜑𝜑 2 ×𝑟𝑟𝑊𝑊

N×𝜑𝜑

2

N×𝜑𝜑 ×𝑟𝑟𝐵𝐵

120.00091

0.00008

0.002

0.00100.00009

0.000

2

T×𝜑𝜑 ×𝑟𝑟𝑊𝑊

120.00360

0.00030

0.008

0.06520.005430.031

80.000270.00003

0.001

0.00060.00007

0.000

T×𝑟𝑟𝐵𝐵 ×𝑟𝑟𝑊𝑊

120.000280.00002

0.001

0.00600.00050

0.003

N×𝑟𝑟𝐵𝐵 ×𝑟𝑟𝑊𝑊

80.000480.00006

0.001

0.00070.00009

0.000

𝜑𝜑 ×𝑟𝑟𝐵𝐵 ×𝑟𝑟𝑊𝑊

80.000470.00006

0.001

0.00140.00017

0.001

T×N×ω

12

0.000800.00007

0.002

0.01030.000860.005

T×𝜑𝜑 ×ω

12

0.000680.00006

0.002

0.01090.00091

0.005

2

2

2

N×𝜑𝜑 2 ×ω

80.000440.00006

0.001

0.00010.00012

0.000

T×𝑟𝑟𝐵𝐵 ×ω

12

0.00400

0.00033

0.009

0.03910.00326

0.019

N×𝑟𝑟𝐵𝐵 ×ω

80.000210.000030.000

0.00080.000100.000

𝜑𝜑 2 ×𝑟𝑟𝐵𝐵 ×ω

8

0.005610.000700.013

0.00120.00015

0.001

T×𝑟𝑟𝑊𝑊 ×ω

12

0.002490.000210.006

0.03120.00260

0.015

N×𝑟𝑟𝑊𝑊 ×ω

8

0.00033

0.00004

0.001

0.00060.00008

0.000

𝜑𝜑 2 ×𝑟𝑟𝑊𝑊 ×ω

80.001130.000140.003

0.00300.00038

0.001

𝑟𝑟𝐵𝐵 ×𝑟𝑟𝑊𝑊 ×ω

80.000700.00009

0.002

0.00080.00009

0.000

T×N×𝜑𝜑 2 ×𝑟𝑟𝐵𝐵

240.00087

0.00004

0.002

0.00130.00005

0.001

T×N×𝜑𝜑 2 ×𝑟𝑟𝑊𝑊

240.00060

0.000020.001

0.00050.00002

0.000

T×N×𝑟𝑟𝐵𝐵 ×𝑟𝑟𝑊𝑊

240.000590.00002

0.001

0.00310.000130.001

T×𝜑𝜑 2 ×𝑟𝑟𝐵𝐵 ×𝑟𝑟𝑊𝑊

240.00101

0.00004

0.002

0.00250.000100.001

N×𝜑𝜑 2 ×𝑟𝑟𝐵𝐵 ×𝑟𝑟𝑊𝑊

160.000690.00004

0.002

0.00100.00006

0.000

T×N×𝜑𝜑 2 ×ω

24

0.003210.000130.007

0.00690.00029

0.003

T×N×𝑟𝑟𝐵𝐵 ×ω

24

0.000970.00004

0.002

0.00190.00008

0.001

T×𝜑𝜑 2 ×𝑟𝑟𝐵𝐵 ×ω

24

0.00157

0.00007

0.004

0.01580.00066

0.007

N×𝜑𝜑 2 ×𝑟𝑟𝐵𝐵 ×ω

16

0.000530.000030.001

0.00180.000110.001

T×N×𝑟𝑟𝑊𝑊 ×ω

24

0.00088

0.00004

0.002

0.00200.00008

0.001

T×𝜑𝜑 2 ×𝑟𝑟𝑊𝑊 ×ω

240.001250.00005

0.003

0.01420.00059

0.007

N×𝜑𝜑 2 ×𝑟𝑟𝑊𝑊 ×ω

160.000530.00003

0.001

0.00080.00005

0.000

T×𝑟𝑟𝐵𝐵 ×𝑟𝑟𝑊𝑊 ×ω

24

0.000630.00003

0.001

0.00420.00018

0.002

N×𝑟𝑟𝐵𝐵 ×𝑟𝑟𝑊𝑊 ×ω

160.00072

0.00005

0.002

0.00280.000170.001

𝜑𝜑 2 ×𝑟𝑟𝐵𝐵 ×𝑟𝑟𝑊𝑊 ×ω

160.000630.00004

0.001

0.00050.00003

0.000

T×N×𝜑𝜑 2 ×𝑟𝑟𝐵𝐵 ×𝑟𝑟𝑊𝑊

48

0.00158

0.00003

0.004

0.00240.00005

0.001

T×N×𝜑𝜑 2 ×𝑟𝑟𝐵𝐵 ×ω

48

0.002100.00004

0.005

0.00660.000140.003

T×N×𝜑𝜑 2 ×𝑟𝑟𝑊𝑊 ×ω

48

0.001140.00002

0.003

0.00430.00009

0.002

T×N×𝑟𝑟𝐵𝐵 ×𝑟𝑟𝑊𝑊 ×ω

48

0.001960.00004

0.005

0.00650.000140.003

T×𝜑𝜑 2 ×𝑟𝑟𝐵𝐵 ×𝑟𝑟𝑊𝑊 ×ω

48

0.001350.00003

0.003

0.00310.00006

0.001

N×𝜑𝜑 2 ×𝑟𝑟𝐵𝐵 ×𝑟𝑟𝑊𝑊 ×ω

32

0.001110.00003

0.003

0.00270.00009

0.001

T×N×𝜑𝜑 2 ×𝑟𝑟𝐵𝐵 ×𝑟𝑟𝑊𝑊 ×ω

96

0.002730.00003

0.006

0.01030.00011

0.005

Total

971

0.43226

2.1121

3

Table S2. Proportions of improper solutions in each condition (second simulation).
T=32

𝜑𝜑 =0.12

𝜑𝜑 =0.4
𝜑𝜑 2 =0.7

T=4

T=6

T=9

RI-CLPMTS-SNMMs

RI-CLPMTS-SNMMs

RI-CLPMTS-SNMMs

RI-CLPMTS-SNMMs

N=200

85.59%

14.75%

76.63%

3.60%

65.26%

0.33%

55.52%

0.01%

N=600

72.79%

1.90%

57.13%

0.07%

44.93%

0.00%

39.11%

0.00%

N=100063.13%

0.40%

46.52%

0.00%

37.11%

0.00%

32.64%

0.00%

N=200

43.33%

2.34%

23.14%

0.13%

9.49%

0.00%

2.59%

0.00%

N=600

13.26%

0.05%

4.52%

0.00%

0.96%

0.00%

0.02%

0.00%

N=1000

6.44%

0.00%

1.68%

0.00%

0.18%

0.00%

0.00%

0.00%

N=200

31.89%

9.46%

13.93%

2.30%

3.50%

0.00%

0.42%

0.00%

N=600

5.67%

0.61%

0.56%

0.04%

0.00%

0.00%

0.00%

0.00%

N=1000

1.32%

0.07%

0.22%

0.02%

0.00%

0.00%

0.00%

0.00%

RI-CLPM, random intercept cross-lagged panel model; TS-SNMMs, two-step estimation approach with astructural nested mean model.

4

Table S3. Proportions of observed values that exceeded prespecified criteria for model fit indices (second simulation).
T=3

𝜑𝜑 2=0.12

𝜑𝜑 =0.42

𝜑𝜑 =0.7

T=4

T=6

T=9

SRMR

CFI

RMSEA

SRMR

CFI

RMSEA

SRMR

CFI

RMSEA

SRMR

CFI

RMSEA

<.05

>.95

<.05

<.05

>.95

<.05

<.05

>.95

<.05

<.05

>.95

<.05

N=200

100.00%99.80%

80.24%

99.89%

99.74%

93.04%

98.07%

99.69%

99.41%

79.24%

99.72%

100.00%N=600

100.00%100.00%94.11%

100.00%100.00%99.98%

100.00%100.00%100.00%100.00%100.00%100.00%N=1000100.00%100.00%98.11%

100.00%100.00%100.00%100.00%100.00%100.00%100.00%100.00%100.00%N=200

100.00%99.98%

80.54%

99.81%

100.00%93.67%

97.11%

100.00%99.31%

78.65%

100.00%100.00%N=600

100.00%100.00%93.78%

100.00%100.00%99.96%

100.00%100.00%100.00%100.00%100.00%100.00%N=1000100.00%100.00%98.37%

100.00%100.00%100.00%100.00%100.00%100.00%100.00%100.00%100.00%N=200

100.00%100.00%83.19%

99.85%

100.00%94.20%

95.59%

100.00%99.20%

81.81%

100.00%100.00%N=600

100.00%100.00%94.54%

100.00%100.00%99.96%

100.00%100.00%100.00%100.00%100.00%100.00%N=1000100.00%100.00%98.56%

100.00%100.00%100.00%100.00%100.00%100.00%100.00%100.00%100.00%

CFI, comparative fit index; RMSEA, root mean square error of approximation; SRMR, standardized root
mean square residual.

5

Table S4. Proportions of improper solutions in each condition (third simulation).
T=32

𝜑𝜑 =0.12

𝜑𝜑 =0.4
𝜑𝜑 2 =0.7

T=4

T=6

T=9

RI-CLPMTS-SNMMs

RI-CLPMTS-SNMMs

RI-CLPMTS-SNMMs

RI-CLPMTS-SNMMs

N=200

87.88%

17.05%

79.12%

3.57%

59.02%

0.18%

33.57%

0.01%

N=600

82.67%

2.62%

66.47%

0.07%

32.26%

0.00%

9.83%

0.00%

N=100080.42%

0.61%

60.35%

0.00%

25.02%

0.00%

5.42%

0.00%

N=200

36.28%

2.15%

10.74%

0.04%

1.54%

0.00%

0.09%

0.00%

N=600

5.37%

0.00%

0.24%

0.00%

0.00%

0.00%

0.00%

0.00%

N=1000

1.10%

0.00%

0.00%

0.00%

0.00%

0.00%

0.00%

0.00%

N=200

26.26%

6.62%

11.61%

1.32%

2.94%

0.00%

0.52%

0.00%

N=600

3.62%

0.16%

0.69%

0.00%

0.00%

0.00%

0.00%

0.00%

N=1000

1.16%

0.04%

0.16%

0.00%

0.00%

0.00%

0.00%

0.00%

RI-CLPM, random intercept cross-lagged panel model; TS-SNMMs, two-step estimation approach with astructural nested mean model.

6

Table S5. Proportions of observed values that exceeded prespecified criteria for model fit indices (third simulation).
T=3

𝜑𝜑 2=0.12

𝜑𝜑 =0.42

𝜑𝜑 =0.7

T=4

T=6

T=9

SRMR

CFI

RMSEA

SRMR

CFI

RMSEA

SRMR

CFI

RMSEA

SRMR

CFI

RMSEA

<.05

>.95

<.05

<.05

>.95

<.05

<.05

>.95

<.05

<.05

>.95

<.05

N=200

99.98%

99.31%

60.00%

96.98%

96.94%

69.81%

87.06%

92.59%

86.20%

33.67%

81.33%

98.54%

N=600

100.00%100.00%52.24%

100.00%99.94%

77.94%

100.00%99.80%

98.74%

100.00%99.76%

100.00%N=1000100.00%100.00%52.46%

100.00%99.98%

80.15%

100.00%100.00%99.63%

100.00%100.00%100.00%N=200

100.00%100.00%67.35%

98.09%

99.78%

76.63%

87.89%

99.35%

88.48%

47.80%

98.56%

98.46%

N=600

100.00%100.00%66.37%

100.00%100.00%89.57%

100.00%100.00%99.65%

100.00%100.00%100.00%N=1000100.00%100.00%64.28%

100.00%100.00%93.26%

100.00%100.00%99.98%

100.00%100.00%100.00%N=200

100.00%100.00%75.41%

99.17%

100.00%85.39%

88.26%

100.00%92.83%

68.19%

100.00%98.93%

N=600

100.00%100.00%81.04%

100.00%100.00%97.20%

99.80%

100.00%99.89%

99.93%

100.00%100.00%N=1000100.00%100.00%84.59%

100.00%100.00%99.00%

100.00%100.00%100.00%100.00%100.00%100.00%

CFI, comparative fit index; RMSEA, root mean square error of approximation; SRMR, standardized root
mean square residual.

7

Table S6. Descriptive statistics.
(a) Mean and standard deviation
BT10

SD10

SMFQ10

BT12

SD12

SMFQ12

BT14

SD14

SMFQ14

BT16

SD16

SMFQ16

Mean

21.71

9.03

4.49

22.25

8.48

3.67

23.10

7.65

3.03

23.72

7.01

3.66

Standard
deviation

0.63

0.63

4.46

0.73

0.74

4.30

0.81

0.86

4.53

0.86

0.90

5.19

(b) Variances, covariances, and correlations
BT10

SD10

SMFQ10

BT12

SD12

SMFQ12

BT14

SD14

SMFQ14

BT16

SD16

SMFQ16

GB10

0.40

-0.75

0.11

0.64

-0.46

0.06

0.38

-0.28

0.09

0.24

-0.16

0.10

SD10

-0.30

0.39

-0.04

-0.44

0.55

-0.04

-0.22

0.34

-0.05

-0.14

0.17

-0.06

SMFQ10

0.30

-0.11

19.90

0.04

0.00

0.41

0.04

0.04

0.27

0.02

0.05

0.18

GB12

0.29

-0.20

0.12

0.53

-0.77

0.07

0.48

-0.40

0.10

0.26

-0.19

0.12

SD12

-0.21

0.25

0.02

-0.42

0.55

-0.06

-0.32

0.48

-0.08

-0.16

0.22

-0.11

SMFQ12

0.16

-0.10

7.87

0.23

-0.19

18.51

0.05

0.00

0.36

0.02

0.06

0.29

GB14

0.19

-0.11

0.16

0.28

-0.19

0.16

0.65

-0.67

0.14

0.47

-0.32

0.12

SD14

-0.15

0.18

0.13

-0.25

0.30

0.00

-0.47

0.74

-0.08

-0.28

0.43

-0.12

SMFQ14

0.26

-0.15

5.39

0.34

-0.27

7.12

0.53

-0.29

20.58

0.11

-0.02

0.44

GB16

0.13

-0.08

0.09

0.16

-0.10

0.06

0.33

-0.20

0.42

0.75

-0.60

0.13

SD16

-0.09

0.09

0.19

-0.12

0.15

0.21

-0.23

0.33

-0.06

-0.47

0.81

-0.07

0.32
-0.19
4.05
0.47
-0.42
6.38
0.49
-0.52
10.33
0.59
-0.31
* BT, bedtime; SD, sleep duration; SMFQ, Short Mood and Feelings Questionnaire.
Variances are in the diagonal elements, covariances are in the lower-left elements, and correlations arein the upper-right elements.

26.92

SMFQ16

8

Table S7. Fit indices of the AR(1) measurement model for each variable.
BedtimeSleep duration

SMFQ

No. of parameters

8

8

8

df

2

2

2

chi-square

2.507

7.325

7.788

p-value

0.285

0.026

0.02

1

0.995

0.992

RMSEA

0.014 [0.000, 0.059]

0.045 [0.014, 0.083]

0.047 [0.016, 0.084]

SRMR

0.011

0.018

0.018

CFI

Model fits are perfect for the AR(2) measurement model.

9Table S8. Estimates from TS-SNMMs with AR(2) measurement models (N=1,294).
TS-SNMMs
Sleep→SMFQ

Estimates

SE

Sleep 3 → SMFQ 4-0.2900.143

Sleep 3 × Bedtime 3 → SMFQ 4

0.033

0.094

Sleep 3 × SMFQ 3 → SMFQ 4-0.0270.066

Sleep 2 → SMFQ 4-0.1980.258

Sleep 2 × Bedtime 2 → SMFQ 4-0.6550.620

Sleep 2 × SMFQ 2 → SMFQ 4-0.161

0.074

Sleep 1 → SMFQ 4

0.382

0.501

Sleep 1 × Bedtime 1 → SMFQ 4

0.223

1.260

Sleep 1 × SMFQ 1 → SMFQ 4

0.101

0.091

Estimates

SE

SMFQ 3 → Sleep 4-0.0030.007

SMFQ 3 × Bedtime 3 → Sleep 4

0.016

0.015

SMFQ 3 × Sleep 3 → Sleep 4

0.008

0.015

SMFQ 2 → Sleep 4

0.015

0.007

SMFQ 2 × Bedtime 2 → Sleep 4-0.0090.017

SMFQ 2 × Sleep 2 → Sleep 4

0.003

0.016

SMFQ 1 → Sleep 4

0.004

0.007

SMFQ 1 × Bedtime 1 → Sleep 4

0.050

0.022

SMFQ 1 × Sleep 1 → Sleep 4

0.043

0.020

SMFQ→Sleep

Bold font indicates statistical significance.

10

Figure S1. RMSEs of estimates for the joint effects of X* on Y* in TS-SNMMs and RI-CLPM.
(N=1000 and 𝑟𝑟𝑊𝑊 = 0.3)

11

Figure S2. Averages of estimates and biases for the joint effects of X* on Y* in TS-SNMMs and RI-CLPM.

(N=200 and 𝑟𝑟𝑊𝑊 = 0.3)

12Figure S3. Averages of estimates and biases for the joint effects of X* on Y* in TS-SNMMs and RI-CLPM.

(N=600 and 𝑟𝑟𝑊𝑊 = 0.3)

13

Figure S4. Averages of estimates and biases for the joint effects of X* on Y* in TS-SNMMs and RI-CLPM.

(N=1000 and 𝑟𝑟𝑊𝑊 = 0.1)

14

Figure S5. Averages of estimates and biases for the joint effects of X* on Y* in TS-SNMMs and RI-CLPM.

(N=1000 and 𝑟𝑟𝑊𝑊 = 0.5)

15Figure S6. RMSEs of estimates for the joint effects of X* on Y* in TS-SNMMs and RI-CLPM.

(N=200 and 𝑟𝑟𝑊𝑊 = 0.3)

16Figure S7. RMSEs of estimates for the joint effects of X* on Y* in TS-SNMMs and RI-CLPM.

(N=600 and 𝑟𝑟𝑊𝑊 = 0.3)

17Figure S8. RMSEs of estimates for the joint effects of X* on Y* in TS-SNMMs and RI-CLPM.

(N=1000 and 𝑟𝑟𝑊𝑊 = 0.1)

18

Figure S9. RMSEs of estimates for the joint effects of X* on Y* in TS-SNMMs and RI-CLPM.

(N=1000 and 𝑟𝑟𝑊𝑊 = 0.5)

19Figure S10. Averages of estimates and biases for the joint effects of Y* on X* in TS-SNMMs and RI-CLPM.

(N=1000 and 𝑟𝑟𝑊𝑊 = 0.3)

20Figure S11. RMSEs of estimates for the joint effects of Y* on X* in TS-SNMMs and RI-CLPM.

(N=1000 and 𝑟𝑟𝑊𝑊 = 0.3)

21

Figure S12. Averages of estimates and biases for the joint effects of X* on Y* in TS-SNMMs and RI-CLPM.

(AR(2) measurement model; N=1000 and 𝑟𝑟𝑊𝑊 = 0.3)

22Figure S13. RMSEs of estimates for the joint effects of X* on Y* in TS-SNMMs and RI-CLPM.

(AR(2) measurement model; N=1000 and 𝑟𝑟𝑊𝑊 = 0.3)

23

Figure S14. Averages of estimates and biases for the joint effects of X* on Y* in TS-SNMMs and RICLPM when there exist unobserved confounders U* that influence observed time-varying
confounders L* and outcomes Y* (AR(2) measurement model; N=1000 and 𝑟𝑟𝑊𝑊 = 0.3).

24

Figure S15. Averages of estimates and biases for the joint effects of Y* on X* in TS-SNMMs and RICLPM when there exist unobserved confounders U* that influence observed time-varying
confounders L* and outcomes Y* (N=1000 and 𝑟𝑟𝑊𝑊 = 0.3).

25

Figure S16. Averages of estimates and biases for the joint effects of X* on Y* in TS-SNMMs and RICLPM when there exist ignored direct higher-order effects of observed time-varying confounders
L* on outcomes Y* (AR(2) measurement model; N=1000 and 𝑟𝑟𝑊𝑊 = 0.3).

26

Figure S17. Averages of estimates and biases for the joint effects of Y* on X* in TS-SNMMs and RICLPM when there exist ignored direct higher-order effects of observed time-varying confounders
L* on outcomes Y* (N=1000 and 𝑟𝑟𝑊𝑊 = 0.3).

27