Estimating the effect of PEG in ALS patients using observational data subject to censoring by death and missing outcomes

TL;DR

This paper develops a causal inference framework to estimate the effect of PEG on BMI in ALS patients using observational registry data, addressing confounding, censoring by death, and missing data.

Contribution

It introduces a novel integrated approach combining propensity scores, principal stratification, and missing data models for causal inference in ALS registry data.

Findings

PEG has a positive effect on BMI at 18 months post-treatment

The framework improves causal effect estimation in observational ALS data

Results provide more definitive evidence on PEG's impact than previous studies.

Abstract

Though they may offer valuable patient and disease information that is impossible to study in a randomized trial, clinical disease registries also require special care and attention in causal inference. Registry data may be incomplete, inconsistent, and subject to confounding. In this paper we aim to address several analytical issues in estimating treatment effects that plague clinical registries such as the Emory amyotrophic lateral sclerosis (ALS) Clinic Registry. When attempting to assess the effect of a surgical insertion of a percutaneous endoscopic gastrostomy (PEG) tube on body mass index (BMI) using the data from the ALS Clinic Registry, one must combat issues of confounding, censoring by death, and missing outcome data that have not been addressed in previous studies of PEG. We propose a causal inference framework for estimating the survivor average causal effect (SACE) of PEG, which incorporates a model for generalized propensity scores to correct for confounding by pre-treatment variables, a model for principal stratification to account for censoring by death, and a model for the missing data mechanism. Applying the proposed framework to the ALS Clinic Registry Data, our analysis shows that PEG has a positive SACE on BMI at month 18 post-baseline; our results likely offer more definitive answers regarding the effect of PEG than previous studies of PEG.

Tables

Table 1. Table 1: Individual likelihood by observed groups if G i subscript 𝐺 𝑖 G_{i} is known

	$G_i$
Observed
Group $O(Z,S,M)$	$LL$	$LD$	$DL$	$DD$
$O(1,1,0)$	$(1-{\varphi }_{LL,1,i})\times$	$(1-{\varphi }_{LD,1,i})\times$	-	-
	${\pi }_{LL,i}{f}_{LL,i}$	${\pi }_{LD,i}{f}_{LD,i}$
$O(1,1,1)$	${\varphi }_{LL,1,i}{\pi }_{LL,i}$	${\varphi }_{LD,1,i}{\pi }_{LD,i}$	-	-
$O(1,0,-)$	-	-	${\pi }_{DL,i}$	${\pi }_{DD,i}$
$O(0,1,0)$	$(1-{\varphi }_{LL,0,i})\times$	-	$(1-{\varphi }_{DL,0,i})\times$	-
	${\pi }_{LL,i}{f}_{LL,i}$		${\pi }_{DL,i}{f}_{DL,i}$
$O(0,1,1)$	${\varphi }_{LL,0,i}{\pi }_{LL,i}$	-	${\varphi }_{DL,0,i}{\pi }_{DL,i}$	-
$O(0,0,-)$	-	${\pi }_{LD,i}$	-	${\pi }_{DD,i}$

Table 2. Table 2: Deviance Information Criterion (DIC) for analysis under Assumption 3 (latent ignorable missingness), Assumption 4 (ignorable missingness) and MCAR, as d 𝑑 d , the degree of polynomial splines for propensity scores, varies.

	d
	0	1	2	3	4	5
Assumption 3	1702.55	1683.98	1668.20	1777.01	1658.00	2034.26
Assumption 4	1769.40	1722.59	1727.66	1733.74	1714.54	1719.51
MCAR	2715.46	2658.06	2655.16	2654.74	2653.52	3200.87

Table 3. Table 3: Estimates of SACE of PEG treatment (with 95% credible intervals) on BMI at 18 months post-baseline Assumption 3 (latent ignorable missingness), Assumption 4 (ignorable missingness) and MCAR with and without adjustment for propensity scores (PS) (N=815).

		Effect Estimates
		Latent Ignorable	Ignorable
		Missingness	Missingness	MCAR
With PS	Time to PEG	0.17	0.16	0.14
		(-0.06,0.39)	(-0.05,0.38)	(-0.01,0.29)
	Binary Indicator	3.09	3.17	2.28
		(0.01,5.83)	(0.14,6.24)	(0.79,3.77)
Without PS	Time to PEG	0.24	0.24	0.20
		(0.03,0.46)	(0.03,0.45)	(0.03,0.35)
	Binary Indicator	0.45	0.56	-0.05
		(-2.42,3.45)	(-2.19,3.31)	(-1.37,1.23)

Equations

$$f(Y_i|S_i,M_i,G_i,Z_i,{\mathbf{D}}_i)=f(Y_i|M_i,G_i,Z_i,{\mathbf{D}}_i)=f(Y_i|G_i,Z_i,{\mathbf{D}}_i)=f(Y_i(Z_i)|G_i,{\mathbf{D}}_i)$$

$$\prod _{i}f(Y_i,S_i,M_i,G_i,Z_i|{\mathbf{D}}_i)$$

$$\prod _{i}f(Y_i(Z_i)|G_i,P{S}_i)f(M_i|G_i,Z_i,{\mathbf{D}}_i)f(G_i|P{S}_i),$$

$$\prod _{i}f(Y_i(Z_i)|G_i,P{S}_i)f(M_i|Z_i,{\mathbf{D}}_i)f(G_i|P{S}_i).$$

$$\prod _{i}f(Y_i(Z_i)|G_i,P{S}_i)f(G_i|P{S}_i).$$

$${\pi }_g=P(G=g|{\mathbf{X}}_2)=\frac{\exp ({\mathbf{X}}_2^T{\mathbf{\alpha }}_g)}{{\sum }_{g^{\prime }}\exp \left({\mathbf{X}}_2^T{\mathbf{\alpha }}_{g^{\prime }}\right)}$$

$${\varphi }_{gz}=Pr(M=1|G=g,Z=z,S=1)=\frac{e^{{\mathbf{X}}_3{\mathbf{\theta }}_{g,z}}}{1+e^{{\mathbf{X}}_3{\mathbf{\theta }}_{g,z}}}$$

$$\prod _{i\in O(1,1,0)}\left\{(1-{\varphi }_{LL,1,i}){\pi }_{LL,i}N({\mathbf{X}}_{1,LL,i}{\eta }_{LL},{\sigma }_{LL}^2)+(1-{\varphi }_{LD,1,i}){\pi }_{LD,i}N({\mathbf{X}}_{1,LD,i}{\eta }_{LD},{\sigma }_{LD}^2)\right\}\times$$

$$\prod _{i\in O(0,1,0)}\left\{(1-{\varphi }_{LL,0,i}){\pi }_{LL,i}N({\mathbf{X}}_{1,LL,i}{\eta }_{LL},{\sigma }_{LL}^2)+(1-{\varphi }_{DL,0,i}){\pi }_{DL,i}N({\mathbf{X}}_{1,DL,i}{\eta }_{DL},{\sigma }_{DL}^2)\right\}\times$$

$$\prod _{i\in O(1,1,1)}\left\{({\varphi }_{LL,1,i}){\pi }_{LL,i}+({\varphi }_{LD,1,i}){\pi }_{LD,i}\right\}\times$$

$$\prod _{i\in O(0,1,1)}\left\{({\varphi }_{LL,0,i}){\pi }_{LL,i}+({\varphi }_{DL,0,i}){\pi }_{DL,i}\right\}\times$$

$$\prod _{i\in O(1,0,-)}\left\{{\pi }_{DD,i}+{\pi }_{DL,i}\right\}\prod _{i\in O(0,0,-)}\left\{{\pi }_{DD,i}+{\pi }_{LD,i}\right\}$$

Full text

Keywords: Censoring by death; Confounding; Generalized propensity scores; Missing data; Observational study; Principal stratification.

1 Introduction

While observational data from disease registries can be used to assess treatment effects, they are often fraught with analytical challenges. This is particularly true of a disease with high disability and mortality rates as issues of unmeasured data often occur due to a patient’s absence or death. Our work is motivated by the analysis of data from an amyotrophic lateral sclerosis (ALS) Clinic Registry, aimed at estimating the effect of a palliative treatment, percutaneous endoscopic gastrostomy (PEG) on body mass index (BMI), a measure of nutritional management and quality of life. ALS is a neurodegenerative disorder with a very poor prognosis (Gelinas and Miller, 2000; Procaccini and Nemergut, 2008). Currently there is no known cure for ALS and existing treatment procedures for ALS patients are mostly palliative; there has been substantial interest in assessing effects of these palliative procedures on quality of life (Miller et al., 2009).

Dysphagia, or difficulty in swallowing, affects almost all patients with ALS, and subsequently along with muscle atrophy and hypermetabolism causes malnutrition amongst the ALS patient population. Nutritional management is key in disease management and palliative care (Desport et al., 2006; Muscaritoli et al., 2012). PEG, a surgically inserted tube providing enteral nutrition, is generally considered when an ALS patient’s nutritional status deteriorates and weight loss is greater than 10% of the baseline weight (Goyal and Mozaffar, 2014; Park et al., 1992), though patients may elect not to receive PEG for other reasons. While some studies have shown that PEG tube insertion offers some benefits in stabilizing or even increasing body weight, the overall results regarding the effect of PEG on weight or BMI are inconclusive, and in particular the long term benefits of the procedure are uncertain (Kasarskis et al., 1996; Katzberg and Benatar, 2011; Mazzini et al., 1995). One challenge for assessing the effect of PEG in ALS patients is that randomized studies are not ethically feasible, since PEG placement has long been recommended as standard of care in ALS clinics for patients with nutritional compromise (Miller et al., 2009). As a result, the effect of PEG has been evaluated using observational data, and previous studies did not adequately address analytical issues including confounding and missing data that were commonly encountered. In addition, the existing studies had small to moderate sample sizes, limiting the generalizability of their findings.

The Emory ALS Clinic Registry has several notable strengths, making it well-suited for assessing the effect of PEG. It includes more than 1000 patients with long-term follow-up and a wide range of relevant clinical variables such as BMI, vital capacity (VC), forced vital capacity (FVC), negative inspiratory force (NIF), and ALS Functional Rating Scale Revised score (ALSFRS-R). In particular, ALSFRS-R is a measure of disease progression and is not always collected in large ALS registries. However, the data present several analytical challenges to estimating a causal effect of PEG. First, insertion of PEG is not randomized, which may lead to confounding by pre-treatment variables. Second, due to the fatal and fast-progressing nature of ALS outcome measurements may be “censored by death,” i.e, a patient dies before the outcome of interest can be measured (Zhang and Rubin, 2003; Rubin, 2006). Third, for those who are not “censored by death,” the outcome may be missing due to various reasons.

In the context of Rubin’s Causal Model, principal stratification, first described in detail in Frangakis and Rubin (2002), offers a framework for causal inference in the presence of a confounding post-treatment variable, such as censoring by death. Subsequently, Zhang and Rubin (2003) extended this approach to cases where a post-treatment variable such as survival or graduation “censors” the outcome of interest. Zhang et al. (2009) further outlined specific modeling approaches for the identification of survivor average causal effect (SACE) within the principal stratification framework. These existing works in this area, however, analyzed data from randomized studies, assuring complete collection of outcome data and, more importantly, avoiding the issues of confounding, an advantage not guaranteed in our retrospective clinic registry data.

More recently, Frumento et al. (2012) addressed post-treatment variables and missing outcome data simultaneously within a principal stratification framework in a randomized study. However, the Emory ALS Clinic Registry data introduce some unique challenges that were not present in the analysis of Frumento et al. (2012) and require new analytical approaches. In particular, our study treatment is not randomized and includes a component of time to treatment, leading to potential confounding by pre-treatment variables. Furthermore, both treatment and outcome measurement are subject to censoring by death in our study. Finally, missing outcome data in our study are due to the observational nature of the study, and have different implications compared to the missing outcomes in Frumento et al. (2012). Therefore, while some of the same tools are utilized in our work, the proposed framework addresses important new issues.

It is particularly important to address potential confounding by pre-treatment variables in a principal stratification model, because Frangakis and Rubin (2002) defined principal stratification in a randomized experiment. When such confounding may be present, either as residual confounding in a randomized clinical trial or due to observational data, Schwartz et al. (2012) showed that the resulting principal effect estimate is likely to be biased. Of the many methods for handling confounding, a popular choice is the propensity score introduced by Rosenbaum and Rubin (1983). The propensity score provides a means of balancing pre-treatment variables across treatment groups, a result that would otherwise be guaranteed if a randomized study design is used. Though propensity scores were initially introduced for balancing treatment assignment groups when treatment is binary, other authors have extended these methods to non-binary treatment assignment models, such as generalized propensity score methods (Imai and Van Dyk, 2004; Hirano and Imbens, 2004). These methods allow for balancing pre-treatment variables when treatment assignment is categorical or continuous. A recent work by Jo and Stuart (2009) used propensity scores within a principal stratification framework, not for the removal of confounding as the study was randomized, but instead to predict principal strata membership in a matched analysis, again different from our setting. Furthermore, propensity score methods for non-binary treatment assignment have not yet been employed for conditional ignorability in a principal stratification framework when outcomes may be missing.

To account for aforementioned complicating analytical issues, we develop a causal inference approach for estimating the survivor average causal effect (SACE) of PEG. Building on the principal stratification framework to account for censoring by death (i.e., confounding by a post-treatment variable), our approach uses generalized propensity scores to correct for confounding by pre-treatment variables and it also accounts for missing outcome data by incorporating a model for the missing data mechanism. Although it is developed for the specific application, our approach can be applied to analysis of other observational studies that are subject to similar issues. The remainder of the article is organized as follows. First, we introduce the data and notation in Section 2. We present our framework for causal inference in Section 3 and a Bayesian inference procedure in Section 4. In Section 5, we apply the proposed approach to the Emory ALS data. We conclude this paper with some discussion remarks in Section 6.

2 Data Structure in the ALS Registry

In the analysis of the ALS Clinic Registry, our goal is to assess the effect of PEG on BMI (denoted by $Y$) at $t^{o}=18$ months post-baseline; without loss of generality we assume that $Y\in\mathcal{R}$. We denote by ${\bf D}$ the set of patient baseline characteristics and by $T_{S}$ the survival time post-baseline confirmed by Social Security database records, noting that all patients in our data set died. The surgical insertion of a PEG tube may be administered at any time post-baseline given patient choice. Let $Z$ be the binary indicator such that $Z=1$ if a patient received PEG prior to the time of outcome measurement $t^{o}$ and $Z=0$ if otherwise; let $T_{Z}$ denote time to treatment. It follows that patients untreated prior to $t^{o}$ have values $Z=0$ and $T_{Z}=T_{S}$ if they did not survive until $t^{o}$, or $Z=0$ and $T_{Z}=t^{o}$ if they survived past $t^{o}$. For those who died prior to $t^{o}$, $Y$ is undefined and is known as “censored by death.” For those who are alive at $t^{o}$, $Y$ is defined but may be missing due to other reasons and let $M$ be the binary indicator such that $M=1$ if $Y$ is missing and $M=0$ if $Y$ is observed. Of note, for patients with $S=0$ who were “censored by death,” both $Y$ and $M$ are undefined, denoted by “*”, following the notation of Zhang and Rubin (2003). This extends the sample space for $Y$ to {$\mathcal{R},*$} and for $M$ to {$0,1,*$}.

3 Causal Inference Framework

As earlier defined, the binary indicator for treatment from baseline until $t^{o}$ (time of outcome measurement) is $Z$ and the outcome of interest is $Y$. A post-treatment indicator of survival past the $t^{o}$ is $S$. Following the Rubin Causal Model (Holland, 1986), we define potential outcomes $\mathcal{Y}=\{\left(Y_{i}(0),Y_{i}(1)\right)\mbox{ for }i=1\ldots n\}$ and $\mathcal{S}=\{\left(S_{i}(0),S_{i}(1)\right)\mbox{ for }i=1...n\}$ under the two treatment options $Z=0$ and 1. For simplicity, we suppress $i$ in what follows if it does not create ambiguity.

Figure 1 illustrates the potential outcomes in the Emory ALS data. For example, patients who are treated prior to time of outcome measurement ($T_{Z}>t^{o}$) may follow one of two paths: survival until $t^{o}$ ($S$ = 1) or death prior to $t^{o}$ ($S$ = 0). Those individuals who survive may have observed $Y$, as designated by the darkest shaded box in Figure 1, or have missing $Y$ as designated in the medium shaded box. Patients who die before $t^{o}$ have both $Y$ and $M$ undefined, i.e., censored by death, as designated by the lightest shaded box.

In the presence of censoring by death, we use post treatment survival status to define four potential principal strata ($G$) by $\left(S_{i}(0),S_{i}(1)\right)$, namely, $LL$, $LD$, $DL$, and $DD$. Individuals who are in the $LL$ stratum are those who would be alive at time $t^{o}$ regardless of treatment. Those individuals who are in the $LD$ stratum would be alive if they received the treatment but would not be alive if they did not and those who are in the $DL$ stratum would not be alive if they receive the treatment, but would be alive if they did not. Finally, individuals in the $DD$ are those who would die prior to time $t^{o}$ regardless of treatment. A causal effect of treatment is only estimable in strata where individuals would be alive and therefore have a measurable outcome, regardless of treatment status. As such, the causal effect is estimable in the $LL$ stratum only. Specifically, the Survivor Average Causal Effect (SACE) is defined as the difference in $Y$ between treated and untreated individuals in the $LL$ stratum.

However, we do not get to observe $G$ since we only know the survival status given the observed treatment status. Based on $(Z,S)$, there are four observed groups denoted by $O(Z,S)$: $O(1,1)$, $O(1,0)$, $O(0,1)$, and $O(0,0)$. Each of these observed groups is a mixture of the principal strata. Specifically, $O$(1,1) is a mixture of individuals from the $LL$ and $LD$ strata, $O$(1,0) is comprised of individuals from the $LL$ and $DL$ strata, $O$(0,1) is comprised of individuals from the $DD$ and $DL$ strata, and $O$(0,0) is comprised of individuals from the $DD$ and $LD$ strata. In addition, based on $(Z,S,M)$, we have six observed groups denoted by $O$($Z,S,M$), namely, $O$(1,1,1), $O$(1,1,0), $O$(1,0,-), $O$(0,1,1), $O$(0,1,0), $O$(0,0,-); see Table LABEL:t3tab:like. The two original observed groups with $S=1$, namely, $O$(1,1) and $O$(0,1), are each divided into two subgroups based on whether $Y$ is observed or missing. The groups with $S=0$ do not need to be divided further, since $M$ is undefined.

3.1 Assumptions

In our causal inference framework, we adopt two standard assumptions, the Stable Unit Treatment Value Assumption (SUTVA) as defined by Cox in 1958 and further elaborated in Rubin (1980) and the Strong Ignorability of Treatment Assignment (Rosenbaum and Rubin, 1983).

Assumption 1

[SUTVA] The potential outcomes of any unit is independent of the treatment status of the other units given ${\bf D}$.

Assumption 2

[Strong Ignorability of Treatment Assignment] $\left(Y_{i}(0),Y_{i}(1)\right)$ and $\left(S_{i}(0),S_{i}(1)\right)$ (or equivalently $G$) are independent of $Z$ given ${\bf D}$.

Regarding the missing data mechanism for $Y$, we consider two assumptions, a latent ignorability (Frangakis and Rubin, 1999) and a standard ignorability (Little and Rubin, 2002).

Assumption 3

[Latent Ignorable Missingness] $M$ does not depend on the missing values or parameters of the data distribution conditional on ${\bf D}$ and $G$; particularly, $M$ is independent of $Y$ given ${\bf D}$ and $G$.

Assumption 4

[Ignorable Missingness] $M$ does not depend on the missing values or parameters of the data distribution conditional on ${\bf D}$; particularly, $M$ is independent of $Y$ given ${\bf D}$.

Since $G$ is not completely known, the latent ignorable missingness in Assumption 3 is not ignorable and it is necessary to model $M$ in the modeling framework.

3.2 Data Likelihood

We define the complete data as $(Y_{i},S_{i},M_{i},G_{i},Z_{i},{\bf D}_{i})$ for $i=1,\ldots,n$, which under Assumption 1 constitute independent and identically distributed observations. Accordingly, the observed data are represented as $(Y_{i},S_{i},M_{i},G_{obs,i},Z_{i},{\bf D}_{i})$, where $G_{obs,i}$ is the set of feasible principal strata that an individual may belong to given $Z_{i}$ and $S_{i}$.

We first relate the distribution of the complete data $(Y_{i},S_{i},M_{i},G_{i},Z_{i},{\bf D}_{i})$ to the distribution of the potential outcomes under Assumptions 1-3. We have

[TABLE]

where the first equality is due to the fact that $S_{i}$ is redundant given $Z_{i}$ and $G_{i}$, the second equality is due to Assumption 3, and the third equality is due to Assumption 2. In addition, $f(S_{i}|M_{i},G_{i},Z_{i},{\bf D}_{i})=1$, since $S_{i}$ is fully determined by $Z_{i}$ and $G_{i}$. From Assumption 2, we have $f(G_{i}|Z_{i},{\bf D}_{i})=f(G_{i}|{\bf D}_{i})$. Using these facts and Assumption 1, the complete data likelihood can be written as

[TABLE]

where $f(Y_{i}(Z_{i})|G_{i},{\bf D}_{i})$ models the potential outcomes, $f(M_{i}|G_{i},{\bf D}_{i})$ models the missing data mechanism, $f(G_{i}|{\bf D}_{i})$ models the principal strata, and $f(Z_{i}|{\bf D}_{i})$ models the treatment assignment. Without loss of clarity, we define $f(Y_{i}|G_{i},{\bf D}_{i})=1$ if $Y_{i}$ is missing; in other words, only observed $Y_{i}$ contributes to the first part of the complete data likelihood (1). Assuming the parameters in $f(Z_{i}|{\bf D}_{i})$ are distinct from the other parameters and $f(Z_{i}=z|{\bf D}_{i})>0$ for $z=0,1$, we can drop $f(Z_{i}|{\bf D}_{i})$ from the complete data likelihood.

The observed data likelihood is considerably more complicated, since it involves mixture distributions since we only observe $G_{obs,i}$ and do not observe $G_{i}$. As shown in Section 4.4, a data augmentation algorithm allows us to work with the complete data likelihood directly and hence avoid the complicated observed data likelihood.

Since the treatment of interest includes a time-to-intervention component $T_{Z}$, we use generalized propensity scores in the spirit of Imai and Van Dyk (2004), Lu (2005) and Hirano and Imbens (2004) to further accommodate this unique feature. Specifically, we model $T_{Z}$ which may be censored by $T_{S}$ using a Cox proportional hazards (PH) model $h(t)=h_{0}(t)\exp\left({\boldsymbol{\beta}}^{T}{\bf D}\right)$ and use its linear predictor ${\boldsymbol{\beta}}^{T}{\bf D}$ as the generalized propensity score, denoted by $PS$. Along the lines of Imai and Van Dyk (2004) and Lu (2005), it follows from Assumptions 1 and 2 that $Z$ is independent of $\left(Y_{i}(0),Y_{i}(1)\right)$ and of $\left(S_{i}(0),S_{i}(1)\right)$ conditional on $PS_{i}$. Then the complete data likelihood (1) can be rewritten as

[TABLE]

where $f(Z_{i}|{\bf D}_{i})$ is dropped.

Alternatively, under Assumptions 1, 2 and 4, i.e., replacing the latent ignorable missingness with the ignorable missingness, the complete data likelihood (2) can be simplified to

[TABLE]

Assuming that its parameters are distinct from parameters in other models, $f(M_{i}|Z_{i},{\bf D}_{i})$ can be dropped from the complete data likelihood, leading to

[TABLE]

In Section 5, we analyze the ALS data using both (2) and (3) under the latent ignorable missingness and the ignorable missingness, respectively.

4 Modeling Framework

In this section, we specify the model for each component of the complete data likelihood (2) and then describe a Bayesian inference approach under Assumptions 1-3. The corresponding inference approach under Assumptions 1-2 and 4 is simpler and can be derived along similar lines.

4.1 Model for Potential Outcomes

For $f(Y_{i}(Z_{i})|G_{i},PS)$ in (2), we assume that $\left(Y(1),Y(2)\right)$ within each principal stratum $G=g$ have a normal distribution, $f_{g}$, with parameters ${\boldsymbol{\eta}}_{g}$ and covariates ${\bf X}_{1,g}$ that may differ by strata. Also, individual likelihoods for $Y$ are dependent on potential strata $g$, and are represented by $f_{g,i}$. Specifically, $Y\sim N({\bf X}_{1,g}{\boldsymbol{\eta}}_{g},\sigma_{g}^{2})$ for $g\in{LL,LD,DL}$, where ${\bf X}_{1,LL}$ includes the column for intercept, $Z$, $T_{Z}$, and $PS$. ${\bf X}_{1,LD}$ and ${\bf X}_{1,DL}$ include intercept and $PS$ but not $Z$ and $T_{Z}$ as patients with an observed outcome in each of these strata are either all treated or all untreated respectively. The regression coefficient for $Z$ in the $LL$ stratum represents SACE and is of primary interest. Inclusion of $T_{Z}$ allows us to assess the relationship between $Y$ and time to treatment within those who receive the treatment.

4.2 Model for Principal Stratification

The probabilities of the four principal strata, denoted by $(\pi_{LL},\pi_{LD},\pi_{DL},\pi_{DD})$, are modeled using a multinomial logit model with a set of covariates denoted by $\bf{X}_{2}$ including $PS$. he probability of a patient being in principal stratum $g$ is given in equation (4), and individual probabilities of principal strata are represented by $\pi_{g,i}$. As in any multinomial logit model, one category must be selected as a reference group.

[TABLE]

Of the assumptions available for sharpening the bounds of effect estimates in the principal stratification framework, a monotonicity assumption is often used, i.e., one of the principal strata, say the $DL$ group does not exist Zhang and Rubin (2003). While this monotonicity assumption may be plausible if treatment is always beneficial on patient survival, it is also possible that a subpopulation of weak patients would be not be able to withstand treatment or recovery and may die shortly after treatment, but could stay alive if left untreated, when evaluating an intrusive treatment. For such a data analysis, it may be appropriate to allow all four principal strata to exist and to relax the assumption of monotonicity. In the ALS data, the SACE for PEG is estimated both with and without the monotonicity assumption to assess its impact.

4.3 Model for Missing Data Mechanism

As discussed before, $M$ for patients with $S=0$ is undefined, so it is only valid to model $M$ for those individuals in the $LL$ stratum with treatment, the $LL$ stratum without treatment, the $LD$ stratum with treatment, and the $DL$ stratum without treatment. The model for $M$ is presented in (5).

[TABLE]

where ${\bf X}_{3}$ is a set of covariates associated with the missing data mechanism. The individual probability of missingness is dependent on strata and treatment value and is represented by $\phi_{gz,i}$. When Assumption 4 is used instead of Assumption 3, this model is not needed.

4.4 Bayesian Inference

Given the models specified in Section 4.1-4.3, the individual likelihoods for all possible combinations of $Z$, $S$, and $M$ are given in Table LABEL:t3tab:like, where each cell value is the likelihood of the complete data if $G$ is known. The conditional probability of $\mathnormal{G_{i}=g}$ given the observed data is the ratio of each cell to the total of that row. Rows $O(1,0,-)$ and $O(0,0,-)$ are included in Table LABEL:t3tab:like for a comprehensive understanding of the possible combinations of treatment and survival, but individuals who fall into these groups do not have outcome data and hence do not contribute to the observed data likelihood. As such, these individuals will only contribute to the model for the probability of principal strata, which is reflected in the observed data likelihood in equation (6).

[TABLE]

Prior distributions for the specified parameters in the observed data likelihood should be chosen carefully, with thought to distributions that may be informative, proper, and conjugate where appropriate. For this analysis, a conjugate multivariate normal prior distribution is assumed for ${\boldsymbol{\eta}}_{g}$, ${\boldsymbol{\eta}}_{g}\sim Normal_{p}\left(\mu_{g},\sigma_{g}^{2}V_{g}\right)$. Similarly, a conjugate prior distribution is assumed for $\sigma_{g}^{2}$, $\sigma_{g}^{2}\sim InverseGamma\left(\nu_{g},\omega_{g}\right)$. We choose diffused priors for ${\boldsymbol{\alpha}}_{g}$ and ${\boldsymbol{\theta}}_{g,z}$.

Though the principal stratum of each individual is unknown, the observed treatment and survival groups may be used to inform imputation of the principal stratum assignments. Applying this idea, we use the Data Augmentation (DA) algorithm (Tanner and Wong, 1987) for posterior computation, in which information about the latent groups, namely, $G$, is imputed and subsequently the posterior parameters distributions are simulated.

The DA algorithm includes two iterative and alternating steps to allow for posterior inference. The first step, the Imputation or I-step, imputes the value of the principal stratum $G$ for each individual. This is accomplished by using the parameter values ${\boldsymbol{\alpha}}_{g}^{\left(k\right)}$, ${\boldsymbol{\eta}}_{g}^{\left(k\right)}$, $\sigma_{g}^{2\left(k\right)}$, and ${\boldsymbol{\theta}}_{g,z}^{\left(k\right)}$ (for approach 2 only) from the current approximation of posterior (from the $k$th iteration) to generate $G^{\left(k+1\right)}$ by using the conditional probabilities that are given by taking the ratio of each cell value to the row total in Table LABEL:t3tab:like. These conditional probabilities, $\rho_{O}$, are used in a Bernoulli distribution that imputes individual memberships to one of the two principal strata that correspond with the observed group $O$. More specifically, at the ($k$+1) iteration, each individual has a probability of being in a stratum that depends on their observed values ($Z$, $S$, $M$, $Y$, $PS$).

The P-step, or Posterior step, is then employed by using the imputed complete data set, and the parameters $\{{\boldsymbol{\theta}}^{\left(k\right)},\left(\pi_{g}^{\left(k\right)},{\boldsymbol{\eta}}_{g}^{\left(k\right)},\sigma_{g}^{2\left(k\right)}\right)\}$ can be updated to $\{{\boldsymbol{\theta}}^{\left(k+1\right)},\left(\pi_{g}^{\left(k+1\right)},{\boldsymbol{\eta}}_{g}^{\left(k+1\right)},\sigma_{g}^{2\left(k+1\right)}\right)\}$ by sampling from the full conditional distributions of each parameter based on the complete data likelihood. Either a Gibbs Sampler or a Metropolis-Hastings (MH) Algorithm can be employed for sampling. The details of the DA algorithm are provided in Web Appendices A and B.

If we use Assumption 4 instead of Assumption 3, we can drop all terms involving $\phi_{g,z}$ from Table LABEL:t3tab:like and use the complete data likelihood (3), leading to a simpler version of the aforementioned DA algorithm.

5 Analysis of ALS Clinic Data

We apply the proposed approach to analysis of the data from patients who visited the Emory ALS clinic at least once between 1997 and 2014. Patients are excluded from the analysis for not having any follow up clinic visits from baseline to outcome measurement at month 18 post-baseline or for having long survival times ($>$5 years post-baseline), resulting in a cohort of $N=815$ for our analyses. It has been noted that ALS patients with long survival times are likely different from the other ALS patients (Mateen et al., 2010), which may not be fully explained by clinical variables collected. Characteristics measured at baseline for each individual include sex, race, age, site of disease onset, months from symptom onset to baseline, months from diagnosis to baseline, body mass index (BMI), vital capacity (VC), forced vital capacity (FVC), negative inspiratory force (NIF), and ALS Functional Rating Scale Revised score (ALSFRS-R). Continuous variables, namely age, months from symptom onset to baseline, months from diagnosis to baseline BMI, VC, FVC, NIF, and ALSFRS-R, are standardized before being included as covariates for the propensity score model and the modesls for principal stratification model, missing data mechanism and outcomes. Since our focus is to handle missing outcome values, we impute missing values in covariates prior to subsequent analysis using the proposed approach.

Of the 275 treated patients, 32% or 89 individuals are alive 18 months from baseline, while of the 540 untreated individuals, 34% or 186 individuals are alive at this time-point (p = 0.606). Baseline measurements of BMI, FVC, and VC, as well as age at diagnosis, survival until $t^{o}$, and proportion of white versus other races are not significantly different between treated and untreated groups. However, clinically relevant differences between treatment groups do exist. NIF and ALSFRS-R scores at baseline are significantly higher for treated versus untreated individuals at a level of $\alpha=0.05$. Also the proportions of female individuals and of individuals with spinal disease are significantly differnt for treated individuals than untreated individuals (p = 0.040 and p$<$ 0.001). Most interestingly, the naïve analysis that compares the outcome of interest, BMI at 18 months post-baseline, across treatment groups indicates that treated individuals have a BMI that is 1 unit lower than the untreated group, and this difference in significant (p = 0.021). These results suggest that the effect of PEG may be confounded by variables that are different between treated and untreated groups, hence it is important to use the generalized propensity score. A table of patient characteristics by treatment status is available in the online Appendix (Web Table 1).

The 275 survivors at $t^{o}$ are significantly different from 540 non-surviving individuals by nearly every clinical comparison, with only demographic characteristics of sex and race being similar amongst the two groups, suggesting the need to adjust for the post-treatment variable, censoring by death. Additionally there are significant differences amongst patients with observed outcome and those with missing outcome data in, for example, BMI, FVC, VC, NIF, and ALS FRS-R score at baseline. indicating that the missing completely at random (MCAR) assumption is unrealistic and necessitates the use of the proposed approach for handling missing data. Web Tables 2 and 3 providing support for these statements are available in the online Appendix.

The SACE of PEG treatment is estimated using three modeling approaches. In addition to the two approaches using Assumptions 3 (latent ignorable missingness) and 4 (ignorable missingness) as described in Section 4.3, we consider a third approach that includes only individuals with an outcome measurement at $t^{o}$, analogous to a complete-case analysis that essentially assumes that the outcomes are MCAR. We also apply the same set of methods without the propensity score adjustment as well as with a monotonicity assumption, i.e., removal of the $DL$ stratum (Zhang and Rubin, 2003). Finally, to allow for flexibility, polynomial splines are used for incorporating the generalized propensity score in the models with $d$ denoting the degree of the polynomial spline. For all analyses, the MCMC algorithm is run for a total of 5,000 iterations, with a burn-in period of 3000 iterations.

To identify the model that best fits the data, we use the deviance information criterion (DIC) (Spiegelhalter et al., 2002). Models with smaller values of DIC are preferred to models with larger values, and the best fit is indicated by a minimum DIC. Table 2 presents DIC for each of the modeling approaches under different assumptions for missing data mechanism with various degrees of polynomial splines for the propensity score. The models for all three approaches with $d=4$ have the lowest DIC, and therefore offer the best fit for each modelling approach. Particularly, the model assuming latent ignorable missingness with $d=4$ has the lowest DIC among all models.

The top two rows of Table 3 present a comparison of the different modeling approaches when propensity scores are included. Specifically, the results from the models with $d=4$ are presented for each modelling approach, as these have been chosen best fitting models according to DIC. In contrast to the naïve analysis in Web Table 1 that shows a negative effect of PEG, analyses under both Assumption 3 and Assumption 4 show that there is a significant, positive SACE of PEG surgery on the outcome of BMI at 18 months post-baseline. Holding all else constant and assuming PEG insertion at or immediately after baseline, BMI at 18 months increases by about 3 units for those individuals who receive treatment when compared to those who are not treated among survivors. The effect of time to treatment is also positive but not statistically significant. Noting that a higher time to PEG ($T_{Z}$) leads to a shorter time from PEG to the time of outcome measurement, these results seem to suggest that the positive effect of PEG on BMI decreases gradually after the treatment, i.e., the possibility for a waning long-term effect. Additionally, the approach assuming MCAR indicates that the estimated effects of time to treatment and binary treatment are smaller in magnitude than those estimated from the two approaches under Assumptions 3 and 4, respectively.

The lower half of Table 3 presents the results from the three approaches without adjustment for propensity scores in the estimation of the SACE of PEG treatment, showing some changes in the magnitude of the effect estimates. On the one hand, removing the propensity score from the models in Approaches 1 and 2, the time to treatment effect remains positive and is now significant. On the other hand, in these same models, the coefficient estimate for the binary treatment indicator is drastically smaller in magnitude and no longer statistically significant. In the model that assumes MCAR and does not use propensity scores, the coefficient estimate for the binary treatment indicator changes direction and is non-significant.

In addition, the results, not provided, are not sensitive to the assumption of monotonicity. This may be due to the small proportion of individuals in $DL$ strata when all four strata are considered. When monotonicity is not assumed, most patients are in the $LL$ and $DD$ strata, with a $LD$ and $DL$ strata comprising less than 20% of the individuals in total. It is conceivable then that reallocating such a small proportion of individuals when removing the $DL$ stratum would not substantially change the effect estimates of the other strata.

Overall, the results in our analysis demonstrate a positive SACE of PEG, in contrast to inconclusive results in magnitude, direction, and significance from previous observational studies of PEG in association with weight or BMI. Given that our methodology is developed to account for the many complicating issues of an observational study that were not addressed in previous studies, it is likely that the positive effect of PEG estimated in this analysis is more definitive. As data collection is still ongoing in the ALS Clinic Registry, the proposed approach can be applied to this registry in the future, potentially further validating our current findings.

6 Discussion

Our current work represents the very first attempt to leverage the rich data in the Emory ALS Clinic Registry for assessing the effect of PEG in a large ALS population. Our proposed framework for causal inference addresses several challenges in the analysis of this data that have not been accounted for in previous studies of PEG. Our results show a positive effect of PEG placement on survivors that may be waning over time after insertion. Our approach can be directly applied or extended to evaluate effects of PEG or other palliative procedures such as the non-invasive positive pressure ventilation (NIPPV) (Miller et al., 2009) on additional outcomes that are of interest to clinicians including quality of life and disease progression.

One limitation of our analysis is that we do not have data on the usage of PEG tube after its insertion. Of note, patients with PEG may elect not to use PEG for nutrition delivery from time to time, thus the level of usage, similar to treatment adherence, could have an impact on patient outcomes such as BMI. Additionally, while the use of propensity scores within the principal stratification framework allows for the estimation of an unbiased principal effect using observational data, the reliability of removing bias via propensity score adjustment is predicated on the assumption of strongly ignorable treatment assignment, which means there must be no unmeasured confounders. Partly based on our experience from this work, a future prospective observational study, currently in the design and planning stage, will collect usage data for palliative procedures such as PEG and NIPPV as well as additional potential confounders that are not collected or validated in standard clinical settings. In addition, future applications of the proposed methods to other data with richer measurements of confounders should further demonstrate the reduction of confounding.

Making adjustments for missing outcome data in the context of causal inference typically requires strong assumptions about the ignorability of the missing data mechanism and creativity in the modeling framework. In our data application, the similar results under Assumptions 3 and 4 indicate that the findings are not sensitive to the latent ignorable assumption for this data, though we cannot empirically test this assumption. When in doubt, further stratification by missingness of outcome data and within the principal stratification framework may allow for additional flexibility in the assumptions imposed on the missing data mechanism.

Future consideration may also be given to jointly modeling the propensity score with the outcome model, missing data model, and principal strata model in a Bayesian framework, similar in spirit with Zigler et al. (2013) and Zigler and Dominici (2014). This would allow the quantities observed by in each of these three models to affect the posterior of propensity score in each MCMC iteration. While this could provide a more robust propensity score adjustment, Zigler et al. (2013) showed that the feedback between model stages in joint modeling can cause biased causal effect estimates if individual covariates are not also adjusted for in the outcome model. This bias should be accounted for if one uses joint modeling of the four models of outcome, missing data mechanism, propensity scores, and principal strata.