# Marginal Structural Models with Latent Class Growth Modeling of Treatment
Trajectories

Awa Diop
<EMAIL_ADDRESS>
Département de médecine sociale et préventive
Université Laval &Caroline Sirois
<EMAIL_ADDRESS>
Faculté de pharmacie
Université Laval &Jason Robert Guertin
<EMAIL_ADDRESS>
Département de médecine sociale et préventive
Université Laval &Denis Talbot
<EMAIL_ADDRESS>
Département de médecine sociale et préventive
Université Laval

###### Abstract

In a real-life setting, little is known regarding the effectiveness of statins
for primary prevention among older adults, and analysis of observational data
can add crucial information on the benefits of actual patterns of use. Latent
class growth models (LCGM) are increasingly proposed as a solution to
summarize the observed longitudinal treatment in a few distinct groups. When
combined with standard approaches like Cox proportional hazards models, LCGM
can fail to control time-dependent confounding bias because of time-varying
covariates that have a double role of confounders and mediators. We propose to
use LCGM to classify individuals into a few latent classes based on their
medication adherence pattern, then choose a working marginal structural model
(MSM) that relates the outcome to these groups. The parameter of interest is
nonparametrically defined as the projection of the true MSM onto the chosen
working model. The combination of LCGM with MSM is a convenient way to
describe treatment adherence and can effectively control time-dependent
confounding. Simulation studies were used to illustrate our approach and
compare it with unadjusted, baseline covariates-adjusted, time-varying
covariates adjusted and inverse probability of trajectory groups weighting
adjusted models. We found that our proposed approach yielded estimators with
little or no bias.

Keywords : Latent class growth models, Marginal stuctural models, Projection-
based approach, Nonparametric, Inverse probability of treatment weighting,
Cardiovascular disease, Statins, time-varying treatment, time-dependent
confounding.

##  1 Introduction

Cardiovascular disease (CVD) is the leading cause of death worldwide ([24]). A
well known treatment to prevent CVDs is HMG-CoA reductase inhibitors
(statins). Several randomized controlled trials (RCTs) have shown the efficacy
of statins to prevent a first event of CVD, that is, for primary prevention
([20]). Despite current evidence, it is still not clear if these conclusions
can be applied to older adults as they are often excluded from RCTs ([18, 1]).
Also, treatment adherence is known to be lower in real-world situations than
in RCTs ([3]). Hence, there is little evidence concerning the effectiveness of
statins for primary prevention among older adults in a real-life setting.
Analysis of observational data could add crucial information on the benefits
of actual statin’s patterns of use among older adults. However, the number of
unique treatment trajectories increases exponentially with the length of
follow-up in longitudinal studies.

Latent class growth models (LCGM), or group-based trajectory models, are
increasingly being proposed as a solution to summarize the observed
trajectories in a few distinct groups ([13, 4]). LCGM defines homogeneous
subgroups of individuals with respect to their patterns of change over time
([13, 12]). Recently [4] applied LCGM to classify patients who initiated a
statin according to their long-term adherence. The authors compared LCGM with
traditional approaches to measure medication adherence such as the proportion
of days covered or medication possession ratio. They concluded that LCGM is a
better method to measure adherence and may facilitate targeting of
interventions as it made describing adherence behaviors easier. [5] have
combined LCGM and Cox models to describe the association between statin
adherence and cardiovascular events.

While LCGM is a promising approach for estimating the effect of treatment
trajectories in real-life settings, adjustment for potential confounders is
essential. Indeed, time-dependent confounders are susceptible to be affected
by previous treatments ([7]). As a consequence, time-varying covariates can
have a double-role of confounders and mediators. In such situations, standard
approaches, such as covariate adjustment in an outcome regression model, fail
to control confounding ([7]). Marginal structural models (MSMs) are well known
in the causal inference literature for their ability to deal correctly with
time-varying covariates and time-varying exposure.

Estimating the parameters of a MSM that links the outcome to all possible
treatment trajectories over long follow-up periods is unlikely to be realistic
in most applications. A common alternative is to model the outcome as a
function of the cumulative number of periods of treatment adherence. However
this approach lacks flexibility since it assumes a linear relationship between
cumulative treatment and outcome. A flexible Cox MSM with such weighted
cumulative exposure was proposed by [25]. The weighted cumulative exposure is
estimated with splines allowing a data-driven specification of the relation
between the outcome and the treatment trajectory. While this alternative is
flexible, it is not as easy to describe treatment adherence with it. An
attempt to combine LCGM with MSM was realized by [11] but the authors did not
use an appropriate estimation method to control for time-varying confounders.
Their method is similar to combining LCGM with standard modeling approachs.

In this paper, we propose a suitable theoretical framework to combine LCGM and
MSMs. We first introduce a finite mixture modeling approach (LCGM) to classify
individuals into a few latent classes. Our approach consists in choosing a
working MSM that relates the outcome to these latent classes. The parameter of
interest is nonparametrically defined as the projection of the true MSM onto
the chosen working causal model ([15]). Hence, we alleviate the burden of
assuming a correctly specified model. As far as we know, such a combination of
MSMs with LCGM that allows to control for time-varying confounders, has never
been investigated. In the remainder, we present the modeling framework
(notation, causal model, mixture model, estimation approach) followed by a
simulation study to illustrate and investigate the finite-sample properties of
our proposal. We end this paper with some recommendations for the practical
application of our proposed approach.

##  2 Notation and nonparametric marginal structural model definition

In the following, we use capital letters to represent random variables and
corresponding lower case letters to represent observed or fixed values these
variables take. Consider a longitudinal study where $t=1,\ldots,K$ is the
follow-up time. Let $\bar{A}_{t}=(A_{1},A_{2},\ldots,A_{t})$ be the treatment
trajectory up to time $t$. Similarly, $\bar{L}_{t}=(L_{1},L_{2},...L_{t})$ is
the covariates’ history up to time $t$. These covariates are assumed to
include all baseline and time-varying confounders. As a notational shortcut,
we define $\bar{A}\equiv\bar{A}_{K}$ and $\bar{L}\equiv\bar{L}_{K}$. Let $Y$
be the observed outcome at the end of follow-up. We denote by $Y^{\bar{a}}$
the counterfactual outcome under a specific treatment trajectory, that is, the
outcome that would have been observed had the treatment trajectory been
$\bar{a}$, possibly contrary to fact. Similarly, $\bar{L}^{\bar{a}}$ denotes a
counterfactual covariate process. We denote by $F_{X}$ the unknown
distribution of all possible counterfactual variables
$X=(\bar{L}^{\bar{a}},Y^{\bar{a}}:{\bar{a}\in\mathcal{A}})$ where
$\mathcal{A}$ is the set of all possible values of $\bar{a}$. We denote by
$P_{F_{X}}$ the distribution of the observed data
$\mathcal{O}=\\{Y,\bar{A}_{K},\bar{L}_{K}\\}$, from which $i=1,...,n$
independent and identically distributed observations are randomly drawn. We
consider the following nonparametric form of a MSM where we define $V$ as a
(possibly empty) subset of the baseline covariates, $V\subseteq L_{1}$:

$E_{F_{X}}(Y^{\bar{a}}|V)=m^{*}(\bar{a},V).$ (1)

Equation (1) represents a model where the distribution of all possible
counterfactual variables, $F_{X}$, is left unknown. Thus, we do not make the
assumption of a correctly specified model as would be the case when
considering a parametric MSM. The nonparametric identification of model (1)
from the observed data can be achieved under the following causal assumptions:
1) sequential conditional exchangeability, 2) no interaction between subjects,
3) consistency, and 4) positivity. Assumption 1) indicates that conditional on
the history of treatment up to time $t-1$ and the history of the covariates up
to time $t$ there are no unmeasured confounders ([17]), that is
$Y^{\bar{a}}\bm{\perp\\!\\!\\!\perp}A_{t}|\bar{A}_{t-1},\bar{L}_{t}$.
Assumption 2) means that the potential outcome of a given individual is not
affected by others’ exposure. Assumption 3) means that, given the observed
treatment trajectory, the observed outcome and the potential outcome under
that given trajectory are the same, formally if $\bar{A}=\bar{a}$ then
$Y^{\bar{a}}=Y$. Assumption 4) requires that for each level
$(\bar{a}_{t},\bar{l}_{t})$,
$P(A_{t}=a_{t}|\bar{A}_{t-1}=\bar{a}_{t-1},\bar{L}_{t}=\bar{l}_{t})>0$. In
other words, in each stratum defined by previous treatment and covariates, we
find both exposed and unexposed individuals at time $t$. Under these
assumptions, we have:

$\displaystyle\mathbb{E}_{F_{X}}\left(Y^{\bar{a}}\right)=\int_{\bm{l_{1}}}\ldots\int_{\bm{l_{K}}}\mathbb{E}(Y|\bar{a}_{K},\bar{\bm{l}}_{K})f(\bm{l_{K}|\bm{l}_{K-1}},\bar{a}_{K-1})\mu(\bm{L_{K}}),\ldots
f(\bm{l_{1}})\mu(\bm{L_{1}}),$ (2)

where $f(\bm{l_{K}|\bm{l}_{K-1}},\bar{a}_{K-1}),\ldots f(\bm{l_{1}})$ are the
densities of $\bm{\bar{L}_{t}}$ relative to appropriate measures
$\mu(\bm{L_{t}}),t=1,\ldots,K$, respectively. Due to the curse of
dimensionality, Equation (2) is rarely estimable with finite sample data. As
mentioned in the introduction, one challenge is related to the number of
treatment trajectories. In the next two sections, we introduce the LCGM and
our LCGM-MSM as a solution to this specific challenge.

##  3 Summarizing treatment trajectories with latent class growth models

We briefly review LCGM below, more details can be found elsewhere (for
example, see [13]). The first step of our LCGM-MSM proposal consists in
summarizing treatment trajectories with a LCGM. Let $C_{i}$ denote the latent
group-membership or trajectory group for the $i^{th}$ individual estimated
from the individual’s observed treatment $\bar{A}_{i},i=1,\ldots,n$. We denote
by $\pi_{j}$ the probability that an individual belongs to the $j^{th}$
trajectory group in the population, $j=1,\ldots,J$ where $J$ is the total
number of groups with the constraints $\pi_{j}\geq 0$ and
$\sum_{j=1}^{J}\pi_{j}=1$. In LCGM, a trajectory is often described by a
polynomial function of time. For example, if a quadratic relation was assumed
for the $j^{th}$ group, we could use the regression model
$logit(P(\bar{A}_{it}|C_{i}=j))=\theta_{0}^{j}+\theta_{1}^{j}t+\theta_{2}^{j}t^{2}$
where the parameters $\theta^{j},j=1,\ldots,J$ describe the treatment
trajectory over time for the $j^{th}$ group. This specification implies that
the following mixture model is the marginal probability of the treatment
trajectories $P(\bar{A}_{i})$:

$P(\bar{A}_{i})=\sum_{j=1}^{J}\pi_{j}\prod_{t=1}^{K}P(A_{it}|C_{i}=j).$ (3)

Model (3) relies on the assumption of local independence, that is treatment
uses at different time-points $A_{i1},\ldots,A_{iK}$ are assumed to be
mutually independent conditional on group-membership ([23]):
$A_{it}\perp\\!\\!\\!\perp A_{it^{{}^{\prime}}}|C_{i}=j,\ \text{for all}\
t\neq t^{{}^{\prime}},j=1,\ldots,J$. The log-likelihood of the $n$
observations is given by:

$\displaystyle\mathcal{L}(\pi,\theta)$
$\displaystyle=\sum_{i=1}^{n}log\left(\sum_{j=1}^{J}\pi_{j}\prod_{t=1}^{K}P(A_{it}|C_{i}=j,\theta^{j})\right).$

A general framework for fitting LCGM is implemented in the R package FlexMix
which is based on an Expectation-Maximization algorithm ([10]). We denote by
$Z_{ij}=P(C_{i}=j|\bar{A}_{i})$, the $i^{th}$ individual’s conditional
(“posterior”) probability of membership in group $j$. These probabilities are
computed postmodel estimation and are used to assign subjects into a
trajectory group, usually based on their largest value ([13]). Therefore, we
can formally define the trajectory groups as:
$C_{i}=\underset{j=1,\ldots,J}{argmax}\ Z_{ij}$. The conditional probability
is given by:

$P(C_{i}=j|\bar{A}_{i})=\dfrac{\pi_{j^{{}^{\prime}}}P(\bar{A}_{i}|C_{i}=j^{{}^{\prime}})}{\sum_{j=1}^{J}\pi_{j}P(\bar{A}_{i}|C_{i}=j)}.$
(4)

In Figure (1), we present a global overview of the LCGM approach. From the
matrix $n\times K$ of treatment trajectories, we determine the $n\times J$
matrix of individual conditional probabilities which in turn yields the vector
$n\times 1$ of trajectory groups. The specification of the order of the
polynomial equation and the number of trajectory groups to include in the
model are discussed in [13]. Note that the local independence assumption of
the LCGM is unlikely to hold in presence of time-dependent confounding.
Indeed, time-dependent confounding is susceptible to create a complex
dependence structure between treatment uses at different time-points. If this
assumption does not hold, LCGM can still be regarded as a convenient
statistical technique to summarize treatment trajectories.

$\ \underbracket{\blockarray{cccccc}A_{i1}&A_{i2}A_{i3}\ldots A_{iK}\\\
\block{(ccccc)c}a_{11}a_{12}a_{13}\ldots a_{1K}\\\ a_{21}a_{22}a_{23}\ldots
a_{2K}\\\ \vdots\ldots\vdots\\\ a_{n1}a_{n2}a_{n3}\ldots a_{nK}\\\
}_{\bm{n}\times\bm{K}}\begin{matrix}\\\
\xrightarrow{\sum_{j=1}^{J}\pi_{j}P(\bar{A}_{i}|C_{i}=j)}\\\
P(C_{i}=j|\bar{A}_{i})\end{matrix}\underbracket{\blockarray{ccccc}Z_{i1}&Z_{i2}\ldots
Z_{iJ}\\\ \block{(cccc)c}z_{11}z_{12}\ldots z_{1J}\\\ z_{21}z_{22}\ldots
z_{2J}\\\ \vdots\ldots\vdots\\\ z_{n1}z_{n2}\ldots z_{nJ}\\\
}_{\bm{n}\times\bm{J}}\begin{matrix}\\\
\xrightarrow{\underset{j=1,\ldots,J}{argmax}\
Z_{ij}}\end{matrix}\underbracket{\blockarray{cc}C&\\\ \block{(c)c}C_{1}\\\
C_{2}\\\ \vdots\\\ C_{n}\\\ }_{\bm{n}\times\bm{1}}$ Figure 1: Summarizing
treatment trajectories with LCGM.The matrix $n\times K$ represents treatment
values for all the $n$ individuals and $K$ follow-up time, the matrix $n\times
J$ represents matrix of $n$ individuals and their conditional probabilities,
the vector $n\times 1$ represents the trajectory groups. The formulas on the
arrows indicate how to obtain these matrices.

##  4 Defining the causal parameter of interest

Once the observed treatment trajectories $\bar{A}$ are summarized using the
LCGM, the next step of our LCGM-MSM approach is to formally define the causal
effect of interest. We propose to approximate the nonparametric MSM
$m^{*}(\bar{a},V)$ for all $\bar{a}\in\mathcal{A}$ with a working MSM that is
a function of the trajectory groups. Following a similar notation to [14], we
denote by $\mathcal{M}^{*}$ the infinite-dimensional space of all possible
functions $m^{*}$ relating the counterfactual expectation with the treatment
trajectories and baseline covariates $V$ as in Equation (1). We define
$\mathcal{M}$ a subset of $\mathcal{M^{*}}$, consisting of all the the family
of functions $m$ that the user is interested in using to characterize how the
counterfactual expectations are related to the treatment trajectories and
baseline covariates $V$ (i.e., the working model). Our MSM-LCGM approach
entails specifying $\mathcal{M}=\\{m(C,V|\beta):\beta\in\mathcal{R}^{p}\\}$,
where $p$ is the number of parameters in $m(C,V|\beta)$. In other words, we
are interested in characterizing the counterfactual outcomes as a function of
the LCGM trajectory groups, for example
$\mathbb{E}[Y|C,V]=\beta_{0}+\beta_{1}C_{1}+...+\beta_{J-1}C_{J-1}$. The model
could feature additional parameters for covariates $V$ and for interactions
terms between $C$ and $V$ if effect modification is of interest. Remark that
if $V=\emptyset$ and $J=2^{K}$, then $\mathcal{M}=\mathcal{M}^{*}$. Indeed,
the number of LCGM trajectory groups would then correspond with the total
number of possible treatment trajectories, and the working model would be
saturated.

Our causal parameter of interest is defined as a projection of the true
nonparametric MSM $m^{*}$ on the working model $m$. This projection is notably
characterized by a weight function $\lambda(\bar{A},V)$ that determines, for
all $(\bar{a},v)\in\mathcal{A}\times\mathcal{S}_{V}$ where $\mathcal{S}_{V}$
is the set of all possible values of $V$ (its support), how well the working
MSM $m$ should approximate the true MSM $m^{*}$ relative to other values of
$(\bar{a},v)$. More formally:

$\displaystyle\beta(\cdot|m,\lambda):\mathcal{M^{*}}\rightarrow\mathcal{R}^{p}$
(5)
$\displaystyle\beta_{m,\lambda}=\underset{\beta\in\mathbf{R}^{p}}{argmin}\lVert
m^{*}-m(C,V|\beta)\rVert_{\lambda}.$

In other words, $\beta\equiv\beta_{m,\lambda}=\beta(F_{X}|m,\lambda)$ is
chosen according to $\lambda$ to minimize the “distance” between the true
model $m^{*}$ and the working model $m$ according to a given loss-function
$||\cdot||_{\lambda}$. When $Y$ is part of an exponential family, a natural
choice for the loss function is derived from its log-likelihood (for more
details see appendix). Note that the trajectory groups $C$ are only a function
of the time-varying treatment ($\bar{A}$) in the LCGM approach. Therefore, we
can write $m(C,V|\beta)\equiv m(\bar{a},V|\beta)$. For example, when the
outcome $Y$ is continuous, a weighted $L_{2}$ loss function is commonly used:

$\displaystyle\beta(\cdot|m,\lambda):\mathcal{M}\rightarrow\mathcal{R}^{p}$
$\displaystyle\beta(F_{X}|m,\lambda)=\underset{\beta\in\mathbf{R}^{k}}{argmin}\text{
}E_{F_{X}}\left[\underset{\bar{a}\in\mathcal{A}}{\sum}\left(Y^{\bar{a}}-m(\bar{a},V|\beta)\right)^{2}\lambda(\bar{a},V)\right].$

The choice of $\lambda(\bar{a},V)$ depends on the modeling goals. Common
choices notably include $\lambda(\bar{a},V)=1$, which would be suitable if it
is desired to get a global summary of the nonparametric causal model, and
$\lambda(\bar{a},V)=\prod_{t=1}^{K}P(A_{t}|\bar{A}_{t-1},V)$ which gives more
weights to trajectories that are frequently observed.

Case of a survival outcome
When $Y$ is a time-to-event outcome, we define our working causal model as a
Cox proportional hazard MSM where $h_{0}(t)$ is the baseline hazard function:

$h(t|\bar{a},V,\beta)=h_{0}(t)exp(\beta_{1}C_{1}+\ldots+\beta_{J-1}C_{J-1}).$
(6)

Indeed,
$E(T^{\bar{a}})=\int_{t}S(t)=\int_{t}exp(-H(t))=\int_{t}exp(-\int_{k}^{t}h(t|\bar{a},V,\beta))\equiv
m(\bar{a},V|\beta)$. As in our illustration, the loss function can be defined
from the negative partial log-likelihood. Let $Y^{\bar{a}}$ denote the minimum
between the counterfactual survival time $T^{\bar{a}}$ and the censoring time
denoted by $U^{\bar{a}}$. The indicator for the event is designated by
$\delta$ where $\delta=1$ if $T^{\bar{a}}<U^{\bar{a}}$ and 0 otherwise and
$\mathcal{RS}(l)=\\{i=1,\ldots n:T^{\bar{a}}>l\\}$ designate the set of
individuals at risk at time $l$. We defined the partial likelihood as a
function of the groups as follows:

$\displaystyle L(\beta)$
$\displaystyle=\prod_{i=1}^{n}\left\\{\dfrac{h(t_{i}|\bar{a}_{i},V_{i},\beta)}{\sum_{l\in\mathcal{RS}(l)}h(t_{l}|\bar{a}_{l},V_{l},\beta)}\right\\}^{\delta_{i}}$
$\displaystyle=\prod_{i=1}^{n}\left\\{\dfrac{exp(C_{i}\beta)}{\sum_{l\in\mathcal{RS}(l)}exp(C_{l}\beta)}\right\\}^{\delta_{i}}$

Thus, the partial log-likelihood is given by:

$\displaystyle log(L(\beta))$
$\displaystyle=\sum_{i=1}^{n}\delta_{i}\left\\{log(h(t_{i}|\bar{a}_{i},V_{i},\beta))-log(\sum_{l\in\mathcal{RS}(l)}h(t_{l}|\bar{a}_{l},V_{l},\beta))\right\\}$
$\displaystyle=\sum_{i=1}^{n}I(T_{i}^{\bar{a}}<U_{i}^{\bar{a}})log\left\\{\sum_{l\in\mathcal{RS}(l)}exp(C_{i}\beta-
C_{l}\beta)\right\\}$

We define our parameter of interest $\beta$ as follows:

$\displaystyle\beta(\cdot|m,\lambda):\mathcal{M}\rightarrow\mathcal{R}^{p}$
$\displaystyle\beta(F_{X}|m,\lambda)=\underset{\beta\in\mathbf{R}^{k}}{argmin}\text{
}E_{F_{X}}\left\\{\bm{-}\sum_{\bar{a}\in\mathcal{\bar{A}}}\left[log(L(\beta))\right]\lambda(\bar{a},V)\right\\}.$

##  5 Estimation of the parameter of the LCGM-MSM

We now describe how the parameter $\beta$ of our MSM-LCGM can be estimated
from the observed data. Because, the distribution of the trajectory groups $C$
does not depend on the parameter of interest $\beta$, $C$ is an ancillary
statistic ([2, 9]). Therefore, in the second-step of our approach, the
trajectory groups are not seen as random but rather as a fixed regressor
([8]). In other words, the fact that the trajectories are estimated from the
data can be ignored when estimating the parameter of our LCGM-MSM and when
estimating the standard error of the estimator and confidence intervals.
However, it should be noted that inferences are conditional on the selected
LCGM.

Inverse probability of treatment weighting

Several approaches exist to estimate the parameters of a MSM and one of the
most popular is the inverse probability of treatment weighting (IPTW). The
IPTW mimics a random assignment of the treatment by creating a pseudo-
population where the treated and the untreated groups are comparable. In
longitudinal settings, IPTW can appropriately adjust for time-varying
covariates affected by prior exposure and selection bias. Following [15], the
IPTW is given by:

$W(\bar{A},V)=\dfrac{\lambda(\bar{A},V)}{g(\bar{A}|X)}$ (7)

where $g(\bar{A}|X)$ is the treatment mechanism which may depend on the
counterfactual data $X$. However, under the sequential conditional
exchangeability assumption,
$g(\bar{A}|X)=\prod_{t=1}^{K}P(A_{t}=1|\bar{A}_{t-1},\bar{L}_{t})$. In the
literature, stabilized weights (i.e when we set
$\lambda(\bar{A},V)=\prod_{t=1}^{K}P(\bar{A}_{t}|\bar{A}_{t-1},V))$ are
largely recommended since they tend to yield less extreme values ([19]). In
other words, each subject receives a weight that is proportional to the
product of the inverse of the probability of their observed treatment at each
time-point conditional on previous confounders and treatments. Note that the
weights are thus computed based on the observed treatment trajectories
($\bar{A}$) and not on the trajectory groups $C$.

IPTW estimating function

To estimate the parameter of interest $\beta_{m,\lambda}$, we have to define
the estimating function based on the IPTW. We first differentiate the loss
function and then we derive from its expression the IPTW estimating function.
We denote by $D(X|\beta)$ its first order derivative. The IPTW estimating
function when $Y$ is continuous is given by:

$D(\mathcal{O}|\beta,g)=\dfrac{D(X|\beta)}{g(\bar{A}|X)}$ (8)

where the first order derivative of the $L_{2}$ loss-function is given by:

$\displaystyle
D(X|\beta)=\underset{\bar{a}\in\mathcal{A}}{\sum}\lambda(\bar{a},V)\dfrac{\partial}{\partial\beta}m(\bar{a},V|\beta)\epsilon_{\bar{a}}(\beta)$

with $\epsilon_{\bar{a}}(\beta)=Y^{\bar{a}}-m(\bar{a},V|\beta)$ is the
residual after projecting the nonparametric MSM onto the chosen working model.
A well known property of the IPTW estimating function is its unbiasness at its
true parameter which satisfy $E(D(\mathcal{O}|g,\beta))=0$ (see for example
[15]). Thus, it can be shown under the positivity assumption, that the IPTW
estimator $\hat{\beta}_{m_{IPTW},\lambda}$ of $\beta_{m,\lambda}$ is the
solution of the estimating equation
$\sum_{i=1}^{n}D(\mathcal{O}_{i}|g_{n},\beta)=0$ where $\lambda_{n}$ and
$g_{n}$ are estimators of $\lambda$ and $g$. To simplify the notation, we
define $\hat{\beta}_{m_{IPTW},\lambda}\equiv\hat{\beta}_{n}$. Under the
positivity assumption, when $g_{n}$ and $\lambda_{n}$ are consistent,
$\hat{\beta}_{n}$ is consistent and asymptotically linear ([15]). An important
property of an asymptotically linear estimator is that it has a unique
influence function ([21]). In other words, we identify $\hat{\beta}_{n}$
through its influence function. From the expression of the influence function
we determine the asymptotic variance of $\hat{\beta}_{n}$ using Slutsky’s
theorem. Thus, we can construct the 95% confidence interval for
$\hat{\beta}_{n}$ given by:

$\hat{\beta}_{jn}\pm 1.96\times se(\hat{\beta}_{jn}),$ (9)

where $se(\hat{\beta}_{jn})=\hat{\sum}_{jn}^{p\times p}$ denotes the $j^{th}$
diagonal element of the sandwich variance estimator matrix of dimension
$p\times p$ where $p$ number of parameters in $m(\bar{a},V|\beta)$. Figure (2)
provides a practical template of the estimation approach. In practice, when
the distribution of $Y$ is part of the exponential family, the parameters of
the working model can be directly estimated using the classical weighted
generalized estimating equations (GEE), with, for example, the function geeglm
from the R package geepack. Indeed, under the causal assumptions described in
Section 2, the estimating Equation (8) are solved by such a weighted GEE
procedure. When $Y$ is a time-to-event outcome, a weighted Cox model can be
fitted using the R package survival. A simple practical solution for inference
is to use the sandwich variance estimator of the GEE routine or of the Cox
model, which treats the weights as known and results in conservative
statistical inferences ([6]).

Figure 2: Practical template for implementation of the LCGM-MSM approach

##  6 Simulation study

We aimed to evaluate the performance of our proposed approach for different
types of outcomes, number of trajectory groups and follow-up times. We also
wished to illustrate how alternative methods that make use of LCGM for causal
inference can yield biased results in presence of time-dependent confounding.
We have considered scenarios where: Y was a continuous, binary or time-to-
event outcome; treatment and covariable were binary and time-varying. We chose
$K=3,5$ and $10$ total number of follow-up times which yield respectively
$8,32$ and $1024$ possible treatment trajectories. These trajectories were
summarized using LCGM in $J=$ 3, 4 and 5 trajectory groups. For all scenarios,
we considered that the exposure and the covariate follow a binomial
distribution (Equations are presented in the supplementary material, appendix
B).

To estimate the effect of the trajectory groups, we fitted a LCGM where the
exposure $A$ is a linear function of the time on the logit scale. For each
outcome, we estimated the working model $m(\bar{a},V|\beta)$ using either a
weighted GEE regression or a weighted Cox regression with the inverse of
probability of treatment as weights. In the case where $Y$ is continuous, we
contrast our approach with a naive one where we adjusted for baseline and
time-varying covariate when $K=3,5$. This was not done for the case where $Y$
is binary or time-to-event outcome because of the problem of non-
collapsibility. We also implemented an approach where after modeling the LCGM,
probabilities of belonging in each trajectory groups are estimated by fitting
a multinomial model with the trajectory group as dependent variable and the
baseline characteristics of individuals as independents variables. In
particular, [11] estimated a MSM with the stabilized inverse of these
probabilities as weights. They stabilized the weights with the marginal
probability of the trajectory groups. We henceforth denote this approach as
the inverse probability of trajectory groups weighting (IPTGW). We also
applied this approach in the case where $Y$ was continuous. To measure the
performance of the LCGM-MSM and alternatives approachs, we determined the true
parameter using counterfactual data. Analytically determining the true value
of the parameter of interest $\beta$ is challenging, particularly since this
parameter is dependent on the chosen LCGM. The parameter of interest is thus
susceptible to change from one simulation replicate to another. Instead, we
determined the true value using a Monte Carlo simulation of counterfactual
data (for more details see supplementary material, appendix B).

##  7 Results

Table 1 presents the LCGM-MSM and crude results for the scenarios with three
trajectory groups. Similar results were obtained in the scenarios with four or
five trajectory groups (see Web Table 1 and 2 in the online supplementary
material, appendix C). The bias for LCGM-MSM when using either unstabilized or
stabilized weights is small in all considered scenarios. However, stabilized
weights yield estimators with a smaller standard deviation. The confidence
intervals also tend to have better coverage when using stabilized weights in
scenarios with multiple time points. Indeed, all confidence interval coverages
were closed to 95% (between 90% and 98%). This is likely because stabilization
tend to avoid extreme outlying weight values. On the other hand, the crude
model estimator has a large bias and low coverage probabilities of confidence
intervals. Table 2 presents the results for the baseline adjusted, time-
varying adjusted and IPTGW adjusted LCGM. Recall that only the continuous
outcome case, with either $K=3$ or $K=5$ time-points were considered for these
alternatives methods. Moreover, regardless of the number of follow-up times
and number of trajectory groups, alternative LCGM analyses were highly biased
with low coverage of their confidence interval (between 2% and 67%).

Table 1: LCGM-MSM with unstabilized and stabilized IPTW estimator when $Y$ is
continuous, binary and time-to-eventt outcome for three trajectories.

| Unstabilized IPTW | Stabilized IPTW
---|---|---
$Y$ continous
t=3
| Bias | SDE | CP | Crude bias | Crude CP |  | Bias | SDE | CP | Crude bias | Crude CP
Intercept | 0.010 | 0.200 | 95% | 0.040 | 6% |  | 0.000 | 0.100 | 95% | 0.150 | 11%
Group 2 | 0.000 | 0.480 | 95% | -0.200 | 16% |  | 0.000 | 0.370 | 95% | -0.250 | 31%
Group 3 | 0.000 | 0.200 | 96% | -0.340 | 8% |  | 0.000 | 0.160 | 97% | -0.370 | 0%
t=5
| Bias | SDE | CP | Crude bias | Crude CP |  | Bias | SDE | CP | Crude bias | Crude CP
Intercept | -0.010 | 0.370 | 92% | -0.090 | 36% |  | -0.010 | 0.210 | 94% | 0.170 | 17%
Group 2 | 0.000 | 0.770 | 93% | -0.050 | 19% |  | 0.000 | 0.520 | 95% | -0.120 | 50%
Group 3 | 0.010 | 0.360 | 94% | -0.420 | 12% |  | 0.030 | 0.230 | 95% | -0.350 | 1%
t=10
| Bias | SDE | CP | Crude bias | Crude CP |  | Bias | SDE | CP | Crude bias | Crude CP
Intercept | -0.024 | 0.308 | 87% | -0.132 | 39% |  | -0.036 | 0.200 | 86% | 0.130 | 28%
Group 2 | -0.008 | 0.638 | 90% | -0.287 | 25% |  | 0.027 | 0.431 | 92% | -0.137 | 50%
Group 3 | -0.037 | 0.375 | 88% | -0.497 | 5% |  | 0.026 | 0.268 | 91% | -0.352 | 1%
$Y$ binary
t=3
| Bias | SDE | CP | Crude bias | Crude CP |  | Bias | SDE | CP | Crude bias | Crude CP
Intercept | 0.002 | 0.195 | 95% | 0.033 | 50% |  | -0.005 | 0.112 | 95% | 0.140 | 48%
Group 2 | 0.003 | 0.462 | 96% | -0.177 | 35% |  | 0.010 | 0.354 | 95% | -0.227 | 43%
Group 3 | 0.001 | 0.218 | 96% | -0.299 | 31% |  | 0.016 | 0.173 | 96% | -0.320 | 21%
t=5
| Bias | SDE | CP | Crude bias | Crude CP |  | Bias | SDE | CP | Crude bias | Crude CP
Intercept | -0.012 | 0.376 | 93% | -0.091 | 54% |  | -0.010 | 0.250 | 94% | 0.144 | 64%
Group 2 | 0.005 | 0.761 | 93% | -0.048 | 42% |  | 0.008 | 0.557 | 94% | -0.108 | 69%
Group 3 | 0.009 | 0.389 | 94% | -0.344 | 33% |  | 0.020 | 0.275 | 94% | -0.295 | 33%
t=10
| Bias | SDE | CP | Crude bias | Crude CP |  | Bias | SDE | CP | Crude bias | Crude CP
Intercept | -0.021 | 0.353 | 91% | -0.125 | 62% |  | -0.036 | 0.218 | 91% | 0.115 | 66%
Group 2 | 0.022 | 0.630 | 92% | -0.221 | 42% |  | 0.033 | 0.417 | 95% | -0.106 | 79%
Group 3 | -0.011 | 0.430 | 90% | -0.360 | 19% |  | 0.029 | 0.328 | 93% | -0.270 | 35%
$Y$ time-to-event
t=3
| Bias | SDE | CP | Crude bias | Crude CP |  | Bias | SDE | CP | Crude bias | Crude CP
Group 2 | 0.000 | 0.33 | 96% | -0.111 | 25% |  | 0.003 | 0.274 | 95% | -0.155 | 37%
Group 3 | 0.003 | 0.14 | 96% | -0.186 | 25% |  | 0.002 | 0.140 | 95% | -0.240 | 4%
t=5
| Bias | SDE | CP | Crude bias | Crude CP |  | Bias | SDE | CP | Crude bias | Crude CP
Group 2 | 0.008 | 0.49 | 94% | -0.012 | 35% |  | 0.006 | 0.353 | 94% | -0.075 | 65%
Group 3 | 0.001 | 0.22 | 94% | -0.165 | 35% |  | 0.018 | 0.163 | 94% | -0.190 | 13%
t=10
| Bias | SDE | CP | Crude bias | Crude CP |  | Bias | SDE | CP | Crude bias | Crude CP
Group 2 | -0.004 | 0.34 | 91% | -0.100 | 40% |  | 0.020 | 0.284 | 96% | -0.069 | 85%
Group 3 | -0.016 | 0.18 | 91% | -0.158 | 35% |  | 0.029 | 0.181 | 91% | -0.166 | 52%

Legend: bias, $CP$ and $SDE$ indicate respectively the bias, the coverage
probability and the standard deviation of the estimate of $\beta_{IPTW}$ for
the adjusted model with IPTW weight, crude %bias and crude $CP$ indicate
respectively the bias and the coverage probability of the estimate of
$\beta_{crude}$ from the unadjusted model.

Table 2: LCGM-MSM with baseline adjustment, time varying covariates adjustment
and IPTGW as an estimation approach of a MSM when $Y$ is continuous.

| baseline adjustment |  | tvc adjustment |  | Stabilized IPTGW adjustment
---|---|---|---|---|---
3 trajectories
K=3
| bias | SDE | CP |  | bias | SDE | CP |  | bias | SDE | CP
Intercept | -0.129 | 0.07 | 59% |  | -0.721 | 0.13 | 0% |  | 0.130 | 0.056 | 24%
Group 2 | -0.182 | 0.17 | 8% |  | -0.062 | 0.43 | 57% |  | -0.211 | 0.123 | 33%
Group 3 | -0.284 | 0.12 | 18% |  | -0.255 | 0.14 | 15% |  | -0.310 | 0.136 | 2%
K=5
| bias | SDE | CP |  | bias | SDE | CP |  | bias | SDE | CP
Intercept | -0.221 | 0.15 | 29% |  | -0.955 | 0.20 | 0% |  | 0.146 | 0.100 | 30%
Group 2 | -0.062 | 0.33 | 15% |  | -0.122 | 0.51 | 27% |  | -0.091 | 0.200 | 46%
Group 3 | -0.398 | 0.15 | 17% |  | -0.517 | 0.18 | 2% |  | -0.305 | 0.130 | 4%
4 trajectories
K=3
| bias | SDE | CP |  | bias | SDE | CP |  | bias | SDE | CP
Intercept | -0.087 | 0.09 | 69% |  | -0.680 | 0.13 | 0% |  | 0.140 | 0.083 | 17%
Group 2 | -0.109 | 0.15 | 31% |  | -0.044 | 0.38 | 60% |  | -0.177 | 0.228 | 41%
Group 3 | -0.288 | 0.18 | 22% |  | -0.182 | 0.34 | 35% |  | -0.279 | 0.181 | 25%
Group 4 | -0.303 | 0.14 | 20% |  | -0.300 | 0.25 | 16% |  | -0.311 | 0.166 | 10%
K=5
| bias | SDE | CP |  | bias | SDE | CP |  | bias | SDE | CP
Intercept | -0.193 | 0.15 | 26% |  | -0.913 | 0.18 | 0% |  | 0.152 | 0.110 | 25%
Group 2 | 0.112 | 0.28 | 25% |  | -0.014 | 0.43 | 51% |  | -0.031 | 0.207 | 60%
Group 3 | -0.248 | 0.24 | 27% |  | -0.333 | 0.34 | 19% |  | -0.243 | 0.197 | 27%
Group 4 | -0.410 | 0.22 | 13% |  | -0.527 | 0.27 | 4% |  | -0.349 | 0.173 | 5%
5 trajectories
K=3
| bias | SDE | CP |  | bias | SDE | CP |  | bias | SDE | CP
Intercept | -0.051 | 0.09 | 77% |  | -0.648 | 0.11 | 0% |  | 0.153 | 0.083 | 11%
Group 2 | -0.162 | 0.14 | 42% |  | -0.104 | 0.34 | 54% |  | -0.185 | 0.235 | 42%
Group 3 | -0.169 | 0.17 | 37% |  | -0.108 | 0.38 | 45% |  | -0.196 | 0.229 | 40%
Group 4 | -0.329 | 0.19 | 23% |  | -0.277 | 0.36 | 20% |  | -0.321 | 0.169 | 25%
Group 5 | -0.335 | 0.13 | 19% |  | -0.320 | 0.24 | 18% |  | -0.355 | 0.154 | 7%
K=5
| bias | SDE | CP |  | bias | SDE | CP |  |  |  |
Intercept | -0.178 | 0.15 | 28% |  | -0.892 | 0.18 | 0% |  | 0.159 | 0.110 | 21%
Group 2 | 0.127 | 0.25 | 39% |  | -0.002 | 0.36 | 57% |  | -0.031 | 0.223 | 67%
Group 3 | -0.104 | 0.28 | 41% |  | -0.154 | 0.39 | 42% |  | -0.157 | 0.243 | 50%
Group 4 | -0.264 | 0.23 | 31% |  | -0.368 | 0.34 | 19% |  | -0.280 | 0.211 | 27%
Group 5 | -0.424 | 0.29 | 14% |  | -0.555 | 0.31 | 5% |  | -0.390 | 0.185 | 4%

Legend: bias, $CP$ and $SDE$ indicate respectively the bias, the coverage
probability and the standard deviation of the estimate of $\beta$ for each
adjustment approach. Abbreviations : tvc=time-varying covariates, IPTGW=
inverse of probability of trajectory groups weights.

##  8 Discussion

We proposed a new approach combining LCGM and MSMs. In a first step, we
propose to use LCGM to classify individuals into a few latent classes based on
their treatment adherence pattern, then in a second step, choose a working MSM
that relates the outcome to these latent classes. The parameter of interest is
nonparametrically defined as the projection of the true MSM onto the chosen
working causal model. This approach has a double benefits. On one hand, we
avoid assuming that we know the functional form of the true model as in a
parametric model. On another hand, we also avoid the curse of dimensionality
that is often encountered when estimating a nonparametric model with finite
sample data. We showed that, the data-driven estimation of the trajectory
groups can be ignored. As such, parameters can be estimated using the IPTW
estimator and conservative inferences can be obtained using a standard robust
variance estimator. Simulation studies were used to illustrate our approach
and compare it with unadjusted, baseline covariates-adjusted, time-varying
covariates adjusted and inverse probability of trajectory groups weighting
adjusted alternatives. Scenarios with varying types of outcomes (continuous,
binary or time-to-event) and number of follow-up times (3, 5, 10) and number
of trajectory classes (3, 4, 5) were considered. In the simulation study, we
found that our proposed approach yielded estimators with little or no bias.
Alternative analyses were biased and had low coverage of their confidence
intervals.

Despite the strengths of our proposal, some limitations need to considered.
For example, in presence of time-varying confounding, the local independence
assumption of the LCGM is unlikely to hold. Indeed, when treatment uses at
different time-points are locally dependent, overlapping or redundant
information is accounted for when it should not be in the estimation of the
LCGM ([22]). Hence, omitting time-varying covariates when modeling trajectory
groups might lead to an overestimation of the “posterior” probability. Also,
we chose arbitrarily the number of trajectory groups in the simulation study
and in the application. As [13] mentioned the groups can be seen as points of
support on the unknown distribution of trajectories. In that sense, the groups
are not meant to represent the true data-generating process. In the case when
$Y$ is binary, one might be cautious regarding the number of available events
per trajectory group to avoid estimation problems.

Our proposed approach can be improved by using a doubly robust estimator
instead of the IPTW estimator. It is clear in our simulation results that as
the number of follow-up time increased, the IPTW estimator failed to
completely eliminate the bias. Indeed, the IPTW estimator heavily relies on a
correct specification of the probability of treatment assignment and is
sensitive to near violations of the positivity assumptions. Furthermore,
relaxing the assumption of local independence might also improve the LCGM-MSM
through a better classification of individuals’ treatment patterns. The number
of trajectory groups can be chose using the usual recommendations in the
literature such as the BIC criterion or by cross-validation ([13, 16]). In
conclusion, we believe that our LCGM-MSM can be useful to correctly estimate
the effect of trajectory groups on an outcome when conducting longitudinal
studies based on administrative data. Our proposed approach might help
practitioners to have a better idea of what type of adherence profile is worth
targeting as LCGM is better to measure treatment adherence pattern than most
other commonly used methods.

## Acknowledgments

This project was funded by a CIHR grant [$\\#$ 165942] and a start-up award
from Faculté de Médecine de l’Université Laval to DT. CS, JRG and DT are
supported by career awards from the Fonds de recherche du Québec – Santé. AD
received a scholarship from VITAM.”

## References

  * Armitage et al. [2019] Jane Armitage, Colin Baigent, Elizabeth Barnes, D John Betteridge, Lisa Blackwell, Michael Blazing, Louise Bowman, Eugene Braunwald, Robert Byington, Christopher Cannon, et al.  Efficacy and safety of statin therapy in older people: a meta-analysis of individual participant data from 28 randomised controlled trials.  _The Lancet_ , 393(10170):407–415, 2019.
  * Basu [1959] D Basu.  The family of ancillary statistics.  _Sankhyā: The Indian Journal of Statistics_ , pages 247–256, 1959\.
  * Chowdhury et al. [2013] Rajiv Chowdhury, Hassan Khan, Emma Heydon, Amir Shroufi, Saman Fahimi, Carmel Moore, Bruno Stricker, Shanthi Mendis, Albert Hofman, Jonathan Mant, et al.  Adherence to cardiovascular therapy: a meta-analysis of prevalence and clinical consequences.  _European heart journal_ , 34(38):2940–2948, 2013\.
  * Franklin et al. [2013] Jessica M Franklin, William H Shrank, Juliana Pakes, Gabriel Sanfélix-Gimeno, Olga S Matlin, Troyen A Brennan, and Niteesh K Choudhry.  Group-based trajectory models: a new approach to classifying and predicting long-term medication adherence.  _Medical care_ , 51(9):789–796, 2013.
  * Franklin et al. [2015] Jessica M Franklin, Alexis A Krumme, Angela Y Tong, William H Shrank, Olga S Matlin, Troyen A Brennan, and Niteesh K Choudhry.  Association between trajectories of statin adherence and subsequent cardiovascular events.  _pharmacoepidemiology and drug safety_ , 24(10):1105–1113, 2015.
  * Hernan and Robins [2020] Miguel A Hernan and James M Robins.  _Causal Inference: What If_.  Chapman and Hall/CRC Boca Raton, FL, 2020.
  * Hernán et al. [2000] Miguel Ángel Hernán, Babette Brumback, and James M Robins.  Marginal structural models to estimate the causal effect of zidovudine on the survival of hiv-positive men.  _Epidemiology_ , pages 561–570, 2000.
  * Lehmann and Romano [2006] Erich L Lehmann and Joseph P Romano.  _Testing statistical hypotheses_.  Springer Science & Business Media, 2006.
  * Lehmann and Scholz [1992] Erich Leo Lehmann and FW Scholz.  Ancillarity.  _Lecture Notes-Monograph Series_ , 17:32–51, 1992.
  * Leisch [2004] Friedrich Leisch.  Flexmix: A general framework for finite mixture models and latent glass regression in r.  2004\.
  * Lim et al. [2015] Sungwoo Lim, Tiffany G Harris, Denis Nash, Mary Clare Lennon, and Lorna E Thorpe.  All-cause, drug-related, and hiv-related mortality risk by trajectories of jail incarceration and homelessness among adults in new york city.  _American journal of epidemiology_ , 181(4):261–270, 2015.
  * Muthen [2001] Bengt Muthen.  Latent variable mixture modeling.  _New developments and techniques in structural equation modeling_ , 2:1–33, 2001.
  * Nagin [2005] Daniel S Nagin.  _Group-based modeling of development_.  Harvard University Press, 2005.
  * Neugebauer and van der Laan [2007] Romain Neugebauer and Mark van der Laan.  Nonparametric causal effects based on marginal structural models.  _Journal of Statistical Planning and Inference_ , 137(2):419–434, 2007.
  * Neugebauer and van der Laan [2003] Romain Neugebauer and Mark J van der Laan.  Locally efficient estimation of nonparametric causal effects on mean outcomes in longitudinal studies.  2003\.
  * Nielsen et al. [2014] Jason D Nielsen, Jeffry S Rosenthal, Ye Sun, David M Day, Irene Bevc, and Thierry Duchesne.  Group-based criminal trajectory analysis using cross-validation criteria.  _Communications in Statistics-Theory and Methods_ , 43(20):4337–4356, 2014.
  * Robins [2000] James M Robins.  _Marginal structural models versus structural nested models as tools for causal inference_.  Springer, 2000.
  * Singh et al. [2018] Sonal Singh, Susan Zieman, Alan S Go, Stephen P Fortmann, Nanette K Wenger, Jerome L Fleg, Barbara Radziszewska, Neil J Stone, Sophia Zoungas, and Jerry H Gurwitz.  Statins for primary prevention in older adults—moving toward evidence-based decision-making.  _Journal of the American Geriatrics Society_ , 66(11):2188–2196, 2018.
  * Talbot et al. [2015] Denis Talbot, Juli Atherton, Amanda M Rossi, Simon L Bacon, and Geneviève Lefebvre.  A cautionary note concerning the use of stabilized weights in marginal structural models.  _Statistics in medicine_ , 34(5):812–823, 2015\.
  * Taylor et al. [2013] F Taylor, MD Huffman, AF Macedo, THM Moore, M Burke, Davey Smith G Ward, and S K Ebrahim.  Statins for the primary prevention of cardiovascular disease.  _Cochrane Database of Systematic Reviews_ , (1), 2013.  ISSN 1465-1858.  doi: 10.1002/14651858.CD004816.pub5.  URL https://doi.org//10.1002/14651858.CD004816.pub5.
  * Tsiatis [2007] Anastasios Tsiatis.  _Semiparametric theory and missing data_.  Springer Science & Business Media, 2007.
  * Vermunt and Magidson [2002] Jeroen K Vermunt and Jay Magidson.  Latent class cluster analysis.  _Applied latent class analysis_ , 11(89-106):60, 2002.
  * Vermunt and Magidson [2004] Jeroen K Vermunt and Jay Magidson.  Local independence.  _Encyclopedia of social sciences research methods_ , pages 732–733, 2004.
  * WHO [2019] WHO.  Cardiovascular diseases (cvds), 2019.  URL https://www.who.int/news-room/fact-sheets/detail/the-top-10-causes-of-death.
  * Xiao et al. [2014] Yongling Xiao, Michal Abrahamowicz, Erica EM Moodie, Rainer Weber, and James Young.  Flexible marginal structural models for estimating the cumulative effect of a time-dependent treatment on the hazard: reassessing the cardiovascular risks of didanosine treatment in the swiss hiv cohort study.  _Journal of the American Statistical Association_ , 109(506):455–464, 2014.