# Prediction of individual outcomes for asthma sufferers

Prediction of individual outcomes for asthma sufferers Summary 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 United States of America. 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. 1. Introduction Asthma is a chronic disease, affecting 8% of the population and costing$37–63 billion/year in the United States of America (Kamble and Bharmal, 2009; Jang and others, 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 health care spending for asthma. In addition, spending on medications has been growing at a rapid rate over the last decade (Rank and others, 2012a). While medication spending continues to increase, this has not resulted in improved asthma outcomes (Moorman and others, 2012; Rank and others, 2012b). 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 and others, 2015; Blakey and others, 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 $$4235$$ 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 five rounds (periods of 4–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 and others (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 6 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, and (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 and others, 2013). Conditional on the latent medication process, the model for the multivariate response over time is a dependent log-Gaussian Cox Process (Møller and others, 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 and others, 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). The remainder of the article 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, 1 week granularity was used in all results below for computational expediency. The response for the $$i$$th individual for period $$n$$ (i.e., week $$n$$) is $${\boldsymbol{y}}_{i,n} = [y_{i,n,1}, y_{i,n,2}, y_{i,n,3}]'$$, with \begin{eqnarray} y_{i,n,1} & = & \mbox{# of ED/IP visits for ith individual during nth time period.} \nonumber\\ y_{i,n,2} & = & \mbox{# of rescue inhaler refills for ith individual during nth time period.} \\ y_{i,n,3} & = & \mbox{# of oral steroid refills for ith individual during nth time period.} \nonumber \end{eqnarray} (2.1) 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$$th round are observed, i.e., \begin{eqnarray*} y^*_{i,m,2} &=& \sum_{n \in R_{i,m}} y_{i,n,2}, \nonumber\\ y^*_{i,m,3} &=& \sum_{n \in R_{i,m}} y_{i,n,3}, \end{eqnarray*} where $$R_{i,m}$$ is the set of time periods that make up the $$m$$th round. The distribution of these responses was assumed to be a function of patient dependent covariates $${\boldsymbol{x}}_{i} = [x_{i,1},\dots,x_{i,p}]'$$ (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 and others (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 $${\boldsymbol{y}}_{i,n}$$ also depends on a time varying random effect for the $$i$$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, $$y_{i,n,l} \sim \mbox{Poisson}(\lambda_{i,n,l}),$$ (2.2) where \begin{equation*} \lambda_{i,n,l} = \int_{t_{n-1}}^{t_n} \eta_{i,l}(s) {\rm d}s, \end{equation*} where $$t_n$$ is the end of the $$n$$th time period and $$\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),$$ (2.3) 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$$th response type, to account for the overall level of severity for the $$i$$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$$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, \begin{eqnarray*} y^*_{i,m,2} \sim \mbox{Poisson}\left(\sum_{n \in R_{i,m}} \! \lambda_{i,n,2}\right)\!, \nonumber \\ y^*_{i,m,3} \sim \mbox{Poisson}\left(\sum_{n \in R_{i,m}} \! \lambda_{i,n,3}\right)\!. \end{eqnarray*} 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}$$, $${\boldsymbol{\alpha}}_i = \left[\alpha_{i,1},\dots,\alpha_{i,L} \right]' \sim N \left( {\boldsymbol{\nu}}, {\boldsymbol{\Sigma}} \right)\!.$$ (2.4) 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 (O–U) process (equivalent to an order one, autoregressive time series) for each $$l$$, i.e., \begin{equation*} K_l(s,t) = \tau_l \exp \{-\theta |s-t|\}. \end{equation*} A product correlation was also assumed between $$\delta_{i,l}$$ and $$\delta_{i,l'}$$ so that the variation in the rate of occurrence for the three response types are correlated, i.e., \begin{equation*} {\text{Cov}} \left( \left[\delta_{i,1}(s), \dots, \delta_{i,L}(s) \right]' \right) = {\boldsymbol{\Phi}}, \end{equation*} or the general covariance function is then $$K\left((l,s),(l',t) \right) = {\text{Cov}} \left( \delta_{i,l}(s), \delta_{i,l'}(t) \right) = {\boldsymbol{\Phi}}_{l,l'} \exp \{-\theta |s-t|\}.$$ (2.5) 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 6 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 $$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\}},$$ (2.6) where (i) $$\tilde{z}_{i,h}$$, $$h=0,1,\dots$$, are iid mean zero 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 (2.6) ensures that $$z_i(t)\geq0$$), (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, \begin{equation*} B_i(t) = C_i(t) + D_i(t), \end{equation*} 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 \begin{equation*} \begin{array}{rcl} \Pr \{ D_i(t^+) = D_i(t)+1 \} & = & \varpi, \\[.1in] \Pr \{ D_i(t^+) = D_i(t) \} & = & 1- \varpi. \end{array} \end{equation*} 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$$th round. In the above model, the latent variable $$z_i$$ is continuous even though the observed value for the $$m$$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 cut-points 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 \begin{equation*} z^*_{i,m} = \min\left\{ 5\;,\; \left\lfloor \int_{t \in R_{i,m}} z_i(t){\rm d}t + 0.5 \right\rfloor \right\}\!. \end{equation*} 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 as 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. Fig. 1. View largeDownload slide Twenty-five posterior realizations of the latent controller medication process, $$z_i$$ for a given patient. When an increase in medication occurs, the model generally prefers to have it occur at ED/IP/OP/OB visits over other times. Fig. 1. View largeDownload slide Twenty-five posterior realizations of the latent controller medication process, $$z_i$$ for a given patient. When an increase in medication occurs, the model generally prefers to have it occur at ED/IP/OP/OB visits over other times. 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 6 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. Table 1. Summary of the model for asthma related events defined in (2.2)–(2.6) and the specification of prior distributions Description Model Prior distributions Specification Individual effects $${\boldsymbol{\alpha}}_{i} \stackrel{\rm iid}{\sim} N({\boldsymbol{\nu}}, {\boldsymbol{\Sigma}})$$, $$i=1,\ldots,n$$ $$\nu_l \stackrel{\rm ind}{\sim} N(M_{\nu,l}, S^2_{\nu,l})$$ $$M_{\nu,l}= 0$$ $$S^2_{\nu,l}= 100$$ $${\boldsymbol{\Sigma}} \sim \mbox{Wishart}({\boldsymbol{P}}_\Sigma, N_\Sigma)$$ $${\boldsymbol{P}}_\Sigma= {\boldsymbol{I}}$$ $$N_\Sigma=4$$ Fixed effects $$\beta_{j,l}$$ in (2.3) $$\beta_{j,l} \stackrel{\rm ind}{\sim} N(B_{j,l}, T^2_{j,l})$$ $$B_{j,l}= 0$$ $$T^2_{j,l}= 100$$ Medication effects $$\gamma_{l}$$ in (2.3) $$-\gamma_{l} \stackrel{\rm ind}{\sim} \log N(C_{l}, U^2_{l})$$ $$C_{l}= -1$$ $$U^2_{l}= 10$$ Individual time effect $${\boldsymbol{\delta}}_{i} {\sim} GP(0, K)$$, $$i=1,\ldots,n$$ $${\rm e}^{-\theta} \sim \mbox{Beta}(A_\theta, B_\theta)$$ $$A_\theta= 1$$ $$B_\theta= 1$$ $$K$$ in (2.5) $${\boldsymbol{\Phi}} \sim \mbox{Wishart}({\boldsymbol{P}}_\Phi, N_\Phi)$$ $${\boldsymbol{P}}_\Phi= {\boldsymbol{I}}$$ $$N_\Phi=4$$ $$\mu_h \sim N(M_\mu, S^2_\mu)$$ $$M_\mu=0$$ $$S^2_\mu= 50$$ Medication level process $$z_{i}$$ in (2.6) $$\sigma^2 \sim \mbox{Gamma}^{-1}(A_\sigma, B_\sigma)$$ $$A_\sigma= 5$$ $$B_\sigma=.05$$ $${\rm e}^{-\phi} \sim \mbox{Beta}(A_\phi, B_\phi)$$ $$A_\phi= 6$$ $$B_\phi= 3$$ $$\rho \sim \mbox{Gamma}(A_\rho, B_\rho)$$ $$A_\rho= 0.25$$ $$B_\rho= 10$$ $$\varpi \sim \mbox{Beta}(A_\varpi, B_\varpi)$$ $$A_\varpi= 3$$ $$B_\varpi= 6$$ Description Model Prior distributions Specification Individual effects $${\boldsymbol{\alpha}}_{i} \stackrel{\rm iid}{\sim} N({\boldsymbol{\nu}}, {\boldsymbol{\Sigma}})$$, $$i=1,\ldots,n$$ $$\nu_l \stackrel{\rm ind}{\sim} N(M_{\nu,l}, S^2_{\nu,l})$$ $$M_{\nu,l}= 0$$ $$S^2_{\nu,l}= 100$$ $${\boldsymbol{\Sigma}} \sim \mbox{Wishart}({\boldsymbol{P}}_\Sigma, N_\Sigma)$$ $${\boldsymbol{P}}_\Sigma= {\boldsymbol{I}}$$ $$N_\Sigma=4$$ Fixed effects $$\beta_{j,l}$$ in (2.3) $$\beta_{j,l} \stackrel{\rm ind}{\sim} N(B_{j,l}, T^2_{j,l})$$ $$B_{j,l}= 0$$ $$T^2_{j,l}= 100$$ Medication effects $$\gamma_{l}$$ in (2.3) $$-\gamma_{l} \stackrel{\rm ind}{\sim} \log N(C_{l}, U^2_{l})$$ $$C_{l}= -1$$ $$U^2_{l}= 10$$ Individual time effect $${\boldsymbol{\delta}}_{i} {\sim} GP(0, K)$$, $$i=1,\ldots,n$$ $${\rm e}^{-\theta} \sim \mbox{Beta}(A_\theta, B_\theta)$$ $$A_\theta= 1$$ $$B_\theta= 1$$ $$K$$ in (2.5) $${\boldsymbol{\Phi}} \sim \mbox{Wishart}({\boldsymbol{P}}_\Phi, N_\Phi)$$ $${\boldsymbol{P}}_\Phi= {\boldsymbol{I}}$$ $$N_\Phi=4$$ $$\mu_h \sim N(M_\mu, S^2_\mu)$$ $$M_\mu=0$$ $$S^2_\mu= 50$$ Medication level process $$z_{i}$$ in (2.6) $$\sigma^2 \sim \mbox{Gamma}^{-1}(A_\sigma, B_\sigma)$$ $$A_\sigma= 5$$ $$B_\sigma=.05$$ $${\rm e}^{-\phi} \sim \mbox{Beta}(A_\phi, B_\phi)$$ $$A_\phi= 6$$ $$B_\phi= 3$$ $$\rho \sim \mbox{Gamma}(A_\rho, B_\rho)$$ $$A_\rho= 0.25$$ $$B_\rho= 10$$ $$\varpi \sim \mbox{Beta}(A_\varpi, B_\varpi)$$ $$A_\varpi= 3$$ $$B_\varpi= 6$$ Table 1. Summary of the model for asthma related events defined in (2.2)–(2.6) and the specification of prior distributions Description Model Prior distributions Specification Individual effects $${\boldsymbol{\alpha}}_{i} \stackrel{\rm iid}{\sim} N({\boldsymbol{\nu}}, {\boldsymbol{\Sigma}})$$, $$i=1,\ldots,n$$ $$\nu_l \stackrel{\rm ind}{\sim} N(M_{\nu,l}, S^2_{\nu,l})$$ $$M_{\nu,l}= 0$$ $$S^2_{\nu,l}= 100$$ $${\boldsymbol{\Sigma}} \sim \mbox{Wishart}({\boldsymbol{P}}_\Sigma, N_\Sigma)$$ $${\boldsymbol{P}}_\Sigma= {\boldsymbol{I}}$$ $$N_\Sigma=4$$ Fixed effects $$\beta_{j,l}$$ in (2.3) $$\beta_{j,l} \stackrel{\rm ind}{\sim} N(B_{j,l}, T^2_{j,l})$$ $$B_{j,l}= 0$$ $$T^2_{j,l}= 100$$ Medication effects $$\gamma_{l}$$ in (2.3) $$-\gamma_{l} \stackrel{\rm ind}{\sim} \log N(C_{l}, U^2_{l})$$ $$C_{l}= -1$$ $$U^2_{l}= 10$$ Individual time effect $${\boldsymbol{\delta}}_{i} {\sim} GP(0, K)$$, $$i=1,\ldots,n$$ $${\rm e}^{-\theta} \sim \mbox{Beta}(A_\theta, B_\theta)$$ $$A_\theta= 1$$ $$B_\theta= 1$$ $$K$$ in (2.5) $${\boldsymbol{\Phi}} \sim \mbox{Wishart}({\boldsymbol{P}}_\Phi, N_\Phi)$$ $${\boldsymbol{P}}_\Phi= {\boldsymbol{I}}$$ $$N_\Phi=4$$ $$\mu_h \sim N(M_\mu, S^2_\mu)$$ $$M_\mu=0$$ $$S^2_\mu= 50$$ Medication level process $$z_{i}$$ in (2.6) $$\sigma^2 \sim \mbox{Gamma}^{-1}(A_\sigma, B_\sigma)$$ $$A_\sigma= 5$$ $$B_\sigma=.05$$ $${\rm e}^{-\phi} \sim \mbox{Beta}(A_\phi, B_\phi)$$ $$A_\phi= 6$$ $$B_\phi= 3$$ $$\rho \sim \mbox{Gamma}(A_\rho, B_\rho)$$ $$A_\rho= 0.25$$ $$B_\rho= 10$$ $$\varpi \sim \mbox{Beta}(A_\varpi, B_\varpi)$$ $$A_\varpi= 3$$ $$B_\varpi= 6$$ Description Model Prior distributions Specification Individual effects $${\boldsymbol{\alpha}}_{i} \stackrel{\rm iid}{\sim} N({\boldsymbol{\nu}}, {\boldsymbol{\Sigma}})$$, $$i=1,\ldots,n$$ $$\nu_l \stackrel{\rm ind}{\sim} N(M_{\nu,l}, S^2_{\nu,l})$$ $$M_{\nu,l}= 0$$ $$S^2_{\nu,l}= 100$$ $${\boldsymbol{\Sigma}} \sim \mbox{Wishart}({\boldsymbol{P}}_\Sigma, N_\Sigma)$$ $${\boldsymbol{P}}_\Sigma= {\boldsymbol{I}}$$ $$N_\Sigma=4$$ Fixed effects $$\beta_{j,l}$$ in (2.3) $$\beta_{j,l} \stackrel{\rm ind}{\sim} N(B_{j,l}, T^2_{j,l})$$ $$B_{j,l}= 0$$ $$T^2_{j,l}= 100$$ Medication effects $$\gamma_{l}$$ in (2.3) $$-\gamma_{l} \stackrel{\rm ind}{\sim} \log N(C_{l}, U^2_{l})$$ $$C_{l}= -1$$ $$U^2_{l}= 10$$ Individual time effect $${\boldsymbol{\delta}}_{i} {\sim} GP(0, K)$$, $$i=1,\ldots,n$$ $${\rm e}^{-\theta} \sim \mbox{Beta}(A_\theta, B_\theta)$$ $$A_\theta= 1$$ $$B_\theta= 1$$ $$K$$ in (2.5) $${\boldsymbol{\Phi}} \sim \mbox{Wishart}({\boldsymbol{P}}_\Phi, N_\Phi)$$ $${\boldsymbol{P}}_\Phi= {\boldsymbol{I}}$$ $$N_\Phi=4$$ $$\mu_h \sim N(M_\mu, S^2_\mu)$$ $$M_\mu=0$$ $$S^2_\mu= 50$$ Medication level process $$z_{i}$$ in (2.6) $$\sigma^2 \sim \mbox{Gamma}^{-1}(A_\sigma, B_\sigma)$$ $$A_\sigma= 5$$ $$B_\sigma=.05$$ $${\rm e}^{-\phi} \sim \mbox{Beta}(A_\phi, B_\phi)$$ $$A_\phi= 6$$ $$B_\phi= 3$$ $$\rho \sim \mbox{Gamma}(A_\rho, B_\rho)$$ $$A_\rho= 0.25$$ $$B_\rho= 10$$ $$\varpi \sim \mbox{Beta}(A_\varpi, B_\varpi)$$ $$A_\varpi= 3$$ $$B_\varpi= 6$$ 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= 0.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= 0.1$$, $$A_\sigma= 0.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 \begin{equation*} \Theta = \left\{ \{{\boldsymbol{\alpha}}_i\}_{i=1}^n, \{{\boldsymbol{\delta}}_i\}_{i=1}^n, \{\tilde{z}_{i,1}, \dots, \tilde{z}_{i,H}\}_{i=1}^n, \{B_i\}_{i=1}^n, \{{\boldsymbol{\beta}}_j\}_{j=1}^J, {\boldsymbol{\gamma}}, {\boldsymbol{\Sigma}}, \theta, {\boldsymbol{\Psi}}, \mu, \sigma^2, \phi, \varpi \right\}, \end{equation*} where $${\boldsymbol{\beta}}_j=[\beta_{j,1},\beta_{j,2},\beta_{j,3}]'$$, $${\boldsymbol{\gamma}}=[\gamma_{1},\gamma_{2},\gamma_{3}]'$$, 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 and others (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})]'$$ 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 (2.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 data set with $$n=4235$$ individuals ($$\sim$$300 000 total time periods), these 20 000 iterations took $$\sim$$40h. 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 $$\{ {\boldsymbol{\alpha}}_i, {\boldsymbol{\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, (respectively). 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. 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 )\!\!$$ 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 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 )\!\!$$ 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 )\!\!$$ 4.2 Prediction of an adverse event Let an adverse event for a 6-month period be defined as any ED/IP visit, any OCS use, or more than three 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 (14th individual in the data set), was at medication levels of 3, 2, 4, 5, and 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. 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 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 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 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 Conversely, Patient 2 (24th 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 6 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. 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 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 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 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 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 to 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 and others (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. Fig. 2. View largeDownload slide (a) Histogram of predicted probabilities for an adverse event in round five for the 500 held out observations. (b) Calibration plot of the predicted probability versus the observed probability (aggregated within different strata of risk). (c) ROC plot of the predicted probability versus observed. Fig. 2. View largeDownload slide (a) Histogram of predicted probabilities for an adverse event in round five for the 500 held out observations. (b) Calibration plot of the predicted probability versus the observed probability (aggregated within different strata of risk). (c) ROC plot of the predicted probability versus observed. 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 Figure 3(a)–(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. Fig. 3. View largeDownload slide Estimates from the proposed latent process approach and a simpler GLMER model obtained from the 100 simulated data sets. (a) Boxplots of estimates for $$y_1\::\:$$ED/IP Visit parameters. (b) Boxplots of estimates for $$y_2\::\:$$Rescue Inhaler Refills parameters. (c) Boxplots of estimates for $$y_3\::\:$$Use of OCS parameters. Fig. 3. View largeDownload slide Estimates from the proposed latent process approach and a simpler GLMER model obtained from the 100 simulated data sets. (a) Boxplots of estimates for $$y_1\::\:$$ED/IP Visit parameters. (b) Boxplots of estimates for $$y_2\::\:$$Rescue Inhaler Refills parameters. (c) Boxplots of estimates for $$y_3\::\:$$Use of OCS parameters. The estimates provided in Figure 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 $${\boldsymbol{\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 me4 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 Figure 3(a)–(c) as well. The true value of each parameter used to generate the 100 data sets 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. Fig. 4. View largeDownload slide (a) Calibration plot of the predicted probability versus the true probability of an adverse event (i.e., either an ED/IP visit, $$\geq 4$$ rescue inhaler fills, or any OCS use) in the next 6 months across all 100 simulations (aggregated within different strata of risk). (b) ROC plot using predicted probability to predict observed events (generated randomly according to the true probabilities). Fig. 4. View largeDownload slide (a) Calibration plot of the predicted probability versus the true probability of an adverse event (i.e., either an ED/IP visit, $$\geq 4$$ rescue inhaler fills, or any OCS use) in the next 6 months across all 100 simulations (aggregated within different strata of risk). (b) ROC plot using predicted probability to predict observed events (generated randomly according to the true probabilities). 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 6 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 6 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. Table 5. $$R^2$$ of the predicted probability of an event in the next 6 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 $$\geq 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) Method Outcome ED/IP OCS $$\geq 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) Table 5. $$R^2$$ of the predicted probability of an event in the next 6 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 $$\geq 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) Method Outcome ED/IP OCS $$\geq 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) 6. Conclusions and further work In this article, 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 health care 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. References Austin, P. C. and Steyerberg, E. W. ( 2014 ). Graphical assessment of internal and external calibration of logistic regression models by using loess smoothers. Statistics in Medicine 33 , 517 – 535 . Google Scholar CrossRef Search ADS PubMed Bateman, E. D. , Buhl, R. , O’Byrne, P. M. , Humbert, M. , Reddel, H. K. , Sears, M. R. , Jenkins, C. , Harrison, T. W. , Quirce, S. , Peterson, S. and others. ( 2015 ). Development and validation of a novel risk score for asthma exacerbations: the risk score for exacerbations. Journal of Allergy and Clinical Immunology 135 , 1457 – 1464 . Google Scholar CrossRef Search ADS PubMed Blakey, J. D. , Woulnough, K. , James, A. C. , Fellows, J. , Obeidat, M. , Navaratnam, V. , Stringfellow, T. , Yeoh, Z. W. , Pavord, I. , Thomas, M. and others. ( 2012 ). S62 a systematic review of factors associated with future asthma attacks to inform a risk assessment questionnaire. Thorax 67 (Suppl 2) , A31 – A32 . Brix, A. and Diggle, P. J. ( 2001 ). Spatiotemporal prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 63 , 823 – 841 . Google Scholar CrossRef Search ADS Carroll, R. J , Ruppert, D. , Stefanski, L. A. and Crainiceanu, C. M. ( 2006 ). Measurement Error in Nonlinear Models: A Modern Perspective , 2nd edition . Boca Raton, FL : CRC Press . Google Scholar CrossRef Search ADS Dunson, D. B. ( 2000 ). Bayesian latent variable models for clustered mixed outcomes. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 62 , 355 – 366 . Google Scholar CrossRef Search ADS Harrell, F. E. , Lee, K. L. and Mark, D. B. ( 1996 ). Tutorial in biostatistics multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Statistics in Medicine 15 , 361 – 387 . Google Scholar CrossRef Search ADS PubMed Jang, J. , Chan, K. C. G. , Huang, H. and Sullivan, S. D. ( 2013 ). Trends in cost and outcomes among adult and pediatric patients with asthma: 2000–2009. Annals of Allergy, Asthma & Immunology 111 , 516 – 522 . Google Scholar CrossRef Search ADS Kamble, S. and Bharmal, M. ( 2009 ). Incremental direct expenditure of treating asthma in the United States. Journal of Asthma 46 , 73 – 80 . Google Scholar CrossRef Search ADS PubMed Møller, J. , Syversveen, A. R. and Waagepetersen, R. P. ( 1998 ). Log Gaussian Cox processes. Scandinavian Journal of Statistics 25 , 451 – 482 . Google Scholar CrossRef Search ADS Moorman, J. E. , Akinbami, L. J. , Bailey, C. M. , Zahran, H. S. , King, M. E. , Johnson, C. A. and Liu, X. ( 2012 ). National surveillance of asthma: United States, 2001–2010. Vital & Health Statistics. Series 3, Analytical and Epidemiological Studies[US Department of Health and Human Services, Public Health Service, National Center for Health Statistics] 35 , 1 – 58 . Muthen, B. ( 1983 ). Latent variable structural equation modeling with categorical data. Journal of Econometrics 22 , 43 – 65 . Google Scholar CrossRef Search ADS Plummer, M. ( 2009 ). rjags: interface to the jags MCMC library. R package version, Vol. 1. https://cran.r-project.org/. Rank, M. A. , Liesinger, J. T. , Branda, M. E. , Gionfriddo, M. R. , Schatz, M. , Zeiger, R. S. and Shah, N. D. ( 2016 ). Comparative safety and costs of stepping down asthma medications in patients with controlled asthma. Journal of Allergy and Clinical Immunology 137 , 1373 – 1379 . Google Scholar CrossRef Search ADS PubMed Rank, M. A. , Liesinger, J. T. , Ziegenfuss, J. Y. , Branda, M. E. , Lim, K. G. , Yawn, B. P. , Li, J. T. and Shah, N. D. ( 2012a ). Asthma expenditures in the United States comparing 2004 to 2006 and 1996 to 1998. The American Journal of Managed Care 18 , 499 – 504 . Rank, M. A. , Liesinger, J. T. , Ziegenfuss, J. Y. , Branda, M. E. , Lim, K. G. , Yawn, B. P. and Shah, N. D. ( 2012b ). The impact of asthma medication guidelines on asthma controller use and on asthma exacerbation rates comparing 1997–1998 and 2004–2005. Annals of Allergy, Asthma & Immunology 108 , 9 – 13 . Google Scholar CrossRef Search ADS Storlie, C. B. , Michalak, S. E. , Quinn, H. M. , Dubois, A. J. , Wender, S. A. and Dubois, D. H. ( 2013 ). A Bayesian reliability analysis of neutron-induced errors in high performance computing hardware. Journal of the American Statistical Association 108 , 429 – 440 . Google Scholar CrossRef Search ADS Tsiatis, A. A. and Davidian, M. ( 2001 ). A semiparametric estimator for the proportional hazards model with longitudinal covariates measured with error. Biometrika 88 , 447 – 458 . Google Scholar CrossRef Search ADS Wulfsohn, M. S. and Tsiatis, A. A. ( 1997 ). A joint model for survival and longitudinal data measured with error. Biometrics 53 , 330 – 339 . Google Scholar CrossRef Search ADS PubMed © The Author 2017. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Biostatistics Oxford University Press

# Prediction of individual outcomes for asthma sufferers

, Volume Advance Article – Nov 6, 2017
15 pages

Publisher
Oxford University Press
ISSN
1465-4644
eISSN
1468-4357
D.O.I.
10.1093/biostatistics/kxx055
Publisher site
See Article on Publisher Site

### Abstract

Summary 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 United States of America. 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. 1. Introduction Asthma is a chronic disease, affecting 8% of the population and costing$37–63 billion/year in the United States of America (Kamble and Bharmal, 2009; Jang and others, 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 health care spending for asthma. In addition, spending on medications has been growing at a rapid rate over the last decade (Rank and others, 2012a). While medication spending continues to increase, this has not resulted in improved asthma outcomes (Moorman and others, 2012; Rank and others, 2012b). 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 and others, 2015; Blakey and others, 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 $$4235$$ 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 five rounds (periods of 4–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 and others (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 6 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, and (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 and others, 2013). Conditional on the latent medication process, the model for the multivariate response over time is a dependent log-Gaussian Cox Process (Møller and others, 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 and others, 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). The remainder of the article 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, 1 week granularity was used in all results below for computational expediency. The response for the $$i$$th individual for period $$n$$ (i.e., week $$n$$) is $${\boldsymbol{y}}_{i,n} = [y_{i,n,1}, y_{i,n,2}, y_{i,n,3}]'$$, with \begin{eqnarray} y_{i,n,1} & = & \mbox{# of ED/IP visits for ith individual during nth time period.} \nonumber\\ y_{i,n,2} & = & \mbox{# of rescue inhaler refills for ith individual during nth time period.} \\ y_{i,n,3} & = & \mbox{# of oral steroid refills for ith individual during nth time period.} \nonumber \end{eqnarray} (2.1) 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$$th round are observed, i.e., \begin{eqnarray*} y^*_{i,m,2} &=& \sum_{n \in R_{i,m}} y_{i,n,2}, \nonumber\\ y^*_{i,m,3} &=& \sum_{n \in R_{i,m}} y_{i,n,3}, \end{eqnarray*} where $$R_{i,m}$$ is the set of time periods that make up the $$m$$th round. The distribution of these responses was assumed to be a function of patient dependent covariates $${\boldsymbol{x}}_{i} = [x_{i,1},\dots,x_{i,p}]'$$ (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 and others (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 $${\boldsymbol{y}}_{i,n}$$ also depends on a time varying random effect for the $$i$$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, $$y_{i,n,l} \sim \mbox{Poisson}(\lambda_{i,n,l}),$$ (2.2) where \begin{equation*} \lambda_{i,n,l} = \int_{t_{n-1}}^{t_n} \eta_{i,l}(s) {\rm d}s, \end{equation*} where $$t_n$$ is the end of the $$n$$th time period and $$\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),$$ (2.3) 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$$th response type, to account for the overall level of severity for the $$i$$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$$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, \begin{eqnarray*} y^*_{i,m,2} \sim \mbox{Poisson}\left(\sum_{n \in R_{i,m}} \! \lambda_{i,n,2}\right)\!, \nonumber \\ y^*_{i,m,3} \sim \mbox{Poisson}\left(\sum_{n \in R_{i,m}} \! \lambda_{i,n,3}\right)\!. \end{eqnarray*} 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}$$, $${\boldsymbol{\alpha}}_i = \left[\alpha_{i,1},\dots,\alpha_{i,L} \right]' \sim N \left( {\boldsymbol{\nu}}, {\boldsymbol{\Sigma}} \right)\!.$$ (2.4) 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 (O–U) process (equivalent to an order one, autoregressive time series) for each $$l$$, i.e., \begin{equation*} K_l(s,t) = \tau_l \exp \{-\theta |s-t|\}. \end{equation*} A product correlation was also assumed between $$\delta_{i,l}$$ and $$\delta_{i,l'}$$ so that the variation in the rate of occurrence for the three response types are correlated, i.e., \begin{equation*} {\text{Cov}} \left( \left[\delta_{i,1}(s), \dots, \delta_{i,L}(s) \right]' \right) = {\boldsymbol{\Phi}}, \end{equation*} or the general covariance function is then $$K\left((l,s),(l',t) \right) = {\text{Cov}} \left( \delta_{i,l}(s), \delta_{i,l'}(t) \right) = {\boldsymbol{\Phi}}_{l,l'} \exp \{-\theta |s-t|\}.$$ (2.5) 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 6 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 $$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\}},$$ (2.6) where (i) $$\tilde{z}_{i,h}$$, $$h=0,1,\dots$$, are iid mean zero 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 (2.6) ensures that $$z_i(t)\geq0$$), (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, \begin{equation*} B_i(t) = C_i(t) + D_i(t), \end{equation*} 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 \begin{equation*} \begin{array}{rcl} \Pr \{ D_i(t^+) = D_i(t)+1 \} & = & \varpi, \\[.1in] \Pr \{ D_i(t^+) = D_i(t) \} & = & 1- \varpi. \end{array} \end{equation*} 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$$th round. In the above model, the latent variable $$z_i$$ is continuous even though the observed value for the $$m$$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 cut-points 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 \begin{equation*} z^*_{i,m} = \min\left\{ 5\;,\; \left\lfloor \int_{t \in R_{i,m}} z_i(t){\rm d}t + 0.5 \right\rfloor \right\}\!. \end{equation*} 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 as 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. Fig. 1. View largeDownload slide Twenty-five posterior realizations of the latent controller medication process, $$z_i$$ for a given patient. When an increase in medication occurs, the model generally prefers to have it occur at ED/IP/OP/OB visits over other times. Fig. 1. View largeDownload slide Twenty-five posterior realizations of the latent controller medication process, $$z_i$$ for a given patient. When an increase in medication occurs, the model generally prefers to have it occur at ED/IP/OP/OB visits over other times. 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 6 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. Table 1. Summary of the model for asthma related events defined in (2.2)–(2.6) and the specification of prior distributions Description Model Prior distributions Specification Individual effects $${\boldsymbol{\alpha}}_{i} \stackrel{\rm iid}{\sim} N({\boldsymbol{\nu}}, {\boldsymbol{\Sigma}})$$, $$i=1,\ldots,n$$ $$\nu_l \stackrel{\rm ind}{\sim} N(M_{\nu,l}, S^2_{\nu,l})$$ $$M_{\nu,l}= 0$$ $$S^2_{\nu,l}= 100$$ $${\boldsymbol{\Sigma}} \sim \mbox{Wishart}({\boldsymbol{P}}_\Sigma, N_\Sigma)$$ $${\boldsymbol{P}}_\Sigma= {\boldsymbol{I}}$$ $$N_\Sigma=4$$ Fixed effects $$\beta_{j,l}$$ in (2.3) $$\beta_{j,l} \stackrel{\rm ind}{\sim} N(B_{j,l}, T^2_{j,l})$$ $$B_{j,l}= 0$$ $$T^2_{j,l}= 100$$ Medication effects $$\gamma_{l}$$ in (2.3) $$-\gamma_{l} \stackrel{\rm ind}{\sim} \log N(C_{l}, U^2_{l})$$ $$C_{l}= -1$$ $$U^2_{l}= 10$$ Individual time effect $${\boldsymbol{\delta}}_{i} {\sim} GP(0, K)$$, $$i=1,\ldots,n$$ $${\rm e}^{-\theta} \sim \mbox{Beta}(A_\theta, B_\theta)$$ $$A_\theta= 1$$ $$B_\theta= 1$$ $$K$$ in (2.5) $${\boldsymbol{\Phi}} \sim \mbox{Wishart}({\boldsymbol{P}}_\Phi, N_\Phi)$$ $${\boldsymbol{P}}_\Phi= {\boldsymbol{I}}$$ $$N_\Phi=4$$ $$\mu_h \sim N(M_\mu, S^2_\mu)$$ $$M_\mu=0$$ $$S^2_\mu= 50$$ Medication level process $$z_{i}$$ in (2.6) $$\sigma^2 \sim \mbox{Gamma}^{-1}(A_\sigma, B_\sigma)$$ $$A_\sigma= 5$$ $$B_\sigma=.05$$ $${\rm e}^{-\phi} \sim \mbox{Beta}(A_\phi, B_\phi)$$ $$A_\phi= 6$$ $$B_\phi= 3$$ $$\rho \sim \mbox{Gamma}(A_\rho, B_\rho)$$ $$A_\rho= 0.25$$ $$B_\rho= 10$$ $$\varpi \sim \mbox{Beta}(A_\varpi, B_\varpi)$$ $$A_\varpi= 3$$ $$B_\varpi= 6$$ Description Model Prior distributions Specification Individual effects $${\boldsymbol{\alpha}}_{i} \stackrel{\rm iid}{\sim} N({\boldsymbol{\nu}}, {\boldsymbol{\Sigma}})$$, $$i=1,\ldots,n$$ $$\nu_l \stackrel{\rm ind}{\sim} N(M_{\nu,l}, S^2_{\nu,l})$$ $$M_{\nu,l}= 0$$ $$S^2_{\nu,l}= 100$$ $${\boldsymbol{\Sigma}} \sim \mbox{Wishart}({\boldsymbol{P}}_\Sigma, N_\Sigma)$$ $${\boldsymbol{P}}_\Sigma= {\boldsymbol{I}}$$ $$N_\Sigma=4$$ Fixed effects $$\beta_{j,l}$$ in (2.3) $$\beta_{j,l} \stackrel{\rm ind}{\sim} N(B_{j,l}, T^2_{j,l})$$ $$B_{j,l}= 0$$ $$T^2_{j,l}= 100$$ Medication effects $$\gamma_{l}$$ in (2.3) $$-\gamma_{l} \stackrel{\rm ind}{\sim} \log N(C_{l}, U^2_{l})$$ $$C_{l}= -1$$ $$U^2_{l}= 10$$ Individual time effect $${\boldsymbol{\delta}}_{i} {\sim} GP(0, K)$$, $$i=1,\ldots,n$$ $${\rm e}^{-\theta} \sim \mbox{Beta}(A_\theta, B_\theta)$$ $$A_\theta= 1$$ $$B_\theta= 1$$ $$K$$ in (2.5) $${\boldsymbol{\Phi}} \sim \mbox{Wishart}({\boldsymbol{P}}_\Phi, N_\Phi)$$ $${\boldsymbol{P}}_\Phi= {\boldsymbol{I}}$$ $$N_\Phi=4$$ $$\mu_h \sim N(M_\mu, S^2_\mu)$$ $$M_\mu=0$$ $$S^2_\mu= 50$$ Medication level process $$z_{i}$$ in (2.6) $$\sigma^2 \sim \mbox{Gamma}^{-1}(A_\sigma, B_\sigma)$$ $$A_\sigma= 5$$ $$B_\sigma=.05$$ $${\rm e}^{-\phi} \sim \mbox{Beta}(A_\phi, B_\phi)$$ $$A_\phi= 6$$ $$B_\phi= 3$$ $$\rho \sim \mbox{Gamma}(A_\rho, B_\rho)$$ $$A_\rho= 0.25$$ $$B_\rho= 10$$ $$\varpi \sim \mbox{Beta}(A_\varpi, B_\varpi)$$ $$A_\varpi= 3$$ $$B_\varpi= 6$$ Table 1. Summary of the model for asthma related events defined in (2.2)–(2.6) and the specification of prior distributions Description Model Prior distributions Specification Individual effects $${\boldsymbol{\alpha}}_{i} \stackrel{\rm iid}{\sim} N({\boldsymbol{\nu}}, {\boldsymbol{\Sigma}})$$, $$i=1,\ldots,n$$ $$\nu_l \stackrel{\rm ind}{\sim} N(M_{\nu,l}, S^2_{\nu,l})$$ $$M_{\nu,l}= 0$$ $$S^2_{\nu,l}= 100$$ $${\boldsymbol{\Sigma}} \sim \mbox{Wishart}({\boldsymbol{P}}_\Sigma, N_\Sigma)$$ $${\boldsymbol{P}}_\Sigma= {\boldsymbol{I}}$$ $$N_\Sigma=4$$ Fixed effects $$\beta_{j,l}$$ in (2.3) $$\beta_{j,l} \stackrel{\rm ind}{\sim} N(B_{j,l}, T^2_{j,l})$$ $$B_{j,l}= 0$$ $$T^2_{j,l}= 100$$ Medication effects $$\gamma_{l}$$ in (2.3) $$-\gamma_{l} \stackrel{\rm ind}{\sim} \log N(C_{l}, U^2_{l})$$ $$C_{l}= -1$$ $$U^2_{l}= 10$$ Individual time effect $${\boldsymbol{\delta}}_{i} {\sim} GP(0, K)$$, $$i=1,\ldots,n$$ $${\rm e}^{-\theta} \sim \mbox{Beta}(A_\theta, B_\theta)$$ $$A_\theta= 1$$ $$B_\theta= 1$$ $$K$$ in (2.5) $${\boldsymbol{\Phi}} \sim \mbox{Wishart}({\boldsymbol{P}}_\Phi, N_\Phi)$$ $${\boldsymbol{P}}_\Phi= {\boldsymbol{I}}$$ $$N_\Phi=4$$ $$\mu_h \sim N(M_\mu, S^2_\mu)$$ $$M_\mu=0$$ $$S^2_\mu= 50$$ Medication level process $$z_{i}$$ in (2.6) $$\sigma^2 \sim \mbox{Gamma}^{-1}(A_\sigma, B_\sigma)$$ $$A_\sigma= 5$$ $$B_\sigma=.05$$ $${\rm e}^{-\phi} \sim \mbox{Beta}(A_\phi, B_\phi)$$ $$A_\phi= 6$$ $$B_\phi= 3$$ $$\rho \sim \mbox{Gamma}(A_\rho, B_\rho)$$ $$A_\rho= 0.25$$ $$B_\rho= 10$$ $$\varpi \sim \mbox{Beta}(A_\varpi, B_\varpi)$$ $$A_\varpi= 3$$ $$B_\varpi= 6$$ Description Model Prior distributions Specification Individual effects $${\boldsymbol{\alpha}}_{i} \stackrel{\rm iid}{\sim} N({\boldsymbol{\nu}}, {\boldsymbol{\Sigma}})$$, $$i=1,\ldots,n$$ $$\nu_l \stackrel{\rm ind}{\sim} N(M_{\nu,l}, S^2_{\nu,l})$$ $$M_{\nu,l}= 0$$ $$S^2_{\nu,l}= 100$$ $${\boldsymbol{\Sigma}} \sim \mbox{Wishart}({\boldsymbol{P}}_\Sigma, N_\Sigma)$$ $${\boldsymbol{P}}_\Sigma= {\boldsymbol{I}}$$ $$N_\Sigma=4$$ Fixed effects $$\beta_{j,l}$$ in (2.3) $$\beta_{j,l} \stackrel{\rm ind}{\sim} N(B_{j,l}, T^2_{j,l})$$ $$B_{j,l}= 0$$ $$T^2_{j,l}= 100$$ Medication effects $$\gamma_{l}$$ in (2.3) $$-\gamma_{l} \stackrel{\rm ind}{\sim} \log N(C_{l}, U^2_{l})$$ $$C_{l}= -1$$ $$U^2_{l}= 10$$ Individual time effect $${\boldsymbol{\delta}}_{i} {\sim} GP(0, K)$$, $$i=1,\ldots,n$$ $${\rm e}^{-\theta} \sim \mbox{Beta}(A_\theta, B_\theta)$$ $$A_\theta= 1$$ $$B_\theta= 1$$ $$K$$ in (2.5) $${\boldsymbol{\Phi}} \sim \mbox{Wishart}({\boldsymbol{P}}_\Phi, N_\Phi)$$ $${\boldsymbol{P}}_\Phi= {\boldsymbol{I}}$$ $$N_\Phi=4$$ $$\mu_h \sim N(M_\mu, S^2_\mu)$$ $$M_\mu=0$$ $$S^2_\mu= 50$$ Medication level process $$z_{i}$$ in (2.6) $$\sigma^2 \sim \mbox{Gamma}^{-1}(A_\sigma, B_\sigma)$$ $$A_\sigma= 5$$ $$B_\sigma=.05$$ $${\rm e}^{-\phi} \sim \mbox{Beta}(A_\phi, B_\phi)$$ $$A_\phi= 6$$ $$B_\phi= 3$$ $$\rho \sim \mbox{Gamma}(A_\rho, B_\rho)$$ $$A_\rho= 0.25$$ $$B_\rho= 10$$ $$\varpi \sim \mbox{Beta}(A_\varpi, B_\varpi)$$ $$A_\varpi= 3$$ $$B_\varpi= 6$$ 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= 0.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= 0.1$$, $$A_\sigma= 0.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 \begin{equation*} \Theta = \left\{ \{{\boldsymbol{\alpha}}_i\}_{i=1}^n, \{{\boldsymbol{\delta}}_i\}_{i=1}^n, \{\tilde{z}_{i,1}, \dots, \tilde{z}_{i,H}\}_{i=1}^n, \{B_i\}_{i=1}^n, \{{\boldsymbol{\beta}}_j\}_{j=1}^J, {\boldsymbol{\gamma}}, {\boldsymbol{\Sigma}}, \theta, {\boldsymbol{\Psi}}, \mu, \sigma^2, \phi, \varpi \right\}, \end{equation*} where $${\boldsymbol{\beta}}_j=[\beta_{j,1},\beta_{j,2},\beta_{j,3}]'$$, $${\boldsymbol{\gamma}}=[\gamma_{1},\gamma_{2},\gamma_{3}]'$$, 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 and others (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})]'$$ 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 (2.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 data set with $$n=4235$$ individuals ($$\sim$$300 000 total time periods), these 20 000 iterations took $$\sim$$40h. 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 $$\{ {\boldsymbol{\alpha}}_i, {\boldsymbol{\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, (respectively). 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. 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 )\!\!$$ 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 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 )\!\!$$ 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 )\!\!$$ 4.2 Prediction of an adverse event Let an adverse event for a 6-month period be defined as any ED/IP visit, any OCS use, or more than three 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 (14th individual in the data set), was at medication levels of 3, 2, 4, 5, and 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. 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 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 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 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 Conversely, Patient 2 (24th 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 6 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. 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 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 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 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 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 to 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 and others (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. Fig. 2. View largeDownload slide (a) Histogram of predicted probabilities for an adverse event in round five for the 500 held out observations. (b) Calibration plot of the predicted probability versus the observed probability (aggregated within different strata of risk). (c) ROC plot of the predicted probability versus observed. Fig. 2. View largeDownload slide (a) Histogram of predicted probabilities for an adverse event in round five for the 500 held out observations. (b) Calibration plot of the predicted probability versus the observed probability (aggregated within different strata of risk). (c) ROC plot of the predicted probability versus observed. 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 Figure 3(a)–(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. Fig. 3. View largeDownload slide Estimates from the proposed latent process approach and a simpler GLMER model obtained from the 100 simulated data sets. (a) Boxplots of estimates for $$y_1\::\:$$ED/IP Visit parameters. (b) Boxplots of estimates for $$y_2\::\:$$Rescue Inhaler Refills parameters. (c) Boxplots of estimates for $$y_3\::\:$$Use of OCS parameters. Fig. 3. View largeDownload slide Estimates from the proposed latent process approach and a simpler GLMER model obtained from the 100 simulated data sets. (a) Boxplots of estimates for $$y_1\::\:$$ED/IP Visit parameters. (b) Boxplots of estimates for $$y_2\::\:$$Rescue Inhaler Refills parameters. (c) Boxplots of estimates for $$y_3\::\:$$Use of OCS parameters. The estimates provided in Figure 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 $${\boldsymbol{\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 me4 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 Figure 3(a)–(c) as well. The true value of each parameter used to generate the 100 data sets 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. Fig. 4. View largeDownload slide (a) Calibration plot of the predicted probability versus the true probability of an adverse event (i.e., either an ED/IP visit, $$\geq 4$$ rescue inhaler fills, or any OCS use) in the next 6 months across all 100 simulations (aggregated within different strata of risk). (b) ROC plot using predicted probability to predict observed events (generated randomly according to the true probabilities). Fig. 4. View largeDownload slide (a) Calibration plot of the predicted probability versus the true probability of an adverse event (i.e., either an ED/IP visit, $$\geq 4$$ rescue inhaler fills, or any OCS use) in the next 6 months across all 100 simulations (aggregated within different strata of risk). (b) ROC plot using predicted probability to predict observed events (generated randomly according to the true probabilities). 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 6 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 6 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. Table 5. $$R^2$$ of the predicted probability of an event in the next 6 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 $$\geq 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) Method Outcome ED/IP OCS $$\geq 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) Table 5. $$R^2$$ of the predicted probability of an event in the next 6 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 $$\geq 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) Method Outcome ED/IP OCS $$\geq 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) 6. Conclusions and further work In this article, 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 health care 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. References Austin, P. C. and Steyerberg, E. W. ( 2014 ). Graphical assessment of internal and external calibration of logistic regression models by using loess smoothers. Statistics in Medicine 33 , 517 – 535 . Google Scholar CrossRef Search ADS PubMed Bateman, E. D. , Buhl, R. , O’Byrne, P. M. , Humbert, M. , Reddel, H. K. , Sears, M. R. , Jenkins, C. , Harrison, T. W. , Quirce, S. , Peterson, S. and others. ( 2015 ). Development and validation of a novel risk score for asthma exacerbations: the risk score for exacerbations. Journal of Allergy and Clinical Immunology 135 , 1457 – 1464 . Google Scholar CrossRef Search ADS PubMed Blakey, J. D. , Woulnough, K. , James, A. C. , Fellows, J. , Obeidat, M. , Navaratnam, V. , Stringfellow, T. , Yeoh, Z. W. , Pavord, I. , Thomas, M. and others. ( 2012 ). S62 a systematic review of factors associated with future asthma attacks to inform a risk assessment questionnaire. Thorax 67 (Suppl 2) , A31 – A32 . Brix, A. and Diggle, P. J. ( 2001 ). Spatiotemporal prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 63 , 823 – 841 . Google Scholar CrossRef Search ADS Carroll, R. J , Ruppert, D. , Stefanski, L. A. and Crainiceanu, C. M. ( 2006 ). Measurement Error in Nonlinear Models: A Modern Perspective , 2nd edition . Boca Raton, FL : CRC Press . Google Scholar CrossRef Search ADS Dunson, D. B. ( 2000 ). Bayesian latent variable models for clustered mixed outcomes. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 62 , 355 – 366 . Google Scholar CrossRef Search ADS Harrell, F. E. , Lee, K. L. and Mark, D. B. ( 1996 ). Tutorial in biostatistics multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Statistics in Medicine 15 , 361 – 387 . Google Scholar CrossRef Search ADS PubMed Jang, J. , Chan, K. C. G. , Huang, H. and Sullivan, S. D. ( 2013 ). Trends in cost and outcomes among adult and pediatric patients with asthma: 2000–2009. Annals of Allergy, Asthma & Immunology 111 , 516 – 522 . Google Scholar CrossRef Search ADS Kamble, S. and Bharmal, M. ( 2009 ). Incremental direct expenditure of treating asthma in the United States. Journal of Asthma 46 , 73 – 80 . Google Scholar CrossRef Search ADS PubMed Møller, J. , Syversveen, A. R. and Waagepetersen, R. P. ( 1998 ). Log Gaussian Cox processes. Scandinavian Journal of Statistics 25 , 451 – 482 . Google Scholar CrossRef Search ADS Moorman, J. E. , Akinbami, L. J. , Bailey, C. M. , Zahran, H. S. , King, M. E. , Johnson, C. A. and Liu, X. ( 2012 ). National surveillance of asthma: United States, 2001–2010. Vital & Health Statistics. Series 3, Analytical and Epidemiological Studies[US Department of Health and Human Services, Public Health Service, National Center for Health Statistics] 35 , 1 – 58 . Muthen, B. ( 1983 ). Latent variable structural equation modeling with categorical data. Journal of Econometrics 22 , 43 – 65 . Google Scholar CrossRef Search ADS Plummer, M. ( 2009 ). rjags: interface to the jags MCMC library. R package version, Vol. 1. https://cran.r-project.org/. Rank, M. A. , Liesinger, J. T. , Branda, M. E. , Gionfriddo, M. R. , Schatz, M. , Zeiger, R. S. and Shah, N. D. ( 2016 ). Comparative safety and costs of stepping down asthma medications in patients with controlled asthma. Journal of Allergy and Clinical Immunology 137 , 1373 – 1379 . Google Scholar CrossRef Search ADS PubMed Rank, M. A. , Liesinger, J. T. , Ziegenfuss, J. Y. , Branda, M. E. , Lim, K. G. , Yawn, B. P. , Li, J. T. and Shah, N. D. ( 2012a ). Asthma expenditures in the United States comparing 2004 to 2006 and 1996 to 1998. The American Journal of Managed Care 18 , 499 – 504 . Rank, M. A. , Liesinger, J. T. , Ziegenfuss, J. Y. , Branda, M. E. , Lim, K. G. , Yawn, B. P. and Shah, N. D. ( 2012b ). The impact of asthma medication guidelines on asthma controller use and on asthma exacerbation rates comparing 1997–1998 and 2004–2005. Annals of Allergy, Asthma & Immunology 108 , 9 – 13 . Google Scholar CrossRef Search ADS Storlie, C. B. , Michalak, S. E. , Quinn, H. M. , Dubois, A. J. , Wender, S. A. and Dubois, D. H. ( 2013 ). A Bayesian reliability analysis of neutron-induced errors in high performance computing hardware. Journal of the American Statistical Association 108 , 429 – 440 . Google Scholar CrossRef Search ADS Tsiatis, A. A. and Davidian, M. ( 2001 ). A semiparametric estimator for the proportional hazards model with longitudinal covariates measured with error. Biometrika 88 , 447 – 458 . Google Scholar CrossRef Search ADS Wulfsohn, M. S. and Tsiatis, A. A. ( 1997 ). A joint model for survival and longitudinal data measured with error. Biometrics 53 , 330 – 339 . Google Scholar CrossRef Search ADS PubMed © The Author 2017. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com.

### Journal

BiostatisticsOxford University Press

Published: Nov 6, 2017

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations

Abstract access only

18 million full-text articles

Print

20 pages / month

PDF Discount

20% off