Prediction of Individual Outcomes for Asthma Sufferers

TL;DR

This paper develops a longitudinal model to predict individual asthma outcomes and optimize medication recommendations, using US survey data and latent variable techniques to improve clinical decision-making.

Contribution

It introduces a novel latent process model for medication effects that accounts for prescription refill data, enhancing personalized asthma treatment predictions.

Findings

Effective prediction of adverse asthma events

Reduced bias in medication effect estimation

Framework generalizable to other conditions

Abstract

We consider the problem of individual-specific medication level recommendation (initiation, removal, increase, or decrease) for asthma sufferers. Asthma is one of the most common chronic diseases in both adults and children, affecting 8% of the US population and costing $37-63 billion/year in the US. Asthma is a complex disease, whose symptoms may wax and wane, making it difficult for clinicians to predict outcomes and prognosis. Improved ability to predict prognosis can inform decision making and may promote conversations between clinician and provider around optimizing medication therapy. Data from the US Medical Expenditure Panel Survey (MEPS) years 2000-2010 were used to fit a longitudinal model for a multivariate response of adverse events (Emergency Department or In-patient visits, excessive rescue inhaler use, and oral steroid use). To reduce bias in the estimation of medication effects, medication level was treated as a latent process which was restricted to be consistent with prescription refill data. This approach is demonstrated to be effective in the MEPS cohort via predictions on a validation hold out set and a synthetic data simulation study. This framework can be easily generalized to medication decisions for other conditions as well.

Tables

Table 1. Table 1 : Summary of the model for asthma related events defined in ( 2 )  − -  ( 6 ) and the specification of prior distributions.

Description	Model	Prior Distributions	Specification
Individual Effects	${𝜶}_i\stackrel{iid}{\sim }N(𝝂,𝚺)$, $i=1,\dots ,n$	${\nu }_l\stackrel{ind}{\sim }N(M_{\nu ,l},S_{\nu ,l}^2)$	$M_{\nu ,l}=0$
			$S_{\nu ,l}^2=100$
		$𝚺\sim \text{Wishart}({𝑷}_{\Sigma },N_{\Sigma })$	${𝑷}_{\Sigma }=𝑰$
			$N_{\Sigma }=4$
Fixed Effects	${\beta }_{j,l}$ in (3)	${\beta }_{j,l}\stackrel{ind}{\sim }N(B_{j,l},T_{j,l}^2)$	$B_{j,l}=0$
			$T_{j,l}^2=100$
Medication Effects	${\gamma }_l$ in (3)	$-{\gamma }_l\stackrel{ind}{\sim }\log N(C_l,U_l^2)$	$C_l=-1$
			$U_l^2=10$
	${𝜹}_i\sim GP(0,K)$, $i=1,\dots ,n$	$e^{-\theta }\sim \text{Beta}(A_{\theta },B_{\theta })$	$A_{\theta }=1$
Individual Time			$B_{\theta }=1$
Effect	$K$ in (5)	$𝚽\sim \text{Wishart}({𝑷}_{\Phi },N_{\Phi })$	${𝑷}_{\Phi }=𝑰$
			$N_{\Phi }=4$
	$z_i$ in (6)	${\mu }_h\sim N(M_{\mu },S_{\mu }^2)$	$M_{\mu }=0$
			$S_{\mu }^2=50$
		${\sigma }^2\sim {\text{Gamma}}^{-1}(A_{\sigma },B_{\sigma })$	$A_{\sigma }=5$
Medication			$B_{\sigma }=.05$
Level Process		$e^{-\varphi }\sim \text{Beta}(A_{\varphi },B_{\varphi })$	$A_{\varphi }=6$
			$B_{\varphi }=3$
		$\rho \sim \text{Gamma}(A_{\rho },B_{\rho })$	$A_{\rho }=0.25$
			$B_{\rho }=10$
		$\varpi \sim \text{Beta}(A_{\varpi },B_{\varpi })$	$A_{\varpi }=3$
			$B_{\varpi }=6$

Table 2. Table 2 : Posterior summary of effects for each of the three outcomes

Effect	ED/IP visit		Rescue Inhaler Usage		OCS Usage
	Estimate	95% CI	Estimate	95% CI	Estimate	95% CI
step_level	$-0.264$	$(-0.388,-0.183)$	$-0.017$	$(-0.023,-0.013)$	$-0.189$	$(-0.228,-0.161)$
sex	$0.270$	$(-0.307,0.917)$	$0.137$	$(-0.216,0.426)$	$0.082$	$(-0.350,0.681)$
age	$-0.028$	$(-0.046,-0.013)$	$-0.009$	$(-0.016,0.001)$	$0.033$	$(0.021,0.042)$
msa	$-0.294$	$(-0.941,0.384)$	$-0.325$	$(-0.692,0.032)$	$0.318$	$(-0.154,0.956)$
smoke1	$-0.244$	$(-1.035,0.529)$	$0.353$	$(-0.124,0.697)$	$-0.029$	$(-0.605,0.525)$
smoke2	$-0.033$	$(-0.794,0.594)$	$0.033$	$(-0.425,0.408)$	$-0.106$	$(-0.917,0.415)$
ch_count	$-0.412$	$(-0.941,0.065)$	$0.116$	$(-0.113,0.323)$	$0.008$	$(-0.279,0.389)$
dep_flag2	$0.400$	$(-0.465,1.285)$	$0.108$	$(-0.400,0.648)$	$-0.615$	$(-1.318,0.067)$
copd_flag2	$0.298$	$(-0.318,0.863)$	$0.467$	$(0.127,0.931)$	$-0.295$	$(-0.759,0.141)$
gerd_flag2	$-0.204$	$(-0.946,0.444)$	$0.017$	$(-0.387,0.426)$	$0.226$	$(-0.290,0.651)$
resp_flag2	$-0.069$	$(-0.647,0.449)$	$0.051$	$(-0.321,0.432)$	$-0.337$	$(-0.818,0.275)$

Table 3. Table 3 : Summary data for patients 1 and 2 over the five rounds of the study.

		Round
		1	2	3	4	5
Patient 1	Step Level	3	2	4	5	0
	ED/IP Visits	0	0	0	0	2
	Rescue Inhaler	1	0	0	0	0
	OCS Use	0	1	0	0	0
Patient 2	Step Level	3	3	0	3	3
	ED/IP Visits	0	0	0	0	0
	Rescue Inhaler	0	0	0	0	0
	OCS Use	0	0	0	0	0

Table 4. Table 4 : Posterior probability of adverse event in next 6 months for patients 1 and 2 at each of the possible medication step levels. The medication level the respective patient was at in the last round of the study is in bold.

	Step Level	Probability of
		Adverse Event
Patient 1	0	0.453
	1	0.369
	2	0.291
	3	0.226
	4	0.175
	5	0.131
Patient 2	0	0.142
	1	0.111
	2	0.086
	3	0.067
	4	0.051
	5	0.040

Table 5. Table 5 : R 2 superscript 𝑅 2 R^{2} of the predicted probability of an event in the next six months to the true probability of event; averaged over the 100 simulations along with (standard error) for each outcome: (i) any ED/IP visit, (ii) four or more rescue inhaler refills, (iii) any oral steroid use, (iv) exacerbation, i.e., any of (i), (ii), or (iii).

Method	Outcome
	ED/IP	OCS	$\ge 4$ rescue	Any
Latent Process	0.713 (0.005)	0.836 (0.002)	0.765 (0.004)	0.811 (0.002)
Latent Indep	0.607 (0.005)	0.818 (0.003)	0.728 (0.004)	0.779 (0.002)
GLMER	0.454 (0.008)	0.696 (0.006)	0.657 (0.007)	0.698 (0.005)

Equations

$$y_{i,n,1}$$

$$y_{i,n,2}$$

$$y_{i,n,3}$$

$$y_{i,m,2}^{\ast }$$

$$y_{i,m,3}^{\ast }$$

$$y_{i,n,l}\sim \text{mbox}Poisson({\lambda }_{i,n,l}),$$

$${\lambda }_{i,n,l}={\int }_{t_{n-1}}^{t_n}{\eta }_{i,l}(s)ds,\text{vspace}-.0in$$

$$\log \left[{\eta }_{i,l}(s)\right]={\alpha }_{i,l}+\sum _{j=1}^J{\beta }_{j,l}{x}_{i,j}+{\gamma }_l{z}_i(s)+{\delta }_{i,l}(s),\text{vspace}-.05in$$

$$y_{i,m,2}^{\ast }\sim \text{mbox}Poisson\left(\sum _{n\in R_{i,m}}{\lambda }_{i,n,2}\right),$$

$$y_{i,m,3}^{\ast }\sim \text{mbox}Poisson\left(\sum _{n\in R_{i,m}}{\lambda }_{i,n,3}\right).$$

$$\mbox{\boldmath$\alpha$}_{i}=\left[\alpha_{i,1},\dots,\alpha_{i,L}\right]^{\prime}\sim N\left(\mbox{\boldmath$\nu$},\mbox{\boldmath$\Sigma$}\right).\vspace{-.1in}$$

$$K_l(s,t)={\tau }_l\exp \{-\theta |s-t|\}.\text{vspace}-.1in$$

$$\mbox{Cov}\left(\left[\delta_{i,1}(s),\dots,\delta_{i,L}(s)\right]^{\prime}\right)=\mbox{\boldmath$\Phi$},\vspace{-.1in}$$

$$K\left((l,s),(l^{\prime},t)\right)=\mbox{Cov}\left(\delta_{i,l}(s),\delta_{i,l^{\prime}}(t)\right)=\mbox{\boldmath$\Phi$}_{l,l^{\prime}}\exp\{-\theta|s-t|\}.\vspace{-.1in}$$

$$z_i(t)=\sum _{k=0}^{\infty }[{\mu }_h+{\tilde{z}}_{i,h}(t)]I_{\{B_i(t)=k\}}{I}_{\{{\mu }_h+{\tilde{z}}_{i,h}(t)>0\}},\text{vspace}-.05in$$

$$B_i(t)=C_i(t)+D_i(t),\text{vspace}-.15in$$

$$\begin{array}[]{rcl}\Pr\{D_{i}(t^{+})=D_{i}(t)+1\}&=&\varpi,\\[7.22743pt] \Pr\{D_{i}(t^{+})=D_{i}(t)\}&=&1-\varpi.\end{array}\vspace{-.0in}$$

$$z_{i,m}^{\ast }=\min \left\{5,\lfloor {\int }_{t\in R_{i,m}}{z}_i(t)dt+0.5\rfloor \right\}.$$

$$\Theta=\left\{\{\mbox{\boldmath$\alpha$}_{i}\}_{i=1}^{n},\{\mbox{\boldmath$\delta$}_{i}\}_{i=1}^{n},\{\tilde{z}_{i,1},\dots,\tilde{z}_{i,H}\}_{i=1}^{n},\{B_{i}\}_{i=1}^{n},\{\mbox{\boldmath$\beta$}_{j}\}_{j=1}^{J},\mbox{\boldmath$\gamma$},\mbox{\boldmath$\Sigma$},\theta,\mbox{\boldmath$\Psi$},\mu,\sigma^{2},\phi,\varpi\right\},$$

Full text

1 Introduction

Asthma is a chronic disease, affecting 8% of the population and costing $37-63 billion/year in the US (Kamble and Bharmal, 2009; Jang, Chan, Huang and Sullivan, 2013). Asthma can lead to death and also includes significant morbidity including hospitalization, emergency department visits, outpatient visits, missed school and work days, decreased productivity, and decreased physical activity. Current strategies for managing asthma include avoidance of triggers, anticipatory guidance for exacerbation treatment, and the scheduled use of various asthma medications as “preventive” medications. Medications account for the highest proportion of healthcare spending for asthma. In addition, spending on medications has been growing at a rapid rate over the last decade (Rank, Liesinger, Ziegenfuss, Branda, Lim, Yawn, Li and Shah, 2012). While medication spending continues to increase, this has not resulted in improved asthma outcomes (Moorman, Akinbami, Bailey, Zahran, King, Johnson and Liu, 2012; Rank, Liesinger, Ziegenfuss, Branda, Lim, Yawn and Shah, 2012). To enhance the management of asthma, we propose to develop a risk prediction model to assist clinicians and patients when deciding how best to manage asthma medication for a given patient.

The current strategy of asthma medication management includes an assessment of current symptom control and future risk for asthma exacerbations which is then matched for treatment intensity. Assessing current symptoms can be accomplished by using standardized patient-reported tools; however, estimating future risk is more challenging. Prior efforts to predict future exacerbation events focused on specific risk factors available in cross-section (e.g. lung function, blood eosinophil levels, prior exacerbations, or a combination of these and other factors) (Bateman, Buhl, O’Byrne, Humbert, Reddel, Sears, Jenkins, Harrison, Quirce, Peterson et al., 2015; Blakey, Woulnough, James, Fellows, Obeidat, Navaratnam, Stringfellow, Yeoh, Pavord, Thomas et al., 2012). However, each individual with asthma may have a differential response to medications over time. Using this individual response data to predict outcomes at differing medication levels has not been previously attempted. Our hypothesis is that using response to medication data at the individual level will provide a meaningful future risk estimate. Individualized estimation of future risk will provide a key step in personalizing asthma medication regimens.

A total of $4,235$ patients with persistent asthma were identified from the US Medical Expenditure Panel Survey (MEPS) years 2000-2010. Each patient had data for $\sim$2 years, and measurement was divided into 5 rounds (periods of 4 to 5 months each). Medication level, taking on ordinal levels of 0 $-$ 5, was defined by type and dose of chronic asthma medication based on current guidelines as described in Rank, Liesinger, Branda, Gionfriddo, Schatz, Zeiger and Shah (2016). This definition of medication level only includes controller medication and does not include rescue inhaler or oral steroid usage which are both considered outcomes here. The primary outcomes of interest are (i) emergency department (ED) and/or in-patient (IP) visits, (ii) Rescue inhaler refills, and (iii) Oral steroid use. The goal of this work is to provide an individualized risk prediction of an adverse event in the next six months for a given patient on a given medication level. An adverse event is defined here to be any of the following: (i) any ED/IP visit, (ii) any oral steroid use, (iii) four or more rescue inhaler refills.

Unfortunately in the MEPS data, medication level is observed only as an aggregate (i.e., the average level) for an entire round. This is problematic because in reality medications are often introduced or increased at an ED/IP visit (which often occur in the middle of a round). If this were to be ignored and medication level were assumed constant throughout a round, then the MEPS data would actually imply that higher medication levels for a round are associated with a higher chance of an ED/IP visit. Thus, some care is needed when modeling the effect of medication level. We propose a latent random process for the medication level for an individual so that it is allowed to change within a round, with a positive probability of changing at the time of an IP, ED, office-based (OB), or out-patient (OP) visit. This process is restricted so that the average medication level for an entire round is consistent with the observed, average medication level for that round in the MEPS data. Latent variable approaches are common in the literature and have been used in similar situations (e.g., Muthen, 1983; Dunson, 2000; Storlie, Michalak, Quinn, Dubois, Wender and Dubois, 2013). Conditional on the latent medication process, the model for the multivariate response over time is a dependent log-Gaussian Cox Process (Møller, Syversveen and Waagepetersen, 1998; Brix and Diggle, 2001) with a mean function that depends on a patient’s covariates and a random intercept. Thus, this model can be considered as a measurement error model (Carroll, Ruppert, Stefanski and Crainiceanu, 2006) on a time-varying covariate (Wulfsohn and Tsiatis, 1997; Tsiatis and Davidian, 2001). The unique features here are that the time varying covariate cannot be assumed to be a smooth process over time and the variable is never measured/known exactly at any particular times. In this case, only the average level is known over several partitions of time and it is quite possible for the process to jump from one level to another. Estimation is performed in a Bayesian framework via Markov Chain Monte Carlo (MCMC) in a similar manner to that for handling multivariate missing data (Storlie, Therneau, Carter, Chia, Bergquist and Romero-Brufau, 2017, 2017).

The remainder of the paper is laid out as follows. Section 2 describes the statistical approach used to model patient outcomes, while computational details are discussed in Section 3. Analysis of the MEPS data and a summary of predictive accuracy is provided in Section 4. Section 5 summarizes the results of applying the model to a simulated test case designed to be similar to the the MEPS data and Section 6 concludes the paper.

2 Statistical Model

The outcomes of rescue inhaler usage and oral steroid usage are aggregated at the round level. However, ED/IP visits are captured at the day granularity, so the model below is defined on a finer granularity of time period within each round. In principle the time period could be as fine as one day, however, one week granularity was used in all results below for computational expediency. The response for the $i^{\mbox{\scriptsize th}}$ individual for period $n$ (i.e., week $n$) is \mbox{\boldmathy}_{i,n}=[y_{i,n,1},y_{i,n,2},y_{i,n,3}]^{\prime}, with

[TABLE]

As discussed above, the $y_{i,n,1}$ are actually observed, but only the sum of $y_{i,n,2}$ and $y_{i,n,3}$ for the $m^{\mbox{\scriptsize th}}$ round are observed, i.e.,

[TABLE]

where $R_{i,m}$ is the set of time periods that make up the $m^{\mbox{\scriptsize th}}$ round.

The distribution of these responses was assumed to be a function of patient dependent covariates \mbox{\boldmathx}_{i}=[x_{i,1},\dots,x_{i,p}]^{\prime} (including physical characteristics, co-morbidities, etc.) and their asthma controller medication level $z_{i}(t)$ at time $t$. The controller medication level takes on ordinal levels of 0 $-$ 5 as defined in Rank, Liesinger, Branda, Gionfriddo, Schatz, Zeiger and Shah (2016). In the MEPS data, $z_{i}(t)$ is observed only as an aggregate (i.e., the average level) for an entire round. Hence, the $z_{i}(t)$ were assumed here to be (latent) random processes so that they can be allowed to change within a round. It is essential to allow for this in order to remove the bias introduced by the fact that medications are often introduced and/or increased at an ED/IP visit. If this were to be ignored in a model, then it would incorrectly imply that higher medication levels increase the chance of an ED/IP visit. Thus, $z_{i}$ was allowed to be a Markov process with the restriction that the average medication level for an entire round must be consistent with the observed (aggregated) medication level for that round in the MEPS data. There is a baseline rate at which increases or decreases to medication level may occur at any given time, but there is also a probability $\varpi$ of changing immediately after a IP/OP/OB/ED visit (i.e., the probability of a change in a given week is allowed higher if there was contact with a provider during that week). The specifics of the model used for controller medication level are described later in this section. Identifiability concerns are addressed in Section 3.

The distribution of \mbox{\boldmathy}_{i,n} also depends on a time varying random effect for the $i^{\mbox{\scriptsize th}}$ individual to account for the overall severity of their condition compared to others in the population and the variation in the severity of an individual’s symptoms over time. A model to account for all of this is the following hierarchical Poisson process,

[TABLE]

where

[TABLE]

where $t_{n}$ is the end of the $n^{\mbox{\scriptsize th}}$ time period and

[TABLE]

where $l=1,2,3$ is indexing the response type (i.e., ED/IP visits, rescue inhaler usage, and oral steroid usage).

• •

$\alpha_{i,l}$ is the random effect on the $l^{\mbox{\scriptsize th}}$ response type, to account for the overall level of severity for the $i^{\mbox{\scriptsize th}}$ individual and for the individual’s response to medication.

• •

$\beta_{j,l}$ are fixed effects, constant across the population.

• •

$\gamma_{l}$ is the fixed effect for controller medication level on the $l^{\mbox{\scriptsize th}}$ response type.

• •

$\delta_{i,l}(s)$ is a time varying random effect to account for the natural fluctuation/trends in an individual’s condition.

This model implies that the round aggregated responses follow,

[TABLE]

In order to complete the model specification a parent distribution must be defined for the random effects $\alpha_{i,l}$ and $\delta_{i,l}$. A multivariate normal distribution was assumed for the $\alpha_{i,l}$,

[TABLE]

The time varying effect $\delta_{i,l}$ was assumed to be a Gaussian process with mean 0 and covariance function $K_{l}(s,t)$; specifically, it is a stationary Ornstein-Uhlenbeck process (equivalent to an order one, autoregressive time series) for each $l$, i.e.,

[TABLE]

A product correlation was also assumed between $\delta_{i,l}$ and $\delta_{i,l^{\prime}}$ so that the variation in the rate of occurrence for the three response types are correlated, i.e.,

[TABLE]

or the general covariance function is then

[TABLE]

With the above model, for a given individual and their response history, predictions can be made for risk of future outcomes at various levels of medication. This prediction is a (correlated) distribution for the three responses so medication decisions can be assessed using the probability that the patient would need fewer than 4 rescue inhalers, not have an ED/IP visit, and not need OCS in the next six months, for example. This probability or some similar measure can be used to assess if this individual is a good candidate for a step down in medication.

Finally, a model for the controller medication level, $z_{i}(t)$, is specified as

[TABLE]

where (i) $\tilde{z}_{i,h}$, $h=0,1,\dots$, are iid mean zero Ornstein-Uhlenbeck (O-U) processes with variance $\sigma$ and correlation parameter $\phi$, (ii) $I_{A}=1$ if $A$, and 0 otherwise (e.g., the last term in (6) ensures that $z_{i}(t)\geq 0$), (iii) $\mu_{h}\sim N(M_{\mu},S^{2}_{\mu})$, and (iv) the $B_{i}(t)$ are independent counting process intended to allow medication level to make large (discontinuous) changes at certain points in time. That is,

[TABLE]

where the $C_{i}$ allow the possibility of abrupt changes in medication level at any time, while the $D_{i}$ allow for a positive (increased) probability of an abrupt change at the time of a provider contact. Specifically, it was assumed that $C_{i}$ is a Poisson process with intensity $\rho$ and $D_{i}(t)$ is a counting process that is allowed to increment by one, only at the time of an IP/OP/OB/ED visit. That is, $D_{i}(0)=0$ and $D_{i}(t^{+})=\lim_{s\downarrow t}D_{i}(s)=D_{i}(t)$ for all $t$, unless $t$ is the time of an IP/OP/OB/ED visit, in which case

[TABLE]

As mentioned above, the aggregated medication level for a round provided by the MEPS data $z^{*}_{i,m}$ was assumed to be the average of $z_{i}(t)$ over the $m^{\mbox{\scriptsize th}}$ round. In the above model, the latent variable $z_{i}$ is continuous even though the observed value for the $m^{\mbox{\scriptsize th}}$ round, $z^{*}_{i,m}$ is ordinal on $0,\dots,5$. This implies that the underlying medication level is continuous, but it is observed as an integer as a result of interval censoring on the cutpoints 0.5, 1.5, 2.5, 3.5, 4.5. For example, an average medication level over round $m$ may be 2.3 or 6.2, which would be “observed” as a medication levels of 2 and 5, respectively. Thus, the relation from $z_{i}(t)$ to $z^{*}_{i,m}$ is

[TABLE]

Several realizations of $z_{i}$ resulting from posterior distribution for an example patient are provided in Figure 1 along with the corresponding $z^{*}_{i,m}$ for reference. The process described above for $z_{i}$ is Markov, which is appropriate for this application. However, there may be situations where the medication level has a strong dependency on the entire history of previous medications. Intuitively this should not have a large influence on the main goal is interpolation of medication level as opposed to prediction into the future. However, there may be situations where this dependence to prior history is very strong and could be leveraged into a better model, perhaps via other covariance functions (e.g., powered exponential or Matern) and/or non-Markovian counting processes, etc.

3 Prior Distributions and Computation

A summary of the model structure and the assumed prior distributions are provided in Table 1. Relatively diffuse priors were used for most of the parameters, with certain exceptions. For example, the prior for $\gamma_{l}$ ensures that the effect of medication level is protective. Also, the parameters governing the medication level process $z_{i}$ were chosen to convey some information, guided by input from clinicians. The prior for $\sigma$ and $\phi$ assumes a relatively consistent medication level in the absence of a true (discontinuous) change; allowing for the possibility of a change in medication level of approximately one step in six months. Meanwhile $\rho$ was chosen to encourage a probability of anywhere from 0 to about a 25% chance of suddenly changing medication in a given week. Finally, the prior for $\varpi$ was chosen to reflect the perception that a patient is anywhere from 5% to 75%, with a best guess of about 33%, likely to make a change in medication shortly after a ED/IP/OP/OB visit.

Since fairly weak priors distributions were used for most parameters, there should be very little prior sensitivity with respect to such parameters. However, parameters governing the medication level process $z_{i}$ had to be chosen more carefully since there is some confounding between the parameters $(\sigma,\phi)$ and $(\rho,\varpi)$. That is, if the medication level needs to change a lot in a given window of time, it can do it either via a big variance for the continuous process or via many discrete jumps. Thus, it is important to assess the extent of the sensitivity of the results to the choices for these priors. Therefore, the MEPS analysis in Section 4 was conducted with prior parameter values specified in Table 1, and also with values that conveyed less prior information. Namely, values of $A_{\sigma}=2$, $A_{\sigma}=.02$, $A_{\phi}=2$, $B_{\phi}=1$, $A_{\rho}=0.1$, $B_{\rho}=4$, $A_{\varpi}=3$, and $B_{\varpi}=6$ were also used. The parameter estimates (and thus also the predictions) did not change in any significant manner. Finally, very diffuse values of $A_{\sigma}=.1$, $A_{\sigma}=.02$, $A_{\phi}=1$, $B_{\phi}=1$, $A_{\rho}=0.1$, $B_{\rho}=1$, $A_{\varpi}=1$, and $B_{\varpi}=1$ were then used. In this case, the estimates were qualitatively similar, but a bit different in magnitude (smaller) for the medication effect. This is because with these very diffuse prior parameter settings, it favors a latent medication process that makes changes continuously as opposes to abruptly (i.e., unlike that in Figure 1. This is a bit in contrast to our knowledge and intuition of medication patterns. Thus, some care is required when specifying prior distributions for this model. Since it can be challenging to translate knowledge about the medication process into the priors directly, a useful sanity check is to do the following. Generate several curves from the proposed prior and make sure that the resulting latent medication process realizations that the prior can produce are consistent with what is known about the process.

The posterior distribution was approximated via MCMC. The complete list of parameters sampled in the MCMC algorithm is

[TABLE]

where \mbox{\boldmath\beta}_{j}=[\beta_{j,1},\beta_{j,2},\beta_{j,3}]^{\prime}, \mbox{\boldmath\gamma}=[\gamma_{1},\gamma_{2},\gamma_{3}]^{\prime}, and the sequence of $\tilde{z}_{i,h}$ is truncated at $h=H$. The upper bound $H$ can be chosen such that it is large enough not to affect the results. A value of $H=10$ was used to generate the results below and was deemed sufficient since all of the posterior samples had $\max_{i,t}B_{i}(t)<10$. As was done in Storlie, Michalak, Quinn, Dubois, Wender and Dubois (2013), the continuous functions $B_{i}$, $\tilde{z}_{i,h}$ and $\delta_{i,l}$ are resolved on a fine time grid, i.e., $B_{i}\equiv[B_{i}(t_{i,1}),\dots,B_{i}(t_{i,N_{i}})]^{\prime}$ and similarly for the $\tilde{z}_{i,h}$ and $\delta_{i,l}$. The corresponding integrals involving these functions needed to evaluate the likelihood are then approximated with quadrature. One week granularity was used for the $t_{i,n}$ to match that used for the response in (1).

The MCMC was carried out within a Gibbs framework via the rjags package in R (Plummer, 2009), where each of the elements of $\Theta$ are updated in turn. Five parallel chains were run for 20,000 iterations, the first 10,000 of which were discarded as burn-in. All chains converged to the same distribution, as gauged by the non-individual specific (not $i$ dependent) parameters, and were thus aggregated. On the MEPS dataset with $n=4,235$ individuals ($\sim 300,000$ total time periods), these 20,000 iterations took $\sim 40$ hours. The Jags model R code is provided in a GitHub repository at https://github.com/cbstorlie/asthma_stepdown. Jags is a very convenient, general framework with which to conduct Gibbs sampling for Bayesian models, however, this convenience can come at the cost of computation time. In this hierarchical model, for each MCMC iteration, the individual specific parameters \{\mbox{\boldmath\alpha}_{i},\mbox{\boldmath\delta}_{i},\tilde{z}_{i,1},\dots,\tilde{z}_{i,H},B_{i}\} could, in principle, be updated for each $i$ independently and thus in parallel. This would substantially decrease computation time, however, this is not possible when using Jags and would require custom MCMC code. Since training is done off-line and individual predictions are quite fast (on the order of $\sim$1 second), Jags seems sufficient for this application.

4 MEPS Analysis Results

4.1 Inference

The patient specific variables considered in addition to medication level were sex, age, residing in a metropolitan statistical area (msa), smoking status (smoke1 = 1 for smoker, smoke2 = 1 for non-smoker but with smoking exposure), Charlson sum of diseases (ch_count), depression/anxiety (dep), chronic obstructive pulmonary disease (copd), gastroesophageal reflux disease (gerd), and rhinitis/sinusitis (resp). These variables were chosen as they were the only variables readily available from MEPS that were thought to potentially have an effect on the outcomes. Other variables such as eosinophil levels, and pulmonary function, were not available. The posterior distribution for the fixed effects is summarized in Table 2. Increased controller medication level is strongly related to a lower rate of ED/IP visits and OCS usage. However, medication level appears to have a much smaller influence on rescue inhaler usage. Age is the only other factor among those investigated that has a significant effect. According to the model, older patients are less likely to have an ED/IP visit, while they are more likely to use OCS.

4.2 Prediction of an Adverse Event

Let an adverse event for a six month period be defined as any ED/IP visit, any OCS use, or more than 3 rescue refills. For a given draw of the parameters from the MCMC sample, the probability of an adverse event for a given patient at a given medication level can be obtained analytically based on the Poisson CDF and the conditional independence of the three outcomes. The posterior predictive probability of an adverse event can then be obtained as the posterior mean of these probabilities. Example patient 1 (14${}^{\mbox{\scriptsize th}}$ individual in the data set), was at medication levels of 3, 2, 4, 5, 0 in his/her five rounds, respectively (see Table 3). This patient had an adverse event in rounds 2 and 5 according to the definition above. The posterior probability of an adverse event in the next 6 months if he/she stays at medication step level 0 is 0.453 (Table 4). It may seem more appropriate for this patient to be on an at least a step level of 3 in order to lower this probability to a more reasonable $\sim$20% level.

Conversely, patient 2 (24${}^{\mbox{\scriptsize th}}$ individual in the data set) was at step level 3 in the last round. The probability of an adverse event at each of the medication step levels for the next six month period is also provided in Table 4. If this patient remains at step level 3, the probability of an adverse event is quite low at 0.067. The provider and patient may decide that it is appropriate to reduce the patient’s medication level in this case.

4.3 Validation

Validating the predictive accuracy of the model using MEPS data is difficult due to the fact that medication level over time is not entirely known within a round. However, just for this validation exercise, it was assumed that medication level remained constant during the last round (round 5) for a set of 500 patients. All observations from round 5 for 500 randomly selected patients were held out of the model fitting process and the model was then used to predict the probability of an adverse event during round 5 for each of these patients. This validation exercise is conservative since the model was given potentially false information for medication level (assumed constant) as input, yet it has to predict the reality (produced by medication level changing during a round). Even still, comparison of such predictions top the observed outcomes should provide a useful sanity check.

The histogram of the distribution of these predicted probabilities is provided in Figure 2(a). The majority of predicted probabilities are below 0.25, but there is a relatively uniform spread of probabilities from 0.25 to 1.0. Figure 2(b) provides a graphical comparison of the 500 predictions to the observed outcomes similar to that used in Harrell, Lee and Mark (1996) and Austin and Steyerberg (2014). Since the outcomes for each individual are 0 or 1, they are aggregated into various bins of predicted risk (in this case the 10 bins resulting from the cut points $0.0,0.1,\dots,0.9,1.0$). The ten plotted points are the average of the predicted probability within the respective bin against the empirical (or observed) probability of an event within that bin. If the model is reasonable, it would be expected that the blue points would fall close to the red $y=x$ line. Individual 95% confidence regions for the predicted versus observed points in each bin were obtained via a nonparametric bootstrap approach. It would be expected that each of these intervals would intersect the red line 95% of the time if the model was correct.

Intuition would suggest that the predictions provided by the model are almost certainly incorrect in this case, since the medication level being provided for prediction is likely not constant during the final round for many of the 500 patients. Nonetheless, the model predictions appear to agree very well with the observations for smaller values of predicted risk. There may be some upward bias in predictions of higher risk individuals, i.e., the model predictions are a bit conservative. However, it is far less critical to have predictive probabilities above 0.5 align well with reality (i.e., anything approaching 0.5 is “high”, and it is not important to distinguish a 0.6 from a 0.7, for example). It is precisely the smaller risk predictions that need to be accurate for the clinical use case. Fortunately, there is no evidence of prediction bias among the lower risk ($\leq 0.5$) predictions.

5 Simulation Study

In this section, a further examination of the model performance is provided on a simulation case designed to be similar to the MEPS data. Specifically, the posterior mean from the MEPS analysis for all parameters was used to generate a reduced data set for the first $n=500$ individuals. The medication level for each patient was set to the average step level he/she had in round 1 in the MEPS data. In order to mimic a key feature of the MEPS data, upon an ED/IP visit (as randomly generated by the model) it was assumed that, with probability 0.5, the medication step level would be increased to a randomly selected higher level (as opposed to a random change in either direction). The same mechanism is applied to an OCS prescription fill as well. There is also 0.3% chance (the posterior mean estimate) that the medication level will change for a simulated patient on any given day (either up or down) in a random fashion with probability of the new level proportional to the frequency that step level was observed in the MEPS data. This process was repeated to create a total of 100 synthetic data sets.

The proposed model (Latent Process) was then fit to each of the 100 data sets and the resulting posterior means of the parameters in Table 2 are summarized as boxplots in Figures 3(a) - 3(c). The “observed” data used to fit the model was kept consistent with the MEPS data. That is, it was assumed that the average medication level for a round was the only information about controller medication available, ED/IP visits were known to the day, but rescue inhaler fills and OCS were only known in aggregate for an entire round. The estimates provided in Figures 3 are the respective posterior means for each of the 100 data sets. The coefficients have been standardized by the standard deviation of the respective variables to make them comparable. A variant of the proposed approach was also used that assumed independent outcomes (Latent Indep), i.e., restricted such that $\Sigma$ and $\Psi$ are diagonal with independent inverse-gamma priors on the diagonal elements. This approach was included to assess how much impact the treatment of multivariate dependence may have on the modeling results. A generalized linear mixed effects regression (GLMER) using the R package lme4 was also fit separately to each of the three responses. The GLMER model has a random intercept for each patient like the proposed model, however, the medication level was assumed fixed at the average value observed for the round. This was intended to provide an evaluation of what might be gained by treating the medication level as a latent process. Boxplots of the corresponding GLMER estimates are provided in Figures 3(a) - 3(c) as well. The true value of each parameter used to generate the 100 datasets is provided as a red line in each boxplot for comparison.

It is clear from Figure 3(a) that GLMER substantially underestimates the magnitude of the effect that medication has on ED/IP visits. The proposed approach on the other hand, has a slight bias for this coefficient (as to be expected with a Bayesian approach) but it provides a far better estimate in general since it allows for the possibility of changes in medication to occur at/after an ED/IP visit. The proposed latent medication approach with independent outcomes has less bias than GLMER, but it is not as close to the true value as the full multivariate approach. The GLMER estimate is also biased toward zero on the effect that medication has on OCS use. The practical ramifications of the substantial bias introduced by GLMER for this problem can be seen in Figure 4.

Figure 4(a) is a calibration plot similar to that in Figure 2(b), but with effectively no error in the estimated points due to the 100 repeated simulations providing enough data to make the area of the confidence regions negligible. Also, the medication level for the patients to be predicted in Figure 4(a) was randomly chosen to be a step down from the patients’ respective medication level at the end of the simulated training data. Restricting to a step down was done to highlight the issues that are present with GLMER predictions on the most important use case. If a simulated patient was at step 0 at the end of the simulated training data, then they could not be stepped down, and were thus not included in this particular prediction exercise (usually about 150/500 were at step level 0 and excluded). Figure 4(a) shows that the predictions from the proposed approach are well calibrated, but that the GLMER predictions are substantially optimistic; an average probability of 0.07 according to GLMER leads to an actual probability of adverse event of over 0.20. This would cause recommendations for medication step down to individuals that should not be lowering their medication level more often than is desirable.

Figure 4(b) provides ROC curves for predicted probability of an event for the three approaches, using a randomly selected medication level for the next six months for each hypothetical subject. The figure demonstrates better overall predictive capability of the proposed approach, above and beyond better calibration for the case of a step down. Table 5 provides the $R^{2}$ measures of the predicted probability of an event in the next six months to the true probability of event for each outcome separately, and for the aggregated adverse event outcome (any). This table demonstrates that the estimated probabilities obtained from the proposed latent process approach are far more accurate than those from GLMER which assumes that medication level is constant between rounds in addition to independent outcomes. The Latent Process approach which is the full multivariate version of the proposed approach is also significantly more accurate than Latent Indep which still uses the latent process for medication, but assumes that the outcomes are dependent. Interestingly, the Latent Process approach performs better prediction of the individual outcomes than Latent Indep due to the ability to borrow strength from each outcome when estimating an individual’s overall and current level of risk for a particular outcome.

6 Conclusions & Further Work

In this paper we developed a general approach to prediction of individual outcomes for patients using prescription data that is ambiguous about the exact medication level a patient is at over time. The statistical approach treated medication level as an uncertain latent process in a measurement error framework, to reduce the bias in outcome prediction. This approach was applied to a population of asthma sufferers using the data available from the MEPS database years 2000-2010. The analysis demonstrated that well-calibrated predictions can be made for the individual-specific probability of an adverse outcome in a future time period. While there is some upward bias of predicted risk for larger risk individuals, the smaller risk predictions agreed very well with the validation set observations. The approach was also compared via a simulation study to a more standard GLMER model with a random intercept, ignoring the uncertainty in medication level. These results showed the potential for GLMER to be significantly biased when predicting outcomes due to the incorrect assumption of constant medication level when training the model.

It would be prudent to validate the model prediction using a separate data set where the prescription refill dates are more precisely known. Such data would also allow for improved model estimation, though the medication level would still need a latent process treatment due to the fact that refill dates do not provide information about daily use, when medication is discontinued, etc. The goal is to ultimately perform pilot trials of embedding the clinical prediction model within the healthcare delivery system using (i) shared-decision making tools used within the patient-provider visit, (ii) electronic medical record alerts directed to providers, and (iii) risk information shared directly with patients using secure messaging.