Modeling treatment effects on absorbing outcomes in clinical trials: Leveraging longitudinal and ordinal data for efficiency gains

Abstract

Randomized clinical trials in critically ill populations often use mortality at the end of follow-up as the primary outcome. When mortality is rare, composite outcomes that capture disease severity over time are frequently employed to improve efficiency. Motivated by the Crystalloid Liberal or Vasopressor Early Resuscitation in Sepsis (CLOVERS) trial—a randomized study comparing two fluid resuscitation strategies in hospitalized patients with sepsis—we evaluate multiple modeling approaches for estimating covariate-adjusted treatment effects on absorbing outcomes such as mortality. While analyses of binary end-of-follow-up outcomes (alive vs. dead) are natural, efficiency can be gained by incorporating low-cost longitudinal and ordinal severity data routinely collected in trials and by imposing reasonable modeling restrictions (constraints). We develop methods that leverage these data within the marginalized model framework, focusing on first-order marginalized transition models, which allows direct estimation of outcome-marginal (population-averaged) treatment effects. We highlight a key property of the marginalized transition model in the presence of absorbing states: the marginal and conditional sub-models correspond to cumulative incidence and discrete hazard models, respectively. Through simulation studies and an application to the CLOVERS data, we compare the validity and efficiency of competing approaches, including standard survival analyses.

Keywords

longitudinal data ordinal responses randomized clinical trial absorbing state marginalized models

1. Introduction

We seek to summarize the treatment effect on disease outcomes in critically ill participants in a randomized clinical trial. A relevant outcome for such studies is mortality at the end of a follow-up period. However, mortality is often sufficiently rare that standard analyses would require excessively large sample sizes or would otherwise lack power. To address concerns about power and efficiency, composite outcomes, such as ventilator-free or oxygen-free days, are often used. These outcomes summarize disease severity over the course of follow-up and have been implemented in multiple studies.^1,2 They are often modeled with ordinal regression methods with mortality being the worst/most severe outcome category. While composite outcomes can improve power to detect associations^3,4 and are pragmatic for studies of critically ill participants, they have limitations. For example, when the treatment effect on mortality is different from the treatment effect on non-fatal severity, interpreting the model’s estimates can be difficult, and power can be reduced even when compared to standard cross-sectional analyses of mortality.⁵

Historically, analyses of baseline covariate-adjusted treatment effects in randomized clinical trials have been conducted with binary logistic regression or with survival analysis models. Binary logistic regression has low power if the prevalence of mortality is low, and a survival analysis can improve power in the presence of dropout. However, a covariate-adjusted hazard ratio may not be of interest because it is based on analyses of a risk set that changes over time and there is an implicit conditioning on a post-randomization, ‘intermediate’ outcome (i.e. mortality has not yet occurred).^6–9 Our approach to improve efficiency over a cross-sectional, binary logistic regression analysis does so by exploiting data that are commonly collected in randomized clinical trials or that could be collected at relatively low costs. Rather than modeling a cross-sectional, end of study, binary response, we can enrich analyses by leveraging daily (longitudinal), severity (ordinal) response data.^10–13 We may then impose constraints or restrictions on the regression model to directly improve efficiency. From an ordinal response model perspective, we may impose proportional odds or partial proportional odds restrictions to improve precision over a more highly parameterized, general cumulative logit regression model. From a longitudinal data perspective, we may impose restrictions on the treatment effect over time to improve estimation efficiency. Notably, like the Cox proportional hazard model,¹⁴ longitudinal data analysis can improve efficiency and robustness of a treatment effect estimator by using the timing of events and the partial information contained in participants who drop out or are lost to follow-up.^15,16 Cross-sectional, binary logistic regression and ordinal regression analyses discard participant’s data prior to dropout.

In this paper, we compare approaches for estimating baseline covariate-adjusted treatment effects on mortality at the end of follow-up. In particular, we focus on longitudinal outcome-marginal models, which estimate covariate-adjusted treatment effects without conditioning on past or future outcome values as in done in standard transition models or mixed effects models.^17–20 Within the broad class of outcome-marginal models, we appeal to marginalized models. Marginalized models were developed for longitudinal binary^21–23 and ordinal outcomes^10,11,13 as likelihood-based alternatives to semiparametric generalized estimating equations²⁴ and as outcome-marginal model alternatives to generalized linear mixed models and transition models.^25,26 Briefly, marginalized models consists of two components: the outcome-marginal model linking the longitudinal outcome to the exposure and covariates and the conditional response, dependence model that captures the outcome dependence or correlation over-time. We will show that when the target outcome state is mortality, these two models correspond to other, more familiar models. Specifically, the outcome-marginal model corresponds to a cumulative incidence model while the conditional, response dependence model corresponds to a discrete hazard model. Thus, unlike in standard survival analysis, marginalized transition models estimate the covariate-adjusted, treatment effect on cumulative incidence of mortality directly. Similar to standard survival models, one can then obtain an end-of-study fully marginal (with respect to both outcome and baseline covariates) treatment effect by integrating over the covariates distribution.

1.1. Organization of the paper

Motivated by the Crystalloid Liberal or Vasopressor Early Resuscitation in Sepsis (CLOVERS) randomized clinical trial, we demonstrate the potential benefit of exploiting longitudinal (daily) and ordinal (health state) data to improve estimation efficiency when interest is in an end-of-follow-up treatment effect on an absorbing outcome state (e.g. death). We extend our previous work¹³ to show that marginalized transition models are well-suited for this setting. We show that the two components of a marginalized transition model naturally correspond to a cumulative incidence model and a discrete hazards model, allowing direct estimation of outcome-marginal, covariate-adjusted treatment effects. Via simulation, we also show the extent to which precision and robustness to missing data can be improved compared to common analyses (e.g. cross-sectional logistic regression) when the models are correctly specified, and the extent of bias incurred when assumption choices are overly restrictive.

In the next section, we describe the CLOVERS trial.² In Section 3, we (i) review the first-order marginalized transition model (MTM) for longitudinal binary response data, (ii) describe its outcome-marginal and conditional, response-dependence components, (iii) discuss restrictions imposed when one of the outcome states is absorbing (e.g. death), and (iv) explore the connection between the MTM components and the more familiar cumulative incidence and discrete hazard models. In Section 4, we show an MTM extension to longitudinal, ordinal outcomes and describe how the ordinal MTM can exploit both the longitudinal and ordinal dimensions of response data to increase efficiency of the outcome-marginal, covariate-adjusted treatment effect estimate on mortality at the end of follow-up. We then conduct a comparative simulation study to examine efficiency and robustness of the estimation strategies in Section 5. In Section 6, we turn back to the CLOVERS trial and conduct analyses using the described methods. We wrap up the paper with a discussion in Section 7.

2. CLOVERS study: data

The CLOVERS study was a multi-center randomized clinical trial on the efficacy of two fluid resuscitation strategies in patients with sepsis.² The trial randomly assigned 1,563 patients to restrictive or liberal fluid resuscitation strategies, and the primary analysis was a survival analysis on mortality before discharge within 90 days after randomization. CLOVERS also conducted numerous secondary analyses that included composite outcomes such as hospital-free days and ICU-free days through 28 days post-randomization.² To calculate hospital-free days and ICU-free days, the study collected daily patient state data during the 28 days post-randomization and summarized them to create the outcomes. We study the setting where interest is in a cross-sectional, outcome-marginal, covariate-adjusted treatment effect on mortality at day 28. To improve precision and power over commonly used analysis approaches, we use the available daily state data to estimate this quantity. Specifically, participants were observed daily to be in one of four ordered states from randomization to day 28: discharged/home, hospital, hospital/ICU, and death. We compare treatment effect estimation strategies using cross-sectional (at day 28) and longitudinal data analysis procedures for binary and ordinal outcomes. The binary outcome combines the states discharged/home, hospital and hospital/ICU into a single alive state while the ordinal outcome considers the four states separately.

We include 1,529 of the 1,563 participants in the CLOVERS study with available daily data. Briefly, the mean age was 59.4 years (standard deviation 15.9 years) and mean sequential organ failure assessment (SOFA) score, indicating level of organ dysfunction, was 3.39 (standard deviation 2.67). Figure 1 shows the proportion of participants in each outcome state from day 0 (randomization) to day 28. As is common in studies of critically ill participants, the proportion of participants who die and/or are discharged/at home increases with time, while the proportions of participants in the hospital and hospital/ICU states decreases. While we are interested in the outcome-marginal, baseline covariate-adjusted, cross-sectional treatment effect on day 28 mortality (yellow vs not yellow on day 28), we want to use all intermediate follow-up days and outcome severity states to improve precision.

Figure 1.

Proportion of patients in each of the four outcome states for the first 28 days of follow-up. Day 0 refers to the day of randomization.

Table 1 shows day-to-day transition probabilities. Death is an absorbing state which is central to our analytical considerations. Discharged/home is nearly absorbing in that several participants died after being in the discharged/home state, although no one was observed to transition from discharged/home to hospital or hospital/ICU. The lack of transitions from discharged/home to hospital or hospital/ICU was part of the study design, that is, study coordinators did not follow hospitalizations after discharge; however, for the purposes of our analyses, we will assume that none of those transitions occurred. While participants transitioned regularly between the hospital and hospital/ICU states, the overall probability of remaining in the same state from day to day was high.

Table 1.

Day to day state transition probabilities across the 1,529 CLOVERS study participants.

	Yesterday [ $Y_{i (j - 1)}$ ]
Today [ $Y_{i j}$ ]	Death	Hospital/ICU	Hospital	Discharged/home	Overall
Death	1.000	0.015	0.005	0.002	0.092
Hospital/ICU	0.000	0.790	0.151	0.000	0.091
Hospital	0.000	0.159	0.835	0.000	0.167
Discharged/Home	0.000	0.036	0.008	0.998	0.650

In studies like CLOVERS, the simplest way to estimate the treatment effect on mortality at the end of follow-up is to fit a binary logistic regression model with the outcome being mortality status (yes/no) at day 28. However, we may gain efficiency by considering alternative strategies that incorporate richer outcome information into the estimation procedure. For example, we may focus analyses exclusively on day 28 but replace the binary mortality status with a four-level ordinal severity outcome (death, hospital/ICU, hospital, discharged/home); alternatively, we may exploit the complete daily trajectory from baseline to day 28 of the binary mortality status (yes/no); finally, we may use daily, severity information to fit a longitudinal, ordinal MTM. One would expect that incorporating richer outcome information into the estimation model could result in efficiency gains. However, to realize the efficiency gains, parameter constraints are likely to be required. For instance, prior work has shown that constraining the treatment effect at baseline to be equal to zero in a randomized clinical trial can improve efficiency.^27,28 This paper aims to explore the potential for efficiency gains from these alternative modeling strategies compared to binary logistic regression.

Note: In studies similar to CLOVERS where interest is in estimating the treatment effect on mortality at day 28 of follow-up, Cox proportional hazard models or other parametric survival models represent valid alternatives to those described herein. Although comparing the efficiency of MTMs with Cox proportional hazard models is not the primary focus of this paper, for completeness, we include a Cox proportional hazard model in our simulation studies and data analysis.

3. Marginalized transition models

Let $i \in {1, \dots, N}$ denote subject, $t_{i j} \in {0, \dots, m_{i} - 1}$ denote follow-up day on visit $j$ (i.e. $t_{i j} = j - 1$ ), $Y Y_{i}$ be a $m_{i} -$ vector of binary response data with $Y_{i j} \in {1, 2}$ , $X_{i}$ be a binary baseline variable indicating the subject $i$ treatment assignment (i.e. $X_{i} = 0$ if subject $i$ is assigned to the control arm or $X_{i} = 1$ if subject $i$ is assigned to the treatment arm), and $Z Z_{i}$ be a vector of baseline covariates. An outcome-marginal probability model, is given by $π_{i j} = P (Y_{i j} = 1 ∣ X_{i}, Z Z_{i})$ . We assume this model is of logistic form so that

logit (π_{i j}) = \underset{α (t_{i j})}{\underset{⏟}{[α_{0} + α_{t} t_{i j}]}} + \underset{β_{x} (t_{i j})}{\underset{⏟}{[β_{x} + β_{x t} t_{i j}]}} X_{i} + Z Z_{i} β β_{z} .

(1)

The

α (t_{i j})

function captures the

t_{i j} -

specific log-odds that the response is at level

1

in the control arm (

X_{i} = 0

) when

Z Z_{i} = 0 0

. It can also be thought of as a

t_{i j} -

specific intercept. The outcome-marginal probability model describes the first moment of the multivariate distribution

[Y Y_{i} ∣ X_{i}, Z Z_{i}]

. To identify higher moments, the marginalized model class specifies a conditional, response dependence structure, and here we appeal to the first order transition model of Azzalini²¹ and Heagerty,²³

logit (μ_{i j}^{c}) = Δ_{i j} (X_{i}, Z Z_{i}) + γ (X_{i}, Z Z_{i}) I (Y_{i (j - 1)} = 1),

(2)

where

μ_{i j}^{c} = P (Y_{i j} = 1 ∣ Y_{i (j - 1)}, X_{i}, Z Z_{i})

is a (conditional) Markov transition probability, and

γ (X_{i}, Z Z_{i})

is a log-odds ratio parameter that describes the strength and nature of the relationship between

Y_{i j}

and

Y_{i (j - 1)}

. The intercept,

Δ_{i j} (X_{i}, Z Z_{i})

, is an implicit function that links the outcome-marginal probability model to the conditional, response dependence model. That is,

Δ_{i j} (X_{i}, Z Z_{i})

is the value that solves

π_{i j} = π_{i (j - 1)} \cdot μ_{i j}^{c, 1} + (1 - π_{i (j - 1)}) \cdot μ_{i j}^{c, 2},

(3)

where

μ_{i j}^{c, y} = P (Y_{i j} = 1 ∣ Y_{i (j - 1)} = y, X_{i}, Z Z_{i})

is a transition probability. Given the vector of parameters

θ θ = (α α, β β, γ γ) = (α_{0}, α_{t}, β_{x}, β_{x t}, β β_{z}, γ)

, and once the values of

Δ_{i j} (X_{i}, Z Z_{i})

have been calculated for all

(i, j)

, subject

i

’s contribution to the likelihood can be identified directly with

\begin{aligned} L_{i} (θ θ; Y Y_{i}, X_{i}, Z Z_{i}) & = \prod_{j = 1}^{m_{i}} L_{i j} (θ θ; Y_{i j}, Y_{i (j - 1)}, X_{i}, Z Z_{i}) \\ = π_{i 1}^{I (Y_{i 1} = 1)} \cdot (1 - π_{i 1})^{I (Y_{i 1} = 2)} \cdot \prod_{j = 2}^{m_{i}} (μ_{i j}^{c})^{I (Y_{i j} = 1)} \cdot (1 - μ_{i j}^{c})^{I (Y_{i j} = 2)} . \end{aligned}

The first term includes the outcome-marginal probability at the first observed time point

π_{i 1}

because the precise value of

Y_{i 0}

is either not available or is included in the linear predictor. The second term includes the conditional model probability,

μ_{i j}^{c}

that explicitly conditions on the lagged response value. Notably, we use

Y_{i j} = 2

to denote a non-event which is not typical notation for binary response models. However, the notation will be consistent with the ordinal response model notation discussed later. The model parameters are estimated using a Newton–Raphson algorithm which alternates between updating

θ θ

and updating

Δ_{i j} (X_{i}, Z Z_{i})

so that equation (3) holds. Score and information calculations for this model are detailed in Heagerty.²³

3.1. Notes about the marginalized transition model when

Y_{i j} = 1

is absorbing

The marginalized transition model can be used for the analysis of longitudinal binary outcome data even when outcome state $Y_{i j} = 1$ is absorbing, that is, one may enter but cannot exit state $1$ . By examining equations (1), (2), and (3), one can see that, in the presence of an absorbing state, $π_{i j}$ is a cumulative incidence, and $μ_{i j}^{c}$ is a discrete hazard model. To explain, when $Y_{i j} = 1$ is absorbing, then $μ_{i j}^{c, 1} \equiv P (Y_{i j} = 1 ∣ X_{i}, Z_{i}, Y_{i (j - 1)} = 1)$ is equal to one and $π_{i j}$ is monotonically increasing with $t_{i j}$ because

π_{i j} = π_{i (j - 1)} \cdot μ_{i j}^{c, 1} + (1 - π_{i (j - 1)}) \cdot μ_{i j}^{c, 2} = π_{i (j - 1)} + (1 - π_{i (j - 1)}) \cdot μ_{i j}^{c, 2} \geq π_{i (j - 1)} .

Further,

π_{i j} = 1 - (1 - π_{i j}) = 1 - P (Y_{i j} = Y_{i (j - 1)} = \dots = Y_{i 1} = 2 ∣ X_{i}, Z Z_{i}, Y_{i 0} = 2)

which is a discrete cumulative incidence function or a failure cumulative distribution function. In contrast, survival analysis methods use post-estimation calculations to obtain cumulative incidence. The cumulative incidence calculation that, in standard survival analyses is done after the hazard model is fitted, is conducted within the marginalized transition model fitting routine, at each step of the Newton–Raphson procedure.

Additionally, because when $Y_{i j} = 1$ , $μ_{i j}^{c, 1} = 1$ , the likelihood contribution for subject $i$ at time $j$ , becomes degenerate after entry into the absorbing state. This implies that all contributions to the likelihood must have been preceded by a non absorbing state day [ $Y_{i (j - 1)} = Y_{i (j - 2)} = \dots = Y_{i 0} = 2$ ] and, implicitly, $μ_{i j}^{c} = P (Y_{i j} = 1 ∣ X_{i}, Z Z_{i}, Y_{i (j - 1)} = 2, Y_{i (j - 2)} = 2, \dots Y_{i 0} = 2)$ . Thus, if we define $W$ to be the entry time into the absorbing state, we can write $μ_{i j}^{c} = P (W = t_{i j} ∣ X_{i}, Z_{i}, W \geq t_{i j})$ which is a discrete hazard.

4. When severity is available: Cumulative probability models

We now incorporate severity into the analysis. That is, we expand state 2 (i.e. alive) into health states with higher values denoting better health. Our interest is still in the outcome $Y_{i j} = 1$ versus $Y_{i j} > 1$ , but we want to include all observed states in the model to potentially improve precision of the treatment effect estimates. We let $Y Y_{i}$ be a $m_{i} -$ vector of ordinal response values with $Y_{i j} \in {1, \dots, K}$ . In settings with an ordinal outcome, a popular model choice is the cumulative probability model (CPM) with the proportional odds (PO) assumption.^12,29 Let $π_{i j, k}^{po} = P (Y_{i j} \leq k ∣ X_{i}, Z Z_{i})$ , then the CPM with the PO assumption for subject $i$ at visit $j$ is given by

logit (π_{i j, k}^{po}) = \underset{α_{k}^{po} (t_{i j})}{\underset{⏟}{[α_{k}^{po} + α_{t}^{po} t_{i j}]}} + \underset{β_{x}^{po} (t_{i j})}{\underset{⏟}{[β_{x}^{po} + β_{x t}^{po} t_{i j}]}} X_{i} + Z Z_{i} β_{z}^{po},

(4)

where the

α_{k}^{po} (t_{i j})

function captures the

t_{i j} -

specific log-odds that

Y_{i j} \leq k

in the control arm (

X_{i} = 0

) when

Z Z_{i} = 0 0

, and the

β_{x}^{po} (t_{i j})

function captures the

t_{i j} -

specific difference in the log-odds between the treatment (

X_{i} = 1

) and the control (

X_{i} = 0

). For ease of exposition, we assume the log cumulative odds association with baseline covariates,

Z Z_{i}

, does not change with

t_{i j}

. The PO assumptions in model (4) imply that (i) the association between the ordinal outcome and

t_{i j}

does not depend on

k

and is captured with

α_{t}^{po}

, and (ii) cross-sectional differences in the log-odds between treatment and control do not depend on

k

. If one believes that the relationship between the log cumulative odds and

t_{i j}

X_{i}

varies according to outcome level

k

, the PO assumption can be relaxed and partial proportional odds (PPO) model assumptions can be made³⁰ for which

π_{i j, k}^{ppo} = P (Y_{i j} \leq k ∣ X_{i}, Z Z_{i})

can be expressed as

logit (π_{i j, k}^{ppo}) = \underset{α α_{k}^{ppo} (t_{i j})}{\underset{⏟}{[α_{k}^{ppo} + {(α α_{t}^{ppo})}^{T} τ τ t_{i j}]}} + \underset{β β_{x}^{ppo} (t_{i j})}{\underset{⏟}{[{(β β_{x}^{ppo})}^{T} τ τ + {(β β_{x t}^{ppo})}^{T} τ τ t_{i j}]}} X_{i} + Z Z_{i} β β_{z}^{ppo} .

(5)

Non-proportionality vector

τ τ

can be expressed, in a

K

-level outcome variable,

[1, I (k = 2), I (k = 3), \dots I (k = K - 1)]^{T}

, so that covariate-outcome associations are allowed to differ across dichotomizations of the log cumulative odds. Model (5) is considered a PPO model because the PO assumption only holds for a subset of the relationships with the

K - 1

log cumulative odds (e.g. for

Z Z_{i}

where

τ τ

is absent). For a four-level outcome variable, the association of

X_{i} t_{i j}

with the log cumulative odds is given by

{(β β_{x t}^{ppo})}^{T} τ τ = [β_{x t, 1}^{ppo}, β_{x t, 2}^{ppo}, β_{x t, 3}^{ppo}] τ τ = β_{x t, 1}^{ppo} + β_{x t, 2}^{ppo} I (k = 2) + β_{x t, 3}^{ppo} I (k = 3) .

To write this marginal CPM out explicitly,

\begin{aligned} logit (π_{i j, 1}^{ppo}) & = [α_{1}^{ppo} + α_{t, 1}^{ppo} t_{i j}] + [β_{x, 1}^{ppo} + β_{x t, 1}^{ppo} t_{i j}] X_{i} + Z Z_{i} β β_{z}^{ppo}, \end{aligned}

(6)

\begin{aligned} logit (π_{i j, 2}^{ppo}) & = [α_{2}^{ppo} + (α_{t, 1}^{ppo} + α_{t, 2}^{ppo}) t_{i j}] + [(β_{x, 1}^{ppo} + β_{x, 2}^{ppo}) + (β_{x t, 1}^{ppo} + β_{x t, 2}^{ppo}) t_{i j}] X_{i} + Z Z_{i} β β_{z}^{ppo}, \end{aligned}

(7)

\begin{aligned} logit (π_{i j, 3}^{ppo}) & = [α_{3}^{ppo} + (α_{t, 1}^{ppo} + α_{t, 3}^{ppo}) t_{i j}] + [(β_{x, 1}^{ppo} + β_{x, 3}^{ppo}) + (β_{x t, 1}^{ppo} + β_{x t, 3}^{ppo}) t_{i j}] X_{i} + Z Z_{i} β β_{z}^{ppo} . \end{aligned}

(8)

$β β_{x t}^{ppo} = [β_{x t, 1}^{ppo}, 0, 0]$ implies PO (i.e. forcing the $β_{x t, 2}^{ppo}$ and $β_{x t, 3}^{ppo}$ to be equal to $0$ induces a PO assumption). Thus, we may reduce the dimensionality of the model by constraining the treatment effect to be common across outcome levels. We note that, while the PO assumption forces there to be a common treatment effect across outcome levels, it does not restrict the possible range of the effect. For example, $β β_{x t}^{ppo} = [β_{x t, 1}^{ppo}, β_{x t, 2}^{ppo}, β_{x t, 2}^{ppo}]$ restricts the treatment effect to be common for $k = 2$ and $k = 3$ but permits a different treatment effect for $k = 1$ . While including these restrictions can improve precision, one must carefully consider whether the restrictions made are reasonable (i.e. whether they will induce appreciable bias). In the presence of one absorbing state, the relationship between $Y_{i j}$ and $t_{i j}$ likely induces a violation of the PO assumption, and with two absorbing states, the PO assumptions is necessarily violated. For example, in the right panel of Figure 1, the PO assumption for the $t_{i j}$ parameter would imply that the lines separating response states are parallel with one another. This assumption is clearly violated in the CLOVERS data. As a consequence, restricting $α α_{t}^{ppo}$ to satisfy the PO assumption will lead to a severely misspecified model.

As in the binary outcome setting, the CPMs in equations (4) and (5) only describe a portion of the multivariate distribution $[Y Y_{i} | X_{i}, Z Z_{i}]$ . To conduct likelihood-based analyses, we need to specify a conditional, response dependence structure that allows us to identify the remaining $[Y Y_{i} | X_{i}, Z Z_{i}]$ . As before, we appeal to the first order transition model which was introduced by Lee and Daniels^10,11 and discussed in Schildcrout:¹³

\log (\frac{μ_{i j, k}^{c}}{μ_{i j, K}^{c}}) = Δ_{i j, k} (X_{i}, Z Z_{i}) + \sum_{l = 1}^{K - 1} γ^{k l} (X_{i}, Z Z_{i}) I {Y_{i (j - 1)} = l},

(9)

where

μ_{i j, k}^{c} = P (Y_{i j} = k | Y_{i (j - 1)}, X_{i}, Z Z_{i})

γ^{k l}

is the log-relative risk ratio parameter for the transition from state

k

to state

l

that could depend on

X_{i}

and

Z Z_{i}

, and

Δ_{i j, k} (X_{i}, Z Z_{i})

is an implicit function of the data and the parameters that links the CPMs to the response dependence model in equation (9). Similar to the binary outcome case, fitting this model requires iterating between calculating

Δ_{i j, k} (X_{i}, Z Z_{i})

and updating model parameter estimates using Newton–Raphson algorithm. Details on the calculation of the score function and Fisher information matrix are described by Lee and Daniels.^10,11 As is the case with binary response data, if state

k = 1

is absorbing then

γ^{11} = \infty

, which induces constraints on the model and fitting algorithm.¹³

5. Simulations

5.1. Aims and data generating mechanism

Our simulation study aims to (i) characterize the validity, efficiency, type I error and power of the described models for estimating an end-of-study treatment effect when one outcome state is absorbing; (ii) highlight the potential benefit of using longitudinal data under various dropout mechanisms; (iii) study robustness of the estimation procedures to model misspecification; and (iv) compare the efficiency of marginalized transition models and Cox proportional hazard models in estimating an end-of-study treatment effect.

We begin by simulating data from a randomized clinical trial where $N = 1, 000$ participants are assigned to either treatment ( $X_{i} = 1$ ) or control ( $X_{i} = 0$ ) with equal probability. Let $Y_{i j} \in {1, 2, 3, 4}$ be a four-level outcome variable with the lowest state, $k = 1$ , being absorbing, $Z_{i}$ be a binary covariate with $P (Z_{i} = 1) = 0.5$ , and $t_{i j} = \log_{28} (j + 1)$ , $j = 0, \dots, 27$ . We generate data according to the CPM in equation (5) with $α α_{k}^{ppo} = (- 4.75, - 1.5, 2)$ , $α α_{t}^{ppo} = (2, - 1.75, - 4)$ , $β β_{x}^{ppo} = (0, 0, 0)$ , $β β_{x t}^{ppo} = (- 0.35, 0, 0)$ and $β β_{z}^{ppo} = (2, 0, 0)$ . We chose model coefficients so that approximately 15% of subjects enter the absorbing state $k = 1$ during follow-up, and the end-of-study treatment effect log-odds ratio, when $j = 27$ , is $β_{x, 1}^{ppo} + β_{x t, 1}^{ppo} \cdot \log_{28} (j + 1) = β_{x, 1}^{ppo} + β_{x t, 1}^{ppo} \cdot 1 = - 0.35$ . Note that with the chosen parameter values, the PO assumption for the treatment effect holds. Let $Γ Γ$ be the $3 \times 3$ dependence model parameter matrix with elements $γ^{k l}$ denoting the log-relative risk ratio for transitioning from state $l$ to state $k$ and with state 4 being the reference state. We set

Γ Γ = (\begin{matrix} 25 & 2 & 1 \\ - 25 & 4 & 2 \\ - 25 & 1 & 4 \end{matrix}) .

The choice of

(γ^{11}, γ^{21}, γ^{31}) = (25, - 25, - 25)

generates the

Y_{i j} = 1

absorbing state. To illustrate the potential benefit of using longitudinal modeling approaches over cross-sectional approaches, we study each estimation procedure under complete follow-up, and under two dropout at random (DAR) mechanisms of different severity: a mild DAR where each subject

i

has a probability of dropping out at time

j

equal to

0.01

if the outcome state at

t_{i (j - 1)}

was observed and

Y_{i (j - 1)} > 1

, and a severe DAR in which we increase the dropout probability to

0.04

. Approximately, 20% and 60% of subjects had missing information on the outcome state at the end of follow-up in the mild and severe DAR scenarios, respectively.

To evaluate and compare power across analysis procedures, we increase the absolute magnitude of the treatment effect from $-$ 0.35 to $-$ 0.50 by setting $β β_{x}^{ppo} = (0, 0, 0)$ and $β β_{x t}^{ppo} = (- 0.5, 0, 0)$ in the main simulation scenario described above. The coefficients were chosen so that a CPM with a cross-sectional ordinal outcome and a single treatment effect across all outcome dichotomizations would yield approximately $80 %$ power. Type I error rates were also evaluated by simulating data under $β β_{x}^{ppo} = β β_{x t}^{ppo} = (0, 0, 0)$ .

To characterize robustness of the estimation procedures to model misspecification, we consider the complete follow-up setting described above and then modify coefficient values to simulate two additional scenarios: one where the PO assumption for the treatment effect does not hold, and one where the PO assumption holds for the treatment effect but not for the covariate parameter. Misspecification of the PO assumption for the treatment effect was simulated by setting $β β_{x t}^{ppo} = (- 0.35, - 0.20, - 0.20)$ , while misspecification of the PO assumption for the covariate parameter was simulated by setting $β β_{z}^{ppo} = (2, 0, - 1)$ .

Finally, because the Cox proportional hazards model is a popular choice in settings like those considered here, we sought to compare marginalized transition models and Cox proportional hazards models. To do so we generated and fit data under 1) the marginalized transition model described at the beginning of the section and 2) a Cox proportional hazard model with a Weibull baseline hazard. Corresponding to 2), survival times were generated from $h (t ∣ X_{i}, Z_{i}) = h_{0} (t) exp (β_{x} + β_{z} Z_{i})$ with $β_{x} = - 0.35$ and $β_{z} = 2$ . The baseline hazard was specified as $h_{0} (t) = λ ν t^{ν - 1}$ , with $ν = 0.5$ and $λ = 0.01$ chosen so that the resulting true risk difference was close to that obtained under the marginalized transition model. The true risk difference parameter value was $- 0.044$ when generating data under MTM $_{PO}$ and $- 0.048$ under when generating data under a survival model with a Weibull hazard function. Since estimates from marginalized transition models are on an outcome-marginal, cumulative incidence scale, while the estimates from a Cox proportional hazard model are on a conditional (on not yet having experienced the event or still being the in risk set) scale, we compare parameter and uncertainty estimates for the two estimation strategies on a common, fully marginal (i.e. marginalized over $Z_{i}$ ) risk difference scale with standard errors calculated using the multivariate delta method. Details on the computation of the risk difference and standard errors are described in the Supplementary Material (Section S1).

5.2. Estimation target and methods

Across all simulations, we are interested in estimating an end of follow-up treatment effect on an absorbing state. For a comprehensive comparison, we estimate the treatment effect with the procedures described above and summarized in Table 2. Each procedure is identified by a distinct combination of data structure (cross-sectional or longitudinal), outcome type (binary or ordinal), model (GLM, CPM, MTM, or Cox proportional hazard models), and, when applicable, a baseline treatment effect constraint (i.e. is the treatment effect constrained to be zero at baseline?). Within the CPM and MTM models we further explore each procedure under different PO assumptions for the treatment effect. Based on these possible combinations, we obtain 13 different modeling approaches.

Table 2.
Estimation procedures compared in the simulation study.

Data structure Outcome type Model Constrained Formula Estimation target

Cross-sectional Binary GLM No $logit (P (Y_{i 27}^{} = 1 | X_{i}, Z_{i})) = β_{0} + β_{x} X_{i} + β_{z} Z_{i}$ $β_{x}$

Ordinal CPM No $logit (P (Y_{i 27} \leq k | X_{i}, Z_{i})) = α_{0, k} + β_{z} Z_{i} + \dots$

CPM $_{PO}$ $X_{i} β_{x}$ $β_{x}$

CPM $_{PPO 1}$ $X_{i} [β_{x, 1} + β_{x, 2} I (Y_{i 27} \in {2, 3})]$ $β_{x, 1}$

CPM $_{PPO 2}$ $X_{i} [β_{x, 1} + β_{x, 2} I (Y_{i 27} = 2) + β_{x, 3} I (Y_{i 27} = 3)]$ $β_{x, 1}$

Longitudinal Binary MTM No $logit (P (Y_{i j}^{} = 1 | X_{i}, Z_{i})) = β_{0} + β_{x} X_{i} + β_{t} t_{i j} + β_{x t} X_{i} t_{i j} + β_{z} Z_{i}$ $β_{x}$ + $β_{x t}$

MTM Yes $logit (P (Y_{i j}^{*} = 1 | X_{i}, Z_{i})) = β_{0} + β_{t} t_{i j} + β_{x t} X_{i} t_{i j} + β_{z} Z_{i}$ $β_{x t}$

Ordinal MTM $logit (P (Y_{i j} \leq k | X_{i}, Z_{i})) = α_{0, k} + t_{i j} [β_{t, 1} + β_{t, 2} I (Y_{i j} = 2) + β_{t, 3} I (Y_{i j} = 3)] + β_{z} Z_{i} + \dots$

MTM $_{PO}$ No $β_{x} X_{i} + X_{i} t_{i j} β_{x t}$ $β_{x}$ + $β_{x t}$

MTM $_{PO}$ Yes $X_{i} t_{i j} β_{x t}$ $β_{x t}$

MTM $_{PPO 1}$ No $β_{x} X_{i} + X_{i} t_{i j} [β_{x t, 1} + β_{x t, 2} I (Y_{i j} \in {2, 3})]$ $β_{x}$ + $β_{x t, 1}$

MTM $_{PPO 1}$ Yes $X_{i} t_{i j} [β_{x t, 1} + β_{x t, 2} I (Y_{i j} \in {2, 3})]$ $β_{x t, 1}$

MTM $_{PPO 2}$ No $β_{x} X_{i} + X_{i} t_{i j} [β_{x t, 1} + β_{x t, 2} I (Y_{i j} = 2) + β_{x t, 3} I (Y_{i j} = 3)]$ $β_{x}$ + $β_{x t, 1}$

MTM $_{PPO 2}$ Yes $X_{i} t_{i j} [β_{x t, 1} + β_{x t, 2} I (Y_{i j} = 2) + β_{x t, 3} I (Y_{i j} = 3)]$ $β_{x t, 1}$

Survival COXPH No $h (t_{i}; X_{i}, Z_{i}) = h_{0} (t_{i}) \exp (β_{x} X_{i} + β_{z} Z_{i})$ $β_{x}$

Data structure	Outcome type	Model	Constrained	Formula	Estimation target
Cross-sectional	Binary	GLM	No	$logit (P (Y_{i 27}^{*} = 1 \| X_{i}, Z_{i})) = β_{0} + β_{x} X_{i} + β_{z} Z_{i}$	$β_{x}$
	Ordinal	CPM	No	$logit (P (Y_{i 27} \leq k \| X_{i}, Z_{i})) = α_{0, k} + β_{z} Z_{i} + \dots$
		CPM $_{PO}$		$X_{i} β_{x}$	$β_{x}$
		CPM $_{PPO 1}$		$X_{i} [β_{x, 1} + β_{x, 2} I (Y_{i 27} \in {2, 3})]$	$β_{x, 1}$
		CPM $_{PPO 2}$		$X_{i} [β_{x, 1} + β_{x, 2} I (Y_{i 27} = 2) + β_{x, 3} I (Y_{i 27} = 3)]$	$β_{x, 1}$
Longitudinal	Binary	MTM	No	$logit (P (Y_{i j}^{*} = 1 \| X_{i}, Z_{i})) = β_{0} + β_{x} X_{i} + β_{t} t_{i j} + β_{x t} X_{i} t_{i j} + β_{z} Z_{i}$	$β_{x}$ + $β_{x t}$
		MTM	Yes	$logit (P (Y_{i j}^{*} = 1 \| X_{i}, Z_{i})) = β_{0} + β_{t} t_{i j} + β_{x t} X_{i} t_{i j} + β_{z} Z_{i}$	$β_{x t}$
	Ordinal	MTM		$logit (P (Y_{i j} \leq k \| X_{i}, Z_{i})) = α_{0, k} + t_{i j} [β_{t, 1} + β_{t, 2} I (Y_{i j} = 2) + β_{t, 3} I (Y_{i j} = 3)] + β_{z} Z_{i} + \dots$
		MTM $_{PO}$	No	$β_{x} X_{i} + X_{i} t_{i j} β_{x t}$	$β_{x}$ + $β_{x t}$
		MTM $_{PO}$	Yes	$X_{i} t_{i j} β_{x t}$	$β_{x t}$
		MTM $_{PPO 1}$	No	$β_{x} X_{i} + X_{i} t_{i j} [β_{x t, 1} + β_{x t, 2} I (Y_{i j} \in {2, 3})]$	$β_{x}$ + $β_{x t, 1}$
		MTM $_{PPO 1}$	Yes	$X_{i} t_{i j} [β_{x t, 1} + β_{x t, 2} I (Y_{i j} \in {2, 3})]$	$β_{x t, 1}$
		MTM $_{PPO 2}$	No	$β_{x} X_{i} + X_{i} t_{i j} [β_{x t, 1} + β_{x t, 2} I (Y_{i j} = 2) + β_{x t, 3} I (Y_{i j} = 3)]$	$β_{x}$ + $β_{x t, 1}$
		MTM $_{PPO 2}$	Yes	$X_{i} t_{i j} [β_{x t, 1} + β_{x t, 2} I (Y_{i j} = 2) + β_{x t, 3} I (Y_{i j} = 3)]$	$β_{x t, 1}$
	Survival	COXPH	No	$h (t_{i}; X_{i}, Z_{i}) = h_{0} (t_{i}) \exp (β_{x} X_{i} + β_{z} Z_{i})$	$β_{x}$

Note: For models with longitudinal outcomes only the marginal (population averaged) model is included in the table as the conditional, response dependence model is the same across procedures. To ease readability, the subscript $p p o$ was omitted from the models with an ordinal outcome. $Y_{i j}^{*}$ is a binary outcome (dead vs alive) for subject $i$ at time $j$ . $Y_{i j}$ is an ordinal outcome for subject $i$ at time $j$ . For the longitudinal ordinal outcome scenarios, the first line indicates the part of the model that is common across models. Constrained (yes/no) indicates whether the treatment effect at baseline was constrained to be equal to zero. $h (\cdot)$ is the hazard function with $h_{0} (\cdot)$ being the baseline hazard.

Abbreviations: GLM = generalized linear model. MTM = marginalized transition model. CPM = cumulative probability model. PO = proportional odds. PPO = partial proportional odds.

First, we consider the binary outcome $Y_{i j} = 1$ versus $Y_{i j} > 1$ (i.e. $Y_{i j} \in {2, 3, 4}$ ) and fit a cross-sectional GLM for binary data at the last day of follow-up ( $j = 27$ or equivalently $t_{i 27} = 1$ ), along with two longitudinal MTMs for binary outcomes that either do or do not constrain the treatment effect to be zero at baseline. Next, we consider the four outcome states separately, and we evaluate three cross-sectional CPM specifications at the end of follow-up: one assuming a common treatment effect across all the dichotomization of $Y_{i j}$ (i.e. PO assumption on the treatment effect, CPM $_{PO}$ ), one allowing the treatment effect for the $k = 1$ dichotomization of the log cumulative odds to be different from the common treatment effect for $k = 2$ and $k = 3$ (CPM $_{PPO 1}$ ), and one allowing the treatment effect to be different for each dichomtoization of $Y_{i j}$ (CPM $_{PPO 2}$ ). Notice that CPM $_{PPO 1}$ and CPM $_{PPO 2}$ relax the PO restriction imposed in CPM $_{PO}$ . We then consider the corresponding longitudinal version of these ordinal models—MTM $_{PO}$ , MTM $_{PPO 1}$ , and MTM $_{PPO 2}$ —each evaluated with and without the additional constraint that the treatment effect at baseline is zero. Finally, for completeness, we include a Cox proportional hazard model in which the event of interest is defined as the entry into the absorbing state.

5.3. Performance measures

For each of the 5,000 replications, we estimated regression coefficients for the treatment effect, their standard error and the $95 %$ confidence interval. We then averaged across all replications to obtain the following performance measures used to assess the validity and efficiency of the procedures described in the previous section:

percent bias in the treatment effect estimate

empirical standard error (SE) of the treatment effect estimate

average standard error estimate (SEE)

root mean square error (RMSE)

coverage probability (CP) of the $95 %$ confidence interval

relative efficiency of each estimation procedure to the cross-sectional generalized linear model for binary outcome (GLM) (i.e. the empirical variance of GLM divided by empirical variances of all other estimation procedures).

5.4. Results

Results on validity of the marginalized transition models are presented in Table 3. Under complete follow-up and mild DAR, all estimation procedures yielded approximately unbiased treatment effect estimates and standard errors, and confidence intervals had appropriate coverage. Under severe DAR, estimates of the association between the treatment and the dichotomized response value $Y_{i 27} = 1$ versus $Y_{i 27} > 1$ tended to be biased (up to 8%) with cross-sectional modeling procedures. This result is observed because even though the dropout mechanism is at random for longitudinal analyses, it is not at random for the cross-sectional analyses. Interestingly, there was effectively no bias for the cross-sectional ordinal regression model with the treatment effect assumed to be common across all dichotomizations (CPM $_{PO}$ ).

Table 3.
Percent bias, empirical standard errors, average estimated standard error, and 95% confidence interval width for the log odds ratio between treatments $X_{i} = 1$ versus $X_{i} = 0$ on dichotomized response $Y_{i 27} = 1$ versus $Y_{i 27} > 1$ across 5,000 replicates.

Complete follow-up DAR (mild) DAR (severe)

Data structure Outcome type Estimation procedure Constrained $%$ Bias SE (SEE) RMSE CP $%$ Bias SE (SEE) RMSE CP $%$ Bias SE (SEE) RMSE CP

Cross-sectional Binary GLM No 0.893 0.182 (0.180) 0.182 0.947 −1.860 0.197 (0.192) 0.197 0.947 −8.207 0.240 (0.237) 0.241 0.947

Longitudinal Binary MTM No 0.935 0.182 (0.180) 0.182 0.949 0.675 0.192 (0.188) 0.192 0.949 0.631 0.219 (0.215) 0.219 0.947

Yes 0.833 0.166 (0.163) 0.166 0.950 0.548 0.179 (0.174) 0.179 0.948 0.546 0.212 (0.207) 0.212 0.944

Cross-sectional Ordinal CPM $_{PO}$ No 0.807 0.120 (0.118) 0.120 0.946 0.700 0.135 (0.133) 0.135 0.948 −0.671 0.192 (0.188) 0.192 0.944

CPM $_{PPO 1}$ No 0.926 0.175 (0.174) 0.175 0.949 −1.729 0.189 (0.185) 0.189 0.946 −7.706 0.232 (0.228) 0.233 0.944

CPM $_{PPO 2}$ No 0.682 0.180 (0.178) 0.180 0.951 −2.088 0.194 (0.190) 0.194 0.948 −8.377 0.238 (0.234) 0.240 0.944

Longitudinal Ordinal MTM $_{PO}$ No 0.460 0.071 (0.071) 0.071 0.948 0.510 0.077 (0.076) 0.077 0.947 0.342 0.095 (0.094) 0.095 0.949

Yes 0.441 0.070 (0.070) 0.070 0.946 0.499 0.075 (0.075) 0.075 0.948 0.318 0.093 (0.093) 0.093 0.950

MTM $_{PPO 1}$ No 0.543 0.154 (0.152) 0.154 0.944 0.407 0.166 (0.162) 0.166 0.945 0.396 0.199 (0.194) 0.199 0.943

Yes 0.567 0.153 (0.151) 0.153 0.944 0.418 0.165 (0.161) 0.165 0.944 0.394 0.199 (0.194) 0.199 0.941

MTM $_{PPO 2}$ No 0.377 0.162 (0.159) 0.162 0.946 0.154 0.173 (0.169) 0.173 0.945 0.200 0.206 (0.201) 0.206 0.943

Yes 0.418 0.160 (0.157) 0.160 0.949 0.180 0.172 (0.167) 0.172 0.946 0.209 0.205 (0.200) 0.205 0.944

				Complete follow-up	DAR (mild)	DAR (severe)
Cross-sectional	Binary	GLM	No	0.893	0.182 (0.180)	0.182	0.947	−1.860	0.197 (0.192)	0.197	0.947	−8.207	0.240 (0.237)	0.241	0.947
Longitudinal	Binary	MTM	No	0.935	0.182 (0.180)	0.182	0.949	0.675	0.192 (0.188)	0.192	0.949	0.631	0.219 (0.215)	0.219	0.947
			Yes	0.833	0.166 (0.163)	0.166	0.950	0.548	0.179 (0.174)	0.179	0.948	0.546	0.212 (0.207)	0.212	0.944
Cross-sectional	Ordinal	CPM $_{PO}$	No	0.807	0.120 (0.118)	0.120	0.946	0.700	0.135 (0.133)	0.135	0.948	−0.671	0.192 (0.188)	0.192	0.944
		CPM $_{PPO 1}$	No	0.926	0.175 (0.174)	0.175	0.949	−1.729	0.189 (0.185)	0.189	0.946	−7.706	0.232 (0.228)	0.233	0.944
		CPM $_{PPO 2}$	No	0.682	0.180 (0.178)	0.180	0.951	−2.088	0.194 (0.190)	0.194	0.948	−8.377	0.238 (0.234)	0.240	0.944
Longitudinal	Ordinal	MTM $_{PO}$	No	0.460	0.071 (0.071)	0.071	0.948	0.510	0.077 (0.076)	0.077	0.947	0.342	0.095 (0.094)	0.095	0.949
			Yes	0.441	0.070 (0.070)	0.070	0.946	0.499	0.075 (0.075)	0.075	0.948	0.318	0.093 (0.093)	0.093	0.950
		MTM $_{PPO 1}$	No	0.543	0.154 (0.152)	0.154	0.944	0.407	0.166 (0.162)	0.166	0.945	0.396	0.199 (0.194)	0.199	0.943
			Yes	0.567	0.153 (0.151)	0.153	0.944	0.418	0.165 (0.161)	0.165	0.944	0.394	0.199 (0.194)	0.199	0.941
		MTM $_{PPO 2}$	No	0.377	0.162 (0.159)	0.162	0.946	0.154	0.173 (0.169)	0.173	0.945	0.200	0.206 (0.201)	0.206	0.943
			Yes	0.418	0.160 (0.157)	0.160	0.949	0.180	0.172 (0.167)	0.172	0.946	0.209	0.205 (0.200)	0.205	0.944

Note: Constrained (yes/no) indicates whether the treatment effect at baseline was constrained to be equal to zero.

Abbreviations: GLM = generalized linear model; MTM = marginalized transition model; CPM = cumulative probability model; PO = proportional odds; PPO = partial proportional odds; SE = empirical standard error; SEE = average estimated standard error; RMSE = root mean squared error; CP = coverage probability of the 95% confidence interval; DAR = dropout at random.

Figure 2 shows the treatment effect relative efficiency of each estimation procedure to the cross-sectional generalized linear model with a binary outcome (GLM). Under complete follow-up, considering a longitudinal binary outcome rather than a cross-sectional binary outcome did not result in any precision gain. However, to the extent that dropout occurs, gains were realized (e.g. $6 %$ under mild DAR and $20 %$ under severe DAR). Constraining the treatment effect to be zero further improved efficiency by 15% to 20% under the complete follow-up and mild DAR scenarios, but there was less of an improvement under severe DAR. We observed that the PO assumption for scalar or longitudinal models yielded enormous efficiency gains. Even though we include these results, it should be noted that the proportional odds assumption is unlikely to be acceptable (even in the absence of evidence that it is wrong) if interest is specifically in the treatment effect on the binary $Y_{i 27} = 1 versus Y_{i 27} > 1$ outcome. Pertaining to the scalar ordinal models, assuming that the treatment effect for the $k = 1$ dichotomization of the log cumulative odds is distinct from the common treatment effect for $k = 2$ and $k = 3$ (MTM $_{PPO 1}$ ), provided an 8% to 9% efficiency gain over cross-sectional generalized linear model with a binary outcome, while assuming that the treatment effect differs for each dichotomization of $Y_{i j}$ (MTM $_{PPO 2}$ ) provided no efficiency improvements. Finally, models that exploit the longitudinal and ordinal outcome data and that impose modest restrictions on the treatment effect (MTM $_{PPO 1}$ and MTM $_{PPO 2}$ ) yielded 25% to 45% percent efficiency gains over the cross-sectional generalized linear model with a binary outcome. In contrast to what we observed in the binary outcome setting, constraining the treatment effect to be zero at baseline had little to no impact on efficiency in the longitudinal ordinal models.

Figure 2.

Relative efficiency of different estimation procedures to a cross-sectional binary logistic regression (GLM) in estimating the association between the treatment $X_{i}$ and the dichotomized response value $Y_{i 27} = 1$ versus $Y_{i 27} > 1$ . DAR, dropout at random. MTM, marginalized transition model. CPM, cumulative probability model. PO, proportional odds. PPO, partial proportional odds.

Consistent with the observed validity when the above models are properly specified, we observed that type 1 error rates maintained their nominal 0.05 levels (see Supplementary Table S3), and consistent with the efficiency gains observed above, we observed power improvements associated with the use of longitudinal and ordinal response models and with the placement of constraints on the treatment effect overall and at baseline (for the longitudinal binary response setting; see Table 4).

Table 4.

Percent bias, empirical standard errors, average estimated standard error, coverage probability of the 95% confidence interval and power for the odds ratio between treatments $X_{i}$ = 1 versus $X_{i}$ = 0 on dichotomized response $Y_{i 27}$ = 1 versus $Y_{i 27} > 1$ across 5,000 replicates.

				Complete follow-up					DAR (mild)					DAR (severe)
Data Structure	Outcome Type	Estimation Procedure	Con-strained	$%$ Bias	SE (SEE)	RMSE	CP	Power	$%$ Bias	SE (SEE)	RMSE	CP	Power	$%$ Bias	SE (SEE)	RMSE	CP	Power
Cross-sectional	Binary	GLM	No	0.203	0.189 (0.185)	0.189	0.947	0.772	−2.157	0.200 (0.197)	0.200	0.951	0.705	−9.715	0.244 (0.241)	0.249	0.943	0.469
Longitudinal	Binary	MTM	No	0.189	0.188 (0.184)	0.188	0.948	0.774	0.272	0.196 (0.192)	0.196	0.950	0.737	0.002	0.223 (0.219)	0.223	0.947	0.629
			Yes	−0.001	0.170 (0.165)	0.170	0.948	0.847	0.082	0.180 (0.176)	0.180	0.948	0.804	−0.217	0.213 (0.209)	0.213	0.947	0.668
Cross-sectional	Ordinal	CPM $_{PO}$	No	0.449	0.120 (0.119)	0.120	0.951	0.991	0.521	0.135 (0.134)	0.135	0.948	0.962	−1.707	0.191 (0.189)	0.192	0.947	0.738
		CPM $_{PPO 1}$	No	0.429	0.182 (0.178)	0.182	0.947	0.802	−1.967	0.192 (0.189)	0.193	0.951	0.736	−8.975	0.236 (0.232)	0.240	0.939	0.505
		CPM $_{PPO 2}$	No	0.101	0.187 (0.183)	0.187	0.948	0.780	−2.274	0.198 (0.194)	0.198	0.951	0.714	−9.746	0.242 (0.238)	0.247	0.941	0.477
Longitudinal	Ordinal	MTM $_{PO}$	No	0.014	0.071 (0.071)	0.071	0.950	1.000	0.027	0.076 (0.076)	0.076	0.951	1.000	−0.337	0.094 (0.094)	0.094	0.953	1.000
			Yes	−0.041	0.070 (0.070)	0.070	0.949	1.000	−0.033	0.075 (0.075)	0.075	0.952	1.000	−0.407	0.092 (0.093)	0.092	0.950	1.000
		MTM $_{PPO 1}$	No	0.064	0.158 (0.154)	0.158	0.943	0.900	0.113	0.168 (0.164)	0.168	0.946	0.857	−0.224	0.200 (0.197)	0.200	0.948	0.717
			Yes	0.151	0.157 (0.153)	0.157	0.943	0.900	0.189	0.167 (0.163)	0.167	0.946	0.860	−0.179	0.200 (0.196)	0.200	0.946	0.719
		MTM $_{PPO 2}$	No	−0.265	0.166 (0.162)	0.166	0.945	0.864	−0.084	0.176 (0.172)	0.176	0.946	0.823	−0.298	0.207 (0.204)	0.207	0.949	0.691
			Yes	−0.123	0.165 (0.160)	0.165	0.945	0.871	0.036	0.175 (0.170)	0.174	0.944	0.832	−0.219	0.206 (0.203)	0.206	0.947	0.697

Note: Constrained (yes/no) indicates whether the treatment effect at baseline was constrained to be equal to zero.

Table 5 summarizes the results when the model is misspecified. When the PO assumption for the treatment effect is misspecified, models that assume a common odds ratio across the three log cumulative odds dichotomizations (CPM $_{PO}$ and MTM $_{PO}$ ) yield substantially biased parameter estimates and inferences are invalid. That is, the PO assumption that yielded enormous efficiency gains can yield large biases if it is not correct. All other estimation procedures yield approximately valid inferences for the treatment effect. When the PO assumption for the covariate effect was violated and not acknowledged, as in the cases of all models tested in the simulation study, most estimation procedure performed well with biases in parameter estimates that tended to be less than 3% and with inferences that were approximately valid. In contrast, marginalized transition models for longitudinal ordinal outcome that assumed a different treatment effect for each dichotomization of $Y_{i j}$ (MTM $_{PPO 2}$ ) exhibited approximately a 5% bias.

Table 5.

Percent bias, empirical standard errors, average estimated standard error, and 95% confidence interval width for the log odds ratio between treatments $X_{i} = 1$ versus $X_{i} = 0$ on dichotomized response $Y_{i 27} = 1$ versus $Y_{i 27} > 1$ across 5,000 replicates in settings where at least some of the models are misspecified.

				Misspecification of $β β_{x t}$				Misspecification of $β β_{z}$
Data Structure	Outcome Type	Estimation Procedure	Constrained	$%$ Bias	SE (SEE)	RMSE	CP	$%$ Bias	SE (SEE)	RMSE	CP
Cross-sectional	Binary	GLM	No	1.16	0.178 (0.180)	0.178	0.957	1.080	0.182 (0.180)	0.182	0.949
Longitudinal	Binary	MTM	No	1.17	0.178 (0.179)	0.178	0.957	1.109	0.181 (0.179)	0.181	0.950
			Yes	1.03	0.163 (0.163)	0.163	0.949	1.134	0.165 (0.163)	0.165	0.949
Cross-sectional	Ordinal	CPM $_{PO}$	No	43.98	0.121 (0.119)	0.196	0.749	−0.561	0.119 (0.119)	0.119	0.950
		CPM $_{PPO 1}$	No	1.12	0.171 (0.173)	0.171	0.955	0.704	0.174 (0.172)	0.174	0.949
		CPM $_{PPO 2}$	No	1.12	0.177 (0.178)	0.177	0.955	−4.530	0.174 (0.174)	0.175	0.953
Longitudinal	Ordinal	MTM $_{PO}$	No	57.80	0.071 (0.071)	0.215	0.186	−0.872	0.071 (0.071)	0.071	0.951
			Yes	58.83	0.071 (0.070)	0.218	0.158	−0.691	0.070 (0.071)	0.070	0.955
		MTM $_{PPO 1}$	No	1.19	0.149 (0.150)	0.149	0.951	3.146	0.152 (0.148)	0.153	0.946
			Yes	1.11	0.148 (0.149)	0.148	0.951	2.696	0.151 (0.147)	0.151	0.946
		MTM $_{PPO 1}$	No	1.15	0.158 (0.158)	0.158	0.951	−5.202	0.155 (0.155)	0.156	0.951
			Yes	1.03	0.156 (0.157)	0.156	0.951	−5.253	0.153 (0.154)	0.154	0.950

Note: Constrained (yes/no) indicates whether the treatment effect at baseline was constrained to be equal to zero.

Table 6 summarizes results of the comparison between the marginalized transition models and the Cox proportional hazard model. When generating the data under a marginalized transition model with an ordinal outcome and a treatment effect assumed to be the same across all dichotomization of $Y_{i j}$ (MTM $_{PO}$ ), the estimated risk difference under a proportional hazards model exhibited a 5% bias under the complete data scenario, and the bias grew under DAR. This is expected since the Cox proportional hazard model is mildly misspecified under the data generating mechanism. Across the scenarios, treatment effect uncertainty for the Cox model tracked with the MTMs for binary data regardless of whether we imposed the constraint on the magnitude of the treatment effect at baseline. We find this result sensible since the three estimation strategies are based on models for longitudinal data with a binary, absorbing outcome. Analogous results were observed when we generated data from the Cox proportional hazard model and fit the data with the marginalized transition models for binary outcomes (i.e. 5−10% bias and similar efficiency under complete follow-up and the mild DAR scenario). The cross-sectional generalized linear model with a binary outcome was biased and inefficient under both DAR mechanisms.

Table 6.

Percent bias, empirical standard errors, average estimated standard error, and 95% confidence interval width for the fully marginal (with respect to $Z_{i}$ ) risk difference between treatments $X_{i} = 1$ versus $X_{i} = 0$ on dichotomized response $Y_{i 27} = 1$ versus $Y_{i 27} > 1$ across 5,000 replicates.

				Complete follow-up				DAR (mild)				DAR (severe)
Data structure	Outcome type	Estimation procedure	Constrained	$%$ Bias	SE (SEE)	RSME	CP	$%$ Bias	SE (SEE)	RSME	CP	$%$ Bias	SE (SEE)	RSME	CP
Data generating process: longitudinal ordinal outcome, MTM $_{PO}$ , no constraint
Cross-sectional	Binary	GLM	No	0.148	0.023 (0.022)	0.023	0.947	6.864	0.027 (0.026)	0.027	0.944	23.749	0.040 (0.040)	0.042	0.940
Longitudinal	Binary	MTM	No	0.209	0.022 (0.022)	0.022	0.947	−0.062	0.024 (0.023)	0.024	0.945	−0.095	0.027 (0.027)	0.027	0.948
			Yes	0.176	0.021 (0.020)	0.021	0.946	−0.136	0.022 (0.022)	0.022	0.946	−0.167	0.026 (0.026)	0.026	0.944
Cross-sectional	Ordinal	CPM $_{PO}$	No	0.421	0.015 (0.015)	0.015	0.947	9.710	0.019 (0.018)	0.019	0.941	31.655	0.032 (0.031)	0.035	0.917
		CPM $_{PPO 1}$	No	0.289	0.022 (0.022)	0.022	0.950	6.849	0.026 (0.025)	0.026	0.944	22.188	0.038 (0.038)	0.040	0.937
		CPM $_{PPO 2}$	No	0.019	0.022 (0.022)	0.022	0.948	6.428	0.026 (0.026)	0.026	0.943	21.243	0.039 (0.039)	0.041	0.938
Longitudinal	Ordinal	MTM $_{PO}$	No	0.245	0.009 (0.009)	0.009	0.949	0.292	0.010 (0.010)	0.010	0.950	0.134	0.012 (0.012)	0.012	0.951
			Yes	0.226	0.009 (0.009)	0.009	0.947	0.280	0.010 (0.010)	0.010	0.949	0.102	0.012 (0.012)	0.012	0.951
		MTM $_{PPO 1}$	No	0.057	0.019 (0.019)	0.019	0.943	−0.093	0.021 (0.020)	0.021	0.946	−0.137	0.025 (0.024)	0.025	0.941
			Yes	0.085	0.019 (0.019)	0.019	0.945	−0.080	0.021 (0.020)	0.021	0.946	−0.138	0.025 (0.024)	0.025	0.939
		MTM $_{PPO 2}$	No	−0.147	0.020 (0.020)	0.020	0.947	−0.382	0.022 (0.021)	0.022	0.944	−0.358	0.026 (0.025)	0.026	0.943
			Yes	−0.099	0.020 (0.020)	0.020	0.946	−0.349	0.021 (0.021)	0.021	0.944	−0.344	0.026 (0.025)	0.026	0.942
	Survival	COXPH		−4.997	0.022 (0.022)	0.022	0.951	−8.528	0.023 (0.023)	0.023	0.947	−17.964	0.025 (0.026)	0.027	0.941
Data generating process: survival model with Weibull baseline hazard
Cross-sectional	Binary	GLM	No	0.744	0.022 (0.022)	0.022	0.950	10.810	0.026 (0.026)	0.027	0.944	39.642	0.040 (0.039)	0.044	0.915
Longitudinal	Binary	MTM	No	4.600	0.023 (0.022)	0.023	0.942	5.936	0.024 (0.024)	0.025	0.938	11.079	0.029 (0.028)	0.029	0.944
			Yes	−7.619	0.021 (0.021)	0.021	0.949	−5.009	0.023 (0.023)	0.023	0.947	4.485	0.028 (0.028)	0.028	0.953
	Survival	COXPH		−1.097	0.021 (0.022)	0.021	0.951	−1.159	0.023 (0.023)	0.023	0.947	−1.026	0.025 (0.025)	0.025	0.953

Note: Constrained (yes/no) indicates whether the treatment effect at baseline was constrained to be equal to zero.

One interesting observation occurred under mild and severe DAR with the cross-sectional modeling procedures (CPM $_{PO}$ , CPM $_{PPO 1}$ and CPM $_{PPO 2}$ ). On the marginalized risk difference scale, the percent biases were substantially larger than they were on the log-odds ratio scale. For example, looking at the ordinal model where the treatment effect is assumed to be the same across all dichotomizations (CPM $_{PO}$ ), the treatment effect estimate was approximately unbiased on the log-odds ratio scale (Table 3) but exhibited a 30% bias on the risk difference scale. Even though the log-odds ratio for the treatment effect was not severely biased under DAR, the DAR mechanism induced biases on the intercept ( $α α_{t}$ ). When we combine small biases in the log-odds ratio with biases in prevalence, we observe substantial percent bias in the estimates of the risk differences.

We note that when severity is available, the treatment effect on mortality can be estimated using a discrete time multi-state survival model.³¹ While we did not include the multi-state survival model in the simulations due to its heavy computational burden, we fit it to estimate the treatment effect on mortality in the CLOVERS data in the next section. A detailed summary of runtime for the currently implemented MTM models is provided in Table A4. Under complete follow-up, average runtime per replicate were 1−2 minutes for longitudinal binary outcomes and 6−8 minutes for longitudinal ordinal outcomes. These runtimes are expected to improve with updates to the software. Computation time was negligible for the GLM, CPM, and Cox proportional hazards models (<0.01 minutes on average).

6. CLOVERS study: Data analysis

For the CLOVERS study analyses, we implemented the estimation procedures described in Section 5 to calculate the effect of a restrictive fluid strategy on mortality at the end of follow-up. The covariate vector $Z Z_{i}$ includes baseline information on patients’ age, sex, and SOFA score. For longitudinal models we included time using a restricted cubic spline with three degrees of freedom on the logarithm of day of follow-up, and added a time by treatment interaction. In addition to log-odds ratios, we estimated the fully marginal (with respect of the covariates) risk differences and compare the results with those from (i) a Cox proportional hazard model that includes the main effect of treatment, age, sex, and SOFA score as covariates, and (ii) a discrete time multi-state survival model with the same set of baseline covariates. In the multi-state survival model, we assumed that baseline covariates could have a different impact on each state transition; in other words, the multi-state survival model estimates a potentially different coefficient for each baseline covariate and each state transition. We estimated standard errors and 95 $%$ confidence intervals for the multi-state survival model via bootstrap sampling with 5,000 replicates. Figure 3 shows the covariate-adjusted log-odds ratio of a restrictive versus liberal fluid resuscitation strategy on mortality at the end of follow-up and the fully marginal (with respect to baseline covariates) risk difference. Overall, there is little to no evidence to suggest that implementing a restrictive fluid strategy impacts mortality. Estimated log-odds ratio and risk differences were largely consistent across all estimation procedures.

Figure 3.

Estimated covariate-adjusted log-odds ratios and fully marginal (with respect to baseline covariate) risk differences (RD) describing the effect of restrictive fluid strategy on mortality at the end of follow-up with $95 %$ confidence intervals (CIs). CS-B, cross-sectional binary outcome. L-B, longitudinal binary outcome. CS-O, cross-sectional ordinal outcome. L-O, longitudinal ordinal outcome. Surv, survival outcome. MTM, marginalized transition model. CPM, cumulative probability model. PO, proportional odds. PPO, partial proportional odds. MSM, multi-state survival model. COXPH, Cox proportional hazard model. SE, standard error.

Similar to the precision gains observed in the simulation studies for the log-odds ratio estimates, model uncertainty decreased when we used an ordinal outcome rather than a binary outcome: compared to the cross-sectional generalized linear model with a binary outcome, we observed a $3 %$ to $20 %$ reduction in the estimated standard error when using cross-sectional, ordinal outcome models. Not surprisingly, we observed the greatest reduction under the PO restriction for all dichotomizations of the outcome (estimated standard errors 0.164 vs 0.137). Focusing on the ordinal outcome models, estimated standard errors were lower with longitudinal models compared to cross-sectional models. For instance, under the PO restriction for all dichotomization of the outcome, the longitudinal model reduced the estimated standard error by $17 %$ (e.g. 0.137 vs 0.116). Consistent with precision comparisons in the simulations, estimated uncertainty was similar for cross-sectional binary and unconstrained longitudinal binary outcome models (standard errors 0.164 vs 0.163). When we constrained the treatment effect to be 0 at baseline, $β β_{x} = 0 0$ , and fit the marginalized transition model, the estimated standard error was reduced to 0.154. Also, consistent with precision comparisons in the simulations, the constrained marginalized transition models for longitudinal ordinal outcomes did not reduce estimated uncertainty when compared to their unconstrained counterparts.

When estimating risk differences, the estimated standard error for the treatment effect from the Cox proportional hazard model tracked with those of the marginalized regression models for longitudinal binary outcome. Similarly, the estimated standard error for the treatment effect from the multi-state survival model tracked with the less constrained MTMs for ordinal outcomes (e.g. estimated standard error was 0.015 in both the multi-state model and the longitudinal ordinal model where the treatment effect was assumed to differ across all outcome dichotomizations).

Likelihood ratio tests can be used to compare many of the model fits. Here, we focus on understanding whether relaxing the PO assumptions and/or the assumption of a null treatment effect at baseline improved model fit. Table 7 shows the maximized log-likelihood for the MTMs with longitudinal binary outcome, the CPMs with cross-sectional ordinal outcome, and the MTMs with longitudinal ordinal outcome. For the longitudinal binary outcome models, we observed insufficient evidence to suggest that relaxing the assumption of a null treatment effect at baseline improved model fit ( $P = 0.138$ ). This is expected since participants were randomized at baseline. Interestingly, under the longitudinal ordinal outcome models, relaxing the no baseline treatment effect assumption seemed to improve model fit ( $P < 0.001$ across all comparisons between the constrained and the unconstrained MTM $_{PO}$ , MTM $_{PPO 1}$ , and MTM $_{PPO 2}$ ). This could imply some degree of model misspecification of the early treatment effect. Also interesting, relaxing the PO assumption for the treatment effect in the cross-sectional and longitudinal ordinal outcome settings did not improve model fit. That is, the treatment effect is basically null, which is consistent with the proportional odds assumption (i.e. the treatment effect is zero for each outcome dichotomization).

Table 7.

Likelihood ratio tests for the CLOVERS study data analysis.

Data Structure	Outcome Type	Estimation Procedure	Constrained	Number of Parameters	log(L)	p-Value
Longitudinal	Binary	MTM	No	12	−1120.492	0.138 (vs constrained MTM)
			Yes	11	−1121.588	−
Cross-sectional	Ordinal	CPM $_{PPO 2}$	No	9	−892.733	0.619 (vs CPM $_{PO}$ )
		CPM $_{PPO 1}$	No	8	−893.166	0.762 (vs CPM $_{PO}$ )
		CPM $_{PO}$	No	7	−893.212	−
Longitudinal	Ordinal	MTM $_{PPO 2}$	No	34	−8077.089	<0.001 (vs constrained MTM $_{PPO 2}$ )
			Yes	33	−8085.270	0.605 (vs constrained MTM $_{PPO 1}$ )
		MTM $_{PPO 1}$	No	31	−8077.941	<0.001 (vs constrained MTM $_{PPO 1}$ )
			Yes	30	−8086.193	0.173 (vs constrained MTM $_{PO}$ )
		MTM $_{PO}$	No	28	−8080.103	<0.001 (vs constrained MTM $_{PO}$ )
			Yes	27	−8088.684	−

Note: Constrained (yes/no) indicates whether the treatment effect at baseline was constrained to be equal to zero.

Abbreviations: GLM = generalized linear model; MTM = marginalized transition model; CPM = cumulative probability model; PO = proportional odds; PPO = partial proportional odds; SE = empirical standard error; SEE = average estimated standard error; CP = coverage probability of the 95% confidence interval.

7. Discussion

We discussed and compared approaches to summarize the outcome-marginal baseline covariate-adjusted treatment effect on mortality at the end of a follow-up period in a randomized clinical trial when health state data are collected daily (i.e. longitudinally). We focused on how and when specific modeling choices, such as the use of longitudinal versus cross-sectional outcomes or of ordinal versus binary outcomes, can improve estimation efficiency compared to the commonly used cross-sectional binary logistic regression and survival analysis.

Together with the choice of outcome, we studied how estimation efficiency can be improved by imposing constraints on the form of regression model. In particular, we considered scenarios where the treatment effect at baseline was forced to be zero and, for the ordinal outcome, we examined scenarios where we imposed and/or relaxed PO assumptions. Constraining the treatment effect at baseline to zero improved estimation efficiency with longitudinal binary outcome models but not with longitudinal ordinal outcome models. This result could be explained by the number of parameters estimated in each procedure: there were 6 parameters in the MTM with a longitudinal binary outcome compared to the 18 to 20 parameters in the MTMs with a longitudinal ordinal outcomes. Hence, the removal of one parameter, known to have no effect on the outcome, may have a greater impact on the model that estimates fewer parameters overall. When we imposed the PO assumption on the treatment effect, we observed additional improvements in efficiency. However, we also note that the PO assumption could lead to invalid inference if misspecified. Thus, when deciding on the form of the regression model, one needs to strike a balance between robustness and efficiency. While we focused on treatment effect estimation in a randomized clinical trial, the proposed approach can be applied to non-randomized studies where one is interested in an end of follow-up exposure effect while adjusting for baseline covariates. However, in these settings, the inclusion of time-varying exposures and/or confounders would require careful consideration.

In the presence of at least one absorbing state, we showed how MTMs can directly model cumulative incidences without having to condition on any post-randomization outcomes. The outcome-marginal and conditional response-dependence components of an MTM correspond to cumulative incidence and discrete hazard models, respectively, and, when interest is in a fully-marginal (i.e. marginal with respect of both outcomes and covariates) treatment effect, the efficiency of MTMs tracks with Cox proportional hazard models and discrete-time multi-state survival models. In contrast to survival analysis methods, first-order MTMs provide a valid way to directly estimate covariate-adjusted treatment effects on cumulative incidence.

The design of CLOVERS did not allow us to observe transitions from the discharged/home state to either the hospital state and/or the hospital/ICU state. While, for the purpose of the analyses in this paper we assumed that none of unobserved transition occurred, this assumption might not be realistic in scenarios where the probability of remaining in the same state from day to day is lower. In future work, we will explore how marginalized transition models with longitudinal ordinal outcomes can be extended to cases where we have interval censored data, for example, when we know participants are alive, but we do not know in which state they are in. The MTMs presented in this paper can be fitted using the R package OMTM1.

Supplemental Material

sj-pdf-1-smm-10.1177_09622802261462925 - Supplemental material for Modeling treatment effects on absorbing outcomes in clinical trials: Leveraging longitudinal and ordinal data for efficiency gains

Supplemental material, sj-pdf-1-smm-10.1177_09622802261462925 for Modeling treatment effects on absorbing outcomes in clinical trials: Leveraging longitudinal and ordinal data for efficiency gains by Chiara Di Gravio, Ran Tao, Yunbi Nam, Eric P Schmidt, Nathan Shapiro and Jonathan S Schildcrout in Statistical Methods in Medical Research

Footnotes

ORCID iD

Chiara Di Gravio

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interest

The authors have no conflicting interests to declare with respect to the research, authorship, and/or publication of this article.

Supplemental material

Supplemental material for this article is available online.

References

Self

Shotwell

Gibbs

, et al. Renin-angiotensin system modulation with synthetic angiotensin (1-7) and angiotensin II type 1 receptor-biased ligand in adults with COVID-19: two randomized clinical trials. JAMA 2023; 329: 1170–1182.

The National Heart, Lung, and Blood Institute Prevention and Early Treatment of Acute Lung Injury Clinical Trials Network. Early restrictive or liberal fluid management for sepsis-induced hypotension. N Engl J Med 2023; 388: 499–510.

Moskowitz

Shotwell

Gibbs

, et al. Oxygen-free days as an outcome measure in clinical trials of therapies for COVID-19 and other causes of new-onset hypoxemia. Chest 2022; 162: 804–814.

Schoenfeld

Bernard

ARDS

Network

. Statistical evaluation of ventilator-free days as an efficacy measure in clinical trials of treatments for acute respiratory distress syndrome. Crit Care Med 2002; 30: 1772–1777.

Yehya

Harhay

Curley

MAQ

, et al. Reappraisal of ventilator-free days in critical care research. Am J Respir Crit Care Med 2019; 200: 828–836.

Aalen

Cook

Røysland

. Does cox analysis of a randomized survival study yield a causal treatment effect?. Lifetime Data Anal 2015; 21: 579–593.

Hernán

. The hazards of hazard ratios. Epidemiology 2010; 21: 13–15.

Post

van den Heuvel

Putter

. The built-in selection bias of hazard ratios formalized using structural causal models. Lifetime Data Anal 2024; 30: 404–438.

Young

Stensrud

Tchetgen Tchetgen

, et al. A causal framework for classical statistical estimands in failure-time settings with competing events. Stat Med 2020; 39: 1199–1236.

10.

Lee

Daniels

. A class of Markov models for longitudinal ordinal data. Biometrics 2007; 63: 1060–1067.

11.

Lee

Daniels

. Marginalized models for longitudinal ordinal data with application to quality of life studies. Stat Med 2008; 27: 4359–4380.

12.

McCullagh

. Regression models for ordinal data. J R Stat Soc Ser B: Methodol 1980; 42: 109–142.

13.

Schildcrout

Harrell

Heagerty

, et al. Model-assisted analyses of longitudinal, ordinal outcomes with absorbing states. Stat Med 2022; 41: 2497–2512.

14.

Cox

. Regression models and life-tables. J R Stat Soc Ser B (Methodol) 1972; 34: 187–220.

15.

Green

Symons

. A comparison of the logistic risk function and the proportional hazards model in prospective epidemiologic studies. J Chronic Dis 1983; 36: 715–723.

16.

Montoya

Daza

Muñoz

, et al. Development of a novel formulation with hypericin to treat cutaneous leishmaniasis based on photodynamic therapy in in vitro and in vivo studies. Antimicrob Agents Chemother 2015; 59: 5804–5813.

17.

Bergsma

Croon

Hagenaars

. Marginal models for dependent, clustered, and longitudinal categorical data. New York: Springer, 2009.

18.

Diggle

Heagerty

Liang

, et al. Analysis of longitudinal data. Oxford University Press, 2002.

19.

Heagerty

Zeger

. Marginal regression models for clustered ordinal measurements. J Am Stat Assoc 1996; 91: 1024–1036.

20.

Zeger

Liang

Albert

. Models for longitudinal data: a generalized estimating equation approach. Biometrics 1988; 44: 1049–1060.

21.

Azzalini

. Logistic regression for autocorrelated data with application to repeated measures. Biometrika 1994; 81: 767–775.

22.

Heagerty

. Marginally specified logistic-normal models for longitudinal binary data. Biometrics 1999; 55: 688–698.

23.

Heagerty

. Marginalized transition models and likelihood inference for longitudinal categorical data. Biometrics 2002; 58: 342–351.

24.

Zeger

Liang

. Longitudinal data analysis for discrete and continuous outcomes. Biometrics 1986; 42: 121–130.

25.

Breslow

Clayton

. Approximate inference in generalized linear mixed models. J Am Stat Assoc 1993; 88: 9–25.

26.

Stiratelli

Laird

Ware

. Random-effects models for serial observations with binary response. Biometrics 1984; 40: 961–971.

27.

Liang

Zeger

. Longitudinal data analysis of continuous and discrete responses for pre-post designs. Sankhyā: The Indian J Stat Ser B 2000; 62: 134–148.

28.

. On efficiency of constrained longitudinal data analysis versus longitudinal analysis of covariance. Biometrics 2010; 66: 891–896.

29.

Agresti

. Categorical data analysis. Vol. 482. John Wiley & Sons, 2003.

30.

Peterson

Harrell

FE Jr

. Partial proportional odds models for ordinal response variables. J R Stat Soc: Ser C (Appl Stat) 1990; 39(2): 205–217.

31.

Cook

Lawless

. Multistate models for the analysis of life history data. Chapman & Hall/CRC, 2018.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.20 MB