# Measuring Effects of Medication Adherence on Time-Varying Health Outcomes
using Bayesian Dynamic Linear Models††Work supported by NIH grant R21 HL121366

Luis F. Campos, Mark E. Glickman, Kristen B. Hunter

###### Abstract

One of the most significant barriers to medication treatment is patients’ non-
adherence to a prescribed medication regimen. The extent of the impact of poor
adherence on resulting health measures is often unknown, and typical analyses
ignore the time-varying nature of adherence. This paper develops a modeling
framework for longitudinally recorded health measures modeled as a function of
time-varying medication adherence or other time-varying covariates. Our
framework, which relies on normal Bayesian dynamic linear models (DLMs),
accounts for time-varying covariates such as adherence and non-dynamic
covariates such as baseline health characteristics. Given the inefficiencies
using standard inferential procedures for DLMs associated with infrequent and
irregularly recorded response data, we develop an approach that relies on
factoring the posterior density into a product of two terms; a marginal
posterior density for the non-dynamic parameters, and a multivariate normal
posterior density of the dynamic parameters conditional on the non-dynamic
ones. This factorization leads to a two-stage process for inference in which
the non-dynamic parameters can be inferred separately from the time-varying
parameters. We demonstrate the application of this model to the time-varying
effect of anti-hypertensive medication on blood pressure levels from a cohort
of patients diagnosed with hypertension. Our model results are compared to
ones in which adherence is incorporated through non-dynamic summaries.

##  1 Introduction

Over 85 million American adults, or about one third of the population over 20
years old, suffer from hypertension (Benjamin et al.,, 2018). Approximately
16% of adults in the United States are unaware that they have hypertension
(Benjamin et al.,, 2018). Left untreated, hypertension can lead to a range of
serious and costly health concerns such as cardiovascular disease, stroke, and
renal disease (Amery et al.,, 1985; Probstfield,, 1991). Among the many
factors associated with uncontrolled BP, poor adherence to prescribed anti-
hypertensive medications is of increasing concern to clinicians, health care
systems, and other stakeholders (Choo et al.,, 1999; Morisky et al.,, 1986;
Osterberg and Blaschke,, 2005). Little doubt exists that patients who adhere
poorly to their prescribed medication are at risk for worse BP outcomes. Given
the large variation in adherence, both across patients and temporally within
patients, it is an open question how to accurately measure the impact of
varying adherence patterns on BP levels. Furthermore, the variation in
adherence patterns creates difficulties in accurately measuring the effects of
socio-demographic and health risk factors on BP levels.

This paper develops a Bayesian dynamic linear model (Durbin and Koopman,,
2001; West and Harrison,, 1997; Petris et al.,, 2009) for health measures
recorded over time as a function of time-varying adherence, with a particular
application to the effects of anti-hypertensive medication on BP levels.
Bayesian dynamic linear models (DLMs) have a long history as a statistical
framework for forecasting and measuring trajectories in many domains,
including real-time missile tracking and financial securities forecasting, but
rarely in the health context. We apply DLMs to describe time-varying health
measures (like blood pressure levels) as a function of detailed adherence (or
other time-varying covariates) and individual demographics and comorbidities
measured typically at study baseline. The application of DLMs to the time-
varying adherence on health measures is novel, but fits naturally into the DLM
framework because these measures can be “tracked” over time as adherence data
accumulate. Because the DLM framework permits the inclusion of patient-level
predictors, the model can be applied to measuring effects of socio-economic or
racial disparities, or the effects of different comorbid conditions. In these
settings, our model can control for differential medication adherence
patterns, resulting in more accurate measurement of the effects of other
covariates.

Data for studies in which the estimated effects of time-varying adherence on
health measures are desired usually consist of three components. First,
adherence data for each study participant are assumed to be collected through
electronic monitoring devices. These devices electronically time-stamp each
time the pill container is opened, and permit an accurate recording of when a
patient took their medication. Second, health and socio-demographic are
typically recorded at the start of medication adherence studies, and are often
non-dynamic. Finally, health measures which might be impacted by differential
medication adherence are recorded over time at intermittent intervals. Such
measures are often recorded at clinic visits, the timing of which may be
determined by the patient. Thus it is quite common for the number of health
measures per study patient to be much fewer than the number of days in which
medication adherence information is recorded.

The DLM framework is challenged by having much fewer days with health measures
compared to the number of days of adherence measurements. In typical uses of
DLMs, both time-varying covariates and responses are measured at every time
point, and inference for the time-varying state parameters can be accomplished
using standard Bayesian updating algorithms (West and Harrison,, 1997, Chapter
16). With adherence data, when the responses are measured at irregular and
infrequent intervals, the usual updating approaches can be demonstrably
inefficient. Instead, we develop an inferential approach in a Bayesian setting
that takes advantage of factoring the posterior density into a product of two
terms. The first term involves the exact marginal likelihood of the DLM,
marginalizing out the dynamic state process. The second is the conditional
posterior density for the state process parameters. These together allow us to
determine the posterior distribution for the non-dynamic parameters using
standard Markov chain Monte Carlo (MCMC) techniques as a first stage and the
state process parameters as a second stage without resorting to needlessly
complex computational tools. With a DLM that has normally distributed
responses and a stochastic process for the latent states the has normal
innovations, the factorization is the product of two multivariate normal
densities. This factored posterior density also easily permits inference for
the non-dynamic parameters, which in the setting of medication adherence is
likely of interest because the researcher may want to understand the effect of
baseline health characteristics on the health measures controlling for
differential adherence.

The remainder of this paper is organized as follows. In Section 2 we introduce
a motivating example and details of the study cohort. We specify the DLM for
our framework in Section 3. The model, which assumes an autoregressive
structure, accounts for possibly multivariate health measures which may or may
not be measured simultaneously. We then develop our computational approach for
inference in Section 4. We apply our methods to modeling BP in Section 5 where
we compare our methods to typical models used to measure the effects of
adherence. We conclude in Section 6 with a discussion on the limitations and
potential extensions of our methods.

##  2 Data

The data we analyzed were obtained from the baseline pre-randomization period
of a trial that studied the effects of a provider-patient communication skill-
building intervention on adherence to anti-hypertensive medication and on BP
control (clinicaltrials.gov ID: NCT00201149). Patients were recruited from
seven outpatient primary care clinics at Boston Medical Center, an inner-city
safety-net hospital affiliated with the Boston University School of Medicine.
Patients enrolled into the study from August 2004 and June 2006, meeting
several eligibility criteria. These included that the patients were of white
or black race (African or Caribbean descent), were at least 21 years old, had
an outpatient diagnosis of hypertension on at least three different occasions
prior to study enrollment, and were currently on anti-hypertensive medication.
The cohort size was 869 patients. The study involved measuring anti-
hypertensive medication-taking using Medication Event Monitoring System (MEMS)
caps, a particular type of electronic pill-top monitoring device. Patients
were given their most frequently taken anti-hypertensive medications in a
bottle with a MEMS cap, and were instructed to open the bottle each time they
took a dose. Each MEMS cap contained a microprocessor that recorded the date
and time whenever the bottle was opened, and the timing information was then
downloaded through a wireless receiver after the patient returned the MEMS
cap. Our study focused on medication-taking behavior during the entire pre-
randomization baseline period of the study. A patient was considered adherent
to the prescribed medication on a day if the MEMS cap was opened as many times
as the daily medication dosing frequency, and not adherent otherwise. These
measurements were recoded on a daily basis, but the duration over which
adherence information was measured varied by patient. Blood pressure
measurements were recorded less frequently as they were obtained as part of
routine clinical care.

The BP readings were assessed using manual or electronic devices, and were
obtained by clinical staff including physicians, nurses and medical
assistants. In cases where multiple readings were obtained on a single day
(typically at the same clinic visit), the individual values were recorded.
Diastolic and systolic blood pressure values (DBP and SBP respectively) were
recorded separately.

In addition to detailed time-varying adherence and BP readings, other patient-
specific baseline information was collected. Race (white versus African
American), gender, income, and age at the start of the study were recorded for
each patient. From electronic medical records, the following comorbidities (as
binary variables) were investigated in our study, given their potential impact
on overall BP levels: Presence of cerebrovascular disease, congestive heart
failure, chronic kidney disease, coronary artery disease, diabetes mellitus,
hyperlipidemia, peripheral vascular disease, and obesity (defined as body mass
index greater than $30kg/m^{2}$).

Figure 1: Example BP and adherence data for a study patient. Red squares are
diastolic BP, red dots are systolic BP, and the black dots near the horizontal
axis indicate whether the patient was adherent on a particular day (filled
dots indicate adherence, open dots indicate non-adherence).

To motivate our modeling framework, consider the data from one of the patients
in the study cohort displayed in Figure 1. The figure shows the DBP and SBP
measures on the three days the patient had their blood pressure measured, and
daily indications of whether the patient was adherent. This patient was
adherent for a total of $29$ of $82$ days. However, this simple summary masks
the time-varying pattern of the patient’s adherence. The patient began the
study by being fully adherent to their prescribed medication, during which
time their blood pressure remained under control. After 20 days, the patient
started becoming less adherent, and starting around day 38 the patient
discontinued taking their medication altogether. Over this latter period, the
patient’s BP increased, and by the end of the study period the patient has DBP
and SBP values that were not under control. Besides suggesting a potential
link between adherence and BP measures, this example motivates a dynamic
linear model for health measures that accounts for time-varying adherence.

##  3 A dynamic model for multivariate health measures

We propose a statistical framework for multivariate time-varying health
measures as a function of medication adherence that is a member of the class
of Bayesian dynamic linear models. Let $y_{itk}$ be the value of the $k^{th}$
health measure, $k=1,\ldots,K$, for patient $i$ at time $t$,
$t=1,\ldots,T_{i}$. We assume that time is discretized and equally spaced
(e.g., days), and that health measures are observed only at times
$t_{1},t_{2},\ldots,t_{m_{i}}$. Let $\bm{y}_{it\bm{\cdot}}$, generally,
indicate a vector collecting all outcomes across the dotted index into a
column vector. In this example we collect the outcomes for patient $i$ at time
$t$, i.e., $\bm{y}_{it\bm{\cdot}}=(\bm{y}_{it1},...,\bm{y}_{itK})^{T}$.
Similarly, $\bm{y}_{i\bm{\cdot}k}$ collects the $k^{th}$ outcome observed at
times $(t_{1},...,t_{m_{i}})$ for patient $i$, i.e.,
$\bm{y}_{i\bm{\cdot}k}=(\bm{y}_{i{t_{1}}k},...,\bm{y}_{i{t_{m}}k})^{T}$. Our
framework assumes that $\bm{y}_{it\bm{\cdot}}$ is an observed measurement
generated from a distribution with mean $\bm{\mu}_{it\bm{\cdot}}$ which
follows a stochastic process. The framework recognizes that actual
measurements on a given day could vary around a mean level due to various
influences including emotional state, activity level, recent alcohol
consumption, ambient temperature, and other unobserved factors that might
affect the outcomes.

The mean process at time $t$ is modeled as the sum of the contributions of a
non-dynamic term and a stochastic term. Specifically, we assume

$\displaystyle\bm{y}_{it\bm{\cdot}}$
$\displaystyle=\bm{\mu}_{it\bm{\cdot}}+\bm{\varepsilon}_{it}$
$\displaystyle=\bm{\beta}\bm{x}_{i}+\bm{\alpha}_{it\bm{\cdot}}+\bm{\varepsilon}_{it}$
(1)

where ${\bm{x}}_{i}$ is a vector of $p$ non-dynamic covariates typically
measured at baseline and $\bm{\beta}$ are the corresponding linear
coefficients, and where $\bm{\alpha}_{it\bm{\cdot}}$ is a stochastic process
that may depend on dynamic covariates, such as time-varying adherence. Here
$\bm{y}_{it\bm{\cdot}},\bm{\mu}_{it\bm{\cdot}},\bm{\alpha}_{it\bm{\cdot}}$ and
$\bm{\varepsilon}_{it}$ are all vectors of length $K$ corresponding to the
observation, mean outcome process, latent state model and sampling error
respectively. The non-dynamic covariate effects $\bm{\beta}$ are of dimension
$K\times p$. We use $\bm{\beta}_{k}$ to denote the $k^{th}$ row of
$\bm{\beta}$, i.e., the non-dynamic covariate effects on the $k^{th}$ outcome.
The error term for subject $i$ at time $t$, $\bm{\varepsilon}_{it}$, accounts
for typical sampling variability, measurement error and possible correlation
of the outcomes within patient and time. We further assume that
$\bm{\varepsilon}_{it}\sim\mbox{N}(0,\;\Sigma_{\varepsilon})$ with possibly
non-diagonal covariance matrix $\Sigma_{\varepsilon}$.

If patient $i$’s adherence is tracked for $T_{i}$ consecutive days, we
indicate this measurement set as $\\{1,...,T_{i}\\}$, and denote the
collection of time-varying adherence measures as $\\{c_{it}\\}_{t=1}^{T_{i}}$.
We let $c_{it}=1$ if patient $i$ was adherent on day $t$ and $-1$ otherwise.
As we describe below, we extend the definition of $c_{it}$ to include other
time-varying covariates. We further let
$\mathcal{T}_{i}=\\{t_{1},...,t_{m_{i}}\\}\subset\\{1,2,...,T_{i}\\}$ denote
the set of $m_{i}$ days on which the outcomes were measured for patient $i$.
The set of outcome measurements for patient $i$ is denoted by
$\bm{y}_{i\bm{\cdot}\bm{\cdot}}=\\{\bm{y}_{it\bm{\cdot}}\\}_{t\in\mathcal{T}_{i}}$.
Crucially, we assume that $m_{i}<<T_{i}$ for all $i$.

Changes in $\bm{\alpha}_{it\bm{\cdot}}$ are assumed to be reflected by whether
a patient takes the prescribed medication along with any other time-varying
covariates. We assume an $AR(1)$ model on the latent states
$\\{\bm{\alpha}_{it\bm{\cdot}}\\}$ and incorporate the $\bm{c}_{it\bm{\cdot}}$
as time-varying covariates. Other forms of stochastic processes are possible,
including autoregressive processes of higher order. More generally, if
$\bm{c}_{it\bm{\cdot}}$ is of dimension $r$, then we assume a process on
$\bm{\alpha}_{it\bm{\cdot}}$ given by

$\displaystyle\bm{\alpha}_{it\bm{\cdot}}$
$\displaystyle=\bm{\rho}\bm{\alpha}_{i,t-1,\bm{\cdot}}+\bm{\phi}\bm{c}_{it\bm{\cdot}}+\bm{\nu}_{it}$
(2)

where $\bm{\rho}$ is the $K$-dimensional diagonal matrix of first-order
autoregressive parameters, that is,
$\bm{\rho}=\text{diag}(\rho_{1},\rho_{2},...,\rho_{K})$. We assume
$|\rho_{k}|<1$ to ensure stationarity. The time-varying outcome effects
$\bm{\phi}$ is a $K\times r$ matrix. We let $\bm{\phi}_{k}$ denote the
$k^{th}$ row of $\bm{\phi}$, i.e., the time-varying covariate effects on the
$k^{th}$ outcome. We further assume in (2) that the innovations
$\bm{\nu}_{it}$ have zero-mean MVN distributions with diagonal covariance
matrix $\Sigma_{\nu}={\text{diag}}(\sigma^{2}_{\nu 1},\sigma^{2}_{\nu
1},...,\sigma^{2}_{\nu K})$. Similarly we assume $\bm{\alpha}_{i1}\sim
N_{K}(0,\Sigma_{0})$ with
$\Sigma_{0}={\text{diag}}(\sigma^{2}_{01},\sigma^{2}_{01},...,\sigma^{2}_{0K})$.
Let $\theta$ denote the non-dynamic parameters
$\theta=(\bm{\beta},\Sigma_{\varepsilon},\bm{\rho},\bm{\phi},\Sigma_{\nu},\Sigma_{0})$.
Equations (3) and (2) define the distributions
$p({\bm{y}}_{it\bm{\cdot}}|\bm{\alpha}_{it\bm{\cdot}},\bm{x}_{i},\theta)$ and
$p(\bm{\alpha}_{it\bm{\cdot}}|\bm{\alpha}_{i,t-1,\bm{\cdot}},\theta,\bm{c}_{it\bm{\cdot}})$
respectively.

Figure 2 contains a simulated example of a patient with a scalar outcome over
a 30-day period corresponding to the model in (3) and (2). The bottom of the
figure displays the adherence indicator simulated independently with a 90%
probability of adherence. The adherence effect is simulated to be large and
negative ($\phi=-0.5$). The contribution of the non-dynamic covariates is
assumed to be $x_{i}^{T}\beta=130$.

As evidenced in Figure 2, the mean process generally declines on days when a
patient is adherent. This is not always the case, and an increase can occur
when the corresponding innovation $\nu_{t}$ is large and positive offsetting
the impact of the patient taking their medication. The observed health measure
is normally distributed around the mean for the day on which the measure is
recorded. On day 1, for example, the health measure is higher than the mean,
and on day 15 the health measure is lower than the mean.

Figure 2:  A 30-day trajectory for a simulated patient from a DLM with a
scalar outcome.

This model is attractive for relating health measures to time-varying
adherence for several reasons. First, the model can account for non-dynamic
baseline variables as well as time-varying adherence. This model feature
permits disentangling the effects of detailed adherence from patient-specific
socio-demographic and health covariates. Second, in settings where multiple
outcome measures are observed at time $t$, inference for the sampling
variability $\Sigma_{\varepsilon}$ can be made more precise, and help to
separate from the innovation variability $\Sigma_{\nu}$.

The assumption of an $AR(1)$ model component has a well-known asymptotic mean
under full adherence and full non-adherence. Specifically, the overall effect
of repeated days of adherence on the outcome measure converges for increasing
values of $t$ as
$\bm{\alpha}_{it\bm{\cdot}}\rightarrow(I-\bm{\rho})^{-1}\bm{\phi}$. We can
therefore predict that if patient $i$ were to continue to be fully adherent,
their mean outcomes would tend to
$\bm{\mu}_{it\bm{\cdot}}\rightarrow\bm{\beta}\bm{x}_{i}+(I-\bm{\rho})^{-1}\bm{\phi}$.
An analogous calculation can be performed when the patient is fully non-
adherent. This property permits estimating the best-case or worst-case health
measure means for perfect adherence or perfect non-adherence, even when
patients’ adherence level is somewhere in between.

##  4 Marginal Dynamic Linear Models

Inference for the model in Section 3 is challenging given the large number of
parameters, both the non-dynamic parameters $\theta$ as well as the time-
varying parameters ${\bm{\alpha}}$. In particular, if the number of recorded
adherence indicators per patient is large, then the size of ${\bm{\alpha}}$ is
similarly large. Standard Bayesian inferential methods through posterior
simulation approaches for such a highly parameterized model can result in slow
convergence and unreliable inferences.

Advances in Bayesian computation have made highly parameterized dynamic models
more tractable. With recent developments in sequential Monte-Carlo (SMC) with
software packages like Libbi (Murray,, 2013), sampling the induced high-
dimensional latent states of the DLM has become computationally feasible and
accessible. The majority of these advances are in situations where the
likelihoods can not be computed directly and are instead approximated.
Marginal sampling schemes, like the particle marginal Metropolis-Hastings
(PMMH) algorithm (Andrieu et al.,, 2010), alternately sample the structural
parameters and latent process parameters and accept or reject them with an
adjusted Metropolis-Hastings step accounting for the approximation of the
likelihood needed in sampling the latent space. Recent work (Bhattacharya and
Wilson,, 2018) approximates the posterior of the structural parameters on a
grid of points, reducing the possible sampling values to a discrete set. SMC
has also been adapted to this situation in methods like SMC2 (Chopin et al.,,
2012), or more aptly named Particle SMC, where one can exploit the known
sequential structure of otherwise intractable likelihoods.

Our approach takes advantage of factoring the joint posterior density of
$\theta$ and $\bm{\alpha}$ as follows

$\displaystyle
p(\bm{\alpha},\theta|\bm{y})=p(\theta|\bm{y})p(\bm{\alpha}|\theta,\bm{y}).$
(3)

where we omit the dependence on both $\\{{\bm{c}}_{it}\\}$ and $\bm{x}_{i}$.
The first factor in (3) is discussed below, while the second factor, the
conditional posterior density of the latent process parameters, can be derived
exactly and is discussed in Section 4.2.

Marginal inference about $\theta$ can be accomplished by integrating (3) with
respect to $\bm{\alpha}$, yielding the first factor in the expression. This
factor can be expanded using Bayes’ Theorem as

$\displaystyle p(\theta|\bm{y})=\frac{p(\bm{y}|\theta)p(\theta)}{p(\bm{y})}.$
(4)

The marginal likelihood, $p(\bm{y}|\theta)$, in the numerator is determined
from

$\displaystyle
p(\bm{y}|\theta)=\int{p(\bm{y},\bm{\alpha}|\theta)}d\bm{\alpha}.$ (5)

As we derive in Section 4.1, the marginal distribution in (5) is multivariate
normal (MVN) with a mean and variance that depend only on the fixed model
parameters $\theta$, the adherence measures $\\{{\bm{c}}_{it}\\}$ and
individual covariates $\bm{x}_{i}$. This marginalization is made possible by
normality of the latent state innovations $\bm{\nu}_{it}$ and sampling error
$\bm{\varepsilon}_{it}$. Inference for $\theta$ can therefore be obtained
directly from this marginal distribution.

###  4.1 Marginal Likelihood of the DLM

We explain our marginalization approach in the case of a DLM with a single
outcome variable ($k=1$) and for one individual, and assume a single baseline
covariate ($p=1$) and single time-varying covariate ($r=1$) which is assumed
to be an adherence indicator. The parameters in $\theta$ are therefore
$(\beta,\sigma^{2}_{\varepsilon},\rho,\phi,\sigma^{2}_{\nu},\sigma^{2}_{0})$.
We assume that both adherence indicators and time-varying outcomes are
measured over $T$ consecutive days; the outcomes are $(y_{1},y_{2},...,y_{T})$
and the adherence indicators are ${\bm{c}}=(c_{1},c_{2},...,c_{{T}})$. We can
write the vector of outcomes as the sum of a shared non-time-varying
component, a time-varying component and an error term as

$\displaystyle\begin{pmatrix}y_{1}\\\ y_{2}\\\ \vdots\\\
y_{T}\end{pmatrix}=\begin{pmatrix}\beta x\\\ \beta x\\\ \vdots\\\ \beta
x\end{pmatrix}+\begin{pmatrix}\alpha_{1}\\\ \alpha_{2}\\\ \vdots\\\
\alpha_{T}\end{pmatrix}+\begin{pmatrix}\varepsilon_{1}\\\ \varepsilon_{2}\\\
\vdots\\\ \varepsilon_{T}\end{pmatrix}=\beta x\
{\mathds{1}}_{T}+\bm{\alpha}+\bm{\varepsilon}.$ (6)

where we use the ${\mathds{1}}_{T}$ notation to indicate a column vector of
length $T$ consisting of all ones. The following development is conditional on
$\theta$ and covariates $(x,\bm{c})$ unless noted otherwise. Thus the first
term on the right hand side of Equation (6), $\beta x\ {\mathds{1}}_{T}$, is
treated as a constant vector. The third term is distributed as

$\displaystyle\bm{\varepsilon}\sim N_{T}({\bm{0}},\sigma_{\varepsilon}^{2}I),$

where ${\bm{0}}$ is a column vector of zeros, and $I$ is a $T$-dimensional
identity matrix.

Because the initial latent variable $\alpha_{1}$ is normally distributed, and
each $\alpha_{t}$ conditional on the previous latent variables is a linear
combination of normal random variables, then $(\alpha_{1},...,\alpha_{T})$ is
MVN.

The AR(1) structure of the latent states admit a recursive mean and variance
calculation

$\displaystyle\mathbb{E}[\alpha_{t}]$
$\displaystyle=\rho\mathbb{E}[\alpha_{t-1}]+\phi\ c_{t},$
$\displaystyle\mathbb{V}ar[\alpha_{t}]$
$\displaystyle=\rho^{2}\mathbb{V}ar[\alpha_{t-1}]+\sigma_{\nu}^{2}.$

This recursion and the initial conditions imply a general formula for the mean
and variance of $\alpha_{t}$

$\displaystyle E_{t}$
$\displaystyle=\mathbb{E}[\alpha_{t}]=\phi\sum_{k=2}^{t}{\rho^{t-k}c_{k}}$ (7)
$\displaystyle V_{t}$
$\displaystyle=\mathbb{V}ar[\alpha_{t}]=\sigma_{\nu}^{2}\sum_{k=2}^{t}{\rho^{2(t-k)}}+\rho^{2(t-1)}\sigma_{0}^{2}$

for $t\in\\{2,...,T\\}$ with initial mean $E_{1}=\mathbb{E}[\alpha_{1}]=0$ and
initial variance $V_{1}=\mathbb{V}ar[\alpha_{1}]=\sigma_{0}^{2}$. We collect
the mean terms into a vector

$\displaystyle{\bf{E}}=(E_{1},...,E_{T}).$ (8)

The covariance follows a similar recursion,

$\displaystyle\mathbb{C}ov(\alpha_{t},\alpha_{t-1})$
$\displaystyle=\mathbb{C}ov(\rho\alpha_{t-1}+\phi\ c_{t}+\nu_{t}\ ,\
\alpha_{t-1})$ $\displaystyle=\rho\ \mathbb{V}ar[\alpha_{t-1}].$

We apply the recursion to obtain a general form of the covariance for
$t\in\\{2,...,T\\}$ and $k\in\\{0,...,t-1\\}$

$\mathbb{C}ov(\alpha_{t},\alpha_{t-k})=\rho^{k}\mathbb{V}ar[\alpha_{t-k}]=\rho^{k}V_{t-k}.$

More compactly, the covariance matrix is

$\displaystyle\Sigma=\begin{pmatrix}\phantom{-}V_{1}&\rho\ V_{1}&\rho^{2}\
V_{1}&\dots&\rho^{T-1}\ V_{1}\\\ \rho\ V_{1}&\phantom{-}V_{2}&\rho\
V_{2}&\dots&\rho^{T-2}\ V_{2}\\\ \rho^{2}\ V_{1}&\rho\
V_{2}&\phantom{-}V_{3}&\dots&\rho^{T-3}\ V_{3}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\ \rho^{T-1}\ V_{1}&\rho^{T-2}\
V_{2}&\rho^{T-3}\ V_{3}&\dots&\phantom{-}V_{T}\end{pmatrix}.$ (9)

The distribution of the terms in Equation (6) can be written explicitly as

$\displaystyle(y_{1},y_{2},...,y_{T})^{T}\sim N_{T}(\beta x\
{\mathds{1}}_{T}+{\bf{E}},\Sigma+\sigma_{\varepsilon}^{2}I).$ (10)

The distribution in (10) is proportional to the marginal likelihood for $T$
sequentially-observed outcomes, marginalizing out the latent parameters.

Let $\mathcal{T}=\\{t_{1},t_{2},...,t_{m}\\}$ be the _actual_ observation
times, and let

$\displaystyle\bm{E}_{\mathcal{T}}=(E_{t_{1}},E_{t_{2}},...,E_{t_{m}})^{T}.$

By properties of the MVN distribution, the outcome vector
$y_{\bm{\cdot}}=(y_{t_{1}},y_{t_{2}},...,y_{t_{m}})^{T}$ is also MVN with mean
and covariance respectively

$\displaystyle\mathbb{E}[y_{\bm{\cdot}}]$ $\displaystyle=$ $\displaystyle\beta
x\ {\mathds{1}}_{|\mathcal{T}|}+\bm{E}_{\mathcal{T}}$
$\displaystyle\mathbb{V}ar[y_{\bm{\cdot}}]$ $\displaystyle=$
$\displaystyle\Sigma_{\mathcal{T}\times\mathcal{T}}+\sigma_{\varepsilon}^{2}I_{|\mathcal{T}|}.$
(11)

The notation $\Sigma_{\mathcal{T}\times\mathcal{T}}$ indicates subsetting the
covariance matrix $\Sigma$ to the corresponding rows and columns designated by
$\mathcal{T}$. If, for example, $\mathcal{T}=\\{1,7,9\\}$ the first, seventh
and ninth rows and columns are taken from $\Sigma$ making
$\Sigma_{\mathcal{T}\times\mathcal{T}}$ a $3\times 3$ matrix. One consequence
of the marginalization is that if we observe multiple measurements of the same
outcome on a single day, they will only vary according to
$\sigma_{\varepsilon}^{2}$.

It is worth noting that even though this marginal distribution does not depend
on the time-varying parameters, the autoregressive structure remains in both
the mean and covariance. For example, the centered marginal mean of $y_{t}$
can be shown with simple algebraic manipulation to be

$\displaystyle\left(\mathbb{E}[y_{t}]-\beta x\right)$
$\displaystyle=\rho\left(\mathbb{E}[y_{t-1}]-\beta x\right)+\phi c_{t}.$

This is the same recursive relationship of the dynamic component of our DLM in
(2). Even though we marginalize out the time-varying parameters, their
sequential structure is retained.

The calculations above were derived for one study participant and a single
outcome variable. These calculations can also be extended to multivariate
outcome measurements. The marginal mean and covariance are both individual-
and outcome- dependent, as both $\bm{\phi}_{k}$ and $\alpha_{itk}$ vary by
outcome and individual. When referring to multiple study participants and
multiple outcome variables, we can make the dependence explicit by denoting
the quantities in (8) and (9) for outcome $k$ of individual $i$ as
${\bf{E}}_{ik}$ and $\Sigma_{ik}$ respectively.

Because the majority of the structure is contained within outcome (recall
$\bm{\rho}$ is a diagonal matrix), we organize the multivariate outcomes as
follows. We collect the outcome $k$ observed through time $t=1,...,T_{i}$ into
a $T_{i}\times 1$ vector
$\bm{y}_{i\bm{\cdot}k}=(y_{i1k},y_{i2k},...,y_{i{T_{i}}k})^{T}$. The marginal
likelihood of a vector of complete outcome measurements can be written as

$\displaystyle\begin{pmatrix}\bm{y}_{i\bm{\cdot}1}\\\ \bm{y}_{i\bm{\cdot}2}\\\
\vdots\\\ \bm{y}_{i\bm{\cdot}k}\\\ \end{pmatrix}\ |\ \theta\sim
N_{kT}\left(\bm{E}_{i}+\bm{\beta}\bm{x}_{i}\otimes{\mathds{1}}_{T_{i}},\;\Sigma_{i}+\Sigma_{\varepsilon}\otimes
I_{T_{i}}\right)$ (12)

where $\otimes$ is the Kronecker product, and where

$\displaystyle{\bm{E}}_{i}=\begin{pmatrix}{\bf{E}}_{i1}\\\ {\bf{E}}_{i2}\\\
\vdots\\\ {\bf{E}}_{ik}\\\
\end{pmatrix}\quad\text{and}\quad\Sigma_{i}=\begin{pmatrix}\Sigma_{i1}&\bm{0}&\cdots&\bm{0}\\\
\bm{0}&\Sigma_{i2}&\cdots&\bm{0}\\\ \vdots&\vdots&\ddots&\vdots\\\
\bm{0}&\bm{0}&\cdots&\Sigma_{ik}\\\ \end{pmatrix}.$ (13)

As in (4.1), we can derive the marginal distribution of the outcome measures
$\bm{y}_{i\bm{\cdot}\bm{\cdot}}|\theta$ by subsetting the appropriate rows and
columns of the quantities in (13).

If we assume that outcomes between study participants are independent
conditional on $\theta$, the marginal distribution of outcomes unconditional
on the time-varying parameters is given by

$\displaystyle
p\left(\bm{y}_{1\bm{\cdot}\bm{\cdot}},...,\bm{y}_{n\bm{\cdot}\bm{\cdot}}|\theta,{\bm{x}},{\bm{c}}\right)=\prod_{i=1}^{n}p\left(\bm{y}_{i\bm{\cdot}\bm{\cdot}}|\theta,x_{i},{\bm{c}}_{i\bm{\cdot}\bm{\cdot}}\right).$
(14)

This simple marginalization drastically reduces the number of parameters in
the model and simplifies computation. We can conduct inference on our non-
dynamic parameters in a Bayesian setting by introducing a prior distribution
for $\theta$ and using the marginal likelihood of $\theta$ which is
proportional to (14). The posterior density for the non-dynamic parameters is
given by

$\displaystyle
p\left(\theta|\left\\{\bm{y}_{i\bm{\cdot}\bm{\cdot}}\right\\}_{i=1}^{n},{\bm{x}},{\bm{c}}\right)\propto
p(\theta)\prod_{i=1}^{n}p\left(\bm{y}_{i\bm{\cdot}\bm{\cdot}}|\theta,x_{i},{\bm{c}}_{i}\right).$
(15)

Inference can be performed in a straightforward and efficient way via MCMC, or
using Hamiltonian Monte Carlo (Homan and Gelman,, 2014) sampling to obtain
draws from the posterior distribution, as implemented in STAN (Carpenter et
al.,, 2017). Alternatively, posterior sampling could be obtained using
Sequential Monte-carlo (SMC) in software packages like Libbi (Murray,, 2013).
SMC works with the sequential likelihoods
$p(y_{it_{j}}|y_{it_{1}},...,y_{it_{j-1}},\theta)$, which are easily available
in our framework because they are conditional distributions of the
Multivariate Normal distribution in (12).

###  4.2 Inference for the Latent Process

Posterior inference for $\bm{\alpha}$ can be determined once
$p(\theta|\bm{y})$ has been obtained by exploiting the factorization of the
posterior density in (3). To do so, we note that the joint distribution of
$\bm{\alpha}$ and $\bm{y}$ conditional on $\theta$ is MVN. For our
development, we consider the case in which outcomes are observed for $T$
consecutive days. With a $k$-dimensional outcome variable, $\bm{\alpha}$
contains $kT$ elements and $\bm{y}$ contains $kT$ (scalar) outcomes. Using the
notation in (13), the joint distribution of outcomes and time-varying
parameters is

$\displaystyle\left.\begin{pmatrix}\bm{\alpha}_{i\bm{\cdot}1}\\\ \vdots\\\
\bm{\alpha}_{i\bm{\cdot}k}\\\ \bm{y}_{i\bm{\cdot}1}\\\ \vdots\\\
\bm{y}_{i\bm{\cdot}k}\\\ \end{pmatrix}\ \right|\ \theta\sim
N_{2kT}\left(\begin{pmatrix}{\bf{E}}_{i}\\\ \\\
{\bf{E}}_{i}+\bm{\beta}\bm{x}_{i}\otimes{\mathds{1}}_{T_{i}}\\\
\end{pmatrix},\begin{pmatrix}\Sigma_{i}&\Sigma_{i}\\\ \\\
\Sigma_{i}&\Sigma_{i}+\Sigma_{\varepsilon}\otimes I_{T_{i}}\\\
\end{pmatrix}\right).$ (16)

The quantities on the diagonal of the covariance matrix in (16) were
determined in Section 4.1. The covariance of $\bm{y}_{i\bm{\cdot}k}$ and
$\bm{\alpha}_{i\bm{\cdot}\ell}$ (for $k\neq\ell$), conditional on $\theta$,
can be shown to be $\bm{0}$ as follows (with conditioning on $\theta$
suppressed).

$\displaystyle\mathbb{C}ov(\bm{\alpha}_{i\bm{\cdot}l},\bm{y}_{i\bm{\cdot}k})$
$\displaystyle=\mathbb{E}\left[\mathbb{C}ov\left(\bm{\alpha}_{i\bm{\cdot}l},\bm{y}_{i\bm{\cdot}k}|\bm{\alpha}_{i\bm{\cdot}k}\right)\right]+\mathbb{C}ov\left(\mathbb{E}[\bm{\alpha}_{i\bm{\cdot}l}|\bm{\alpha}_{i\bm{\cdot}k}],\mathbb{E}[\bm{y}_{i\bm{\cdot}k}|\bm{\alpha}_{i\bm{\cdot}k}]\right)$
$\displaystyle=\mathbb{E}\left[\mathbb{C}ov\left(\bm{\alpha}_{i\bm{\cdot}l},\beta
x_{i}+\bm{\alpha}_{i\bm{\cdot}k}+\bm{\varepsilon}_{i\bm{\cdot}k}|\bm{\alpha}_{i\bm{\cdot}k}\right)\right]+\mathbb{C}ov\left(\bm{\alpha}_{i\bm{\cdot}l},\beta
x_{i}+\bm{\alpha}_{i\bm{\cdot}k}\right)$
$\displaystyle=\mathbb{E}\left[\mathbb{C}ov\left(\bm{\alpha}_{i\bm{\cdot}l},\bm{\varepsilon}_{i\bm{\cdot}k}\right)\right]+\mathbb{C}ov\left(\bm{\alpha}_{i\bm{\cdot}l},\bm{\alpha}_{i\bm{\cdot}k}\right)$
$\displaystyle=\bm{0}.$

We can subset the mean vector and covariance matrix of (16) as in (4.1) to
account for the sporadically observed outcomes, but we continue to assume that
the response values are observed at each time period. The conditional
posterior distribution of the latent parameters is determined from the joint
distribution in (16) by conditioning on $\bm{y}$ as

$\displaystyle\left.\begin{pmatrix}\bm{\alpha}_{i\bm{\cdot}1}\\\ \vdots\\\
\bm{\alpha}_{i\bm{\cdot}k}\\\
\end{pmatrix}\right|\bm{y}_{i\bm{\cdot}1},...,\bm{y}_{i\bm{\cdot}K},\theta\sim
N_{kT}\left(\tilde{\alpha}_{i},\tilde{\Sigma}_{i}\right)$ (17)

where

$\displaystyle\tilde{\alpha}_{i}=\left({\bf{E}}_{i}\right)+\Sigma_{i}\left(\Sigma_{i}+\Sigma_{\varepsilon}\otimes
I_{T_{i}}\right)^{-1}\left[\begin{pmatrix}\bm{y}_{i\bm{\cdot}1}\\\ \vdots\\\
\bm{y}_{i\bm{\cdot}k}\\\
\end{pmatrix}-\left({\bf{E}}_{i}+\bm{\beta}\bm{x}_{i}\otimes{\mathds{1}}_{T_{i}}\right)\right]$

and

$\displaystyle\tilde{\Sigma}_{i}=\Sigma_{i}-\Sigma_{i}\left(\Sigma_{i}+\Sigma_{\varepsilon}\otimes
I_{T_{i}}\right)^{-1}\Sigma_{i}.$

We can use the joint distribution specified in (17) along with the posterior
distribution on the structural parameters shown in (15) to sample from the
latent process parameters by first simulating $\theta^{*}\sim
p\left(\theta|\bm{y}_{1\bm{\cdot}\bm{\cdot}},...,\bm{y}_{n\bm{\cdot}\bm{\cdot}},{\bm{x}},{\bm{c}}\right)$
and then sampling from
$p(\bm{\alpha}_{i\bm{\cdot}1},...,\bm{\alpha}_{i\bm{\cdot}1}|\bm{y}_{1\bm{\cdot}\bm{\cdot}},...,\bm{y}_{n\bm{\cdot}\bm{\cdot}},\theta^{*})$.

##  5 Application to BP Study

We return to an application of our framework to modeling time-varying BP as a
function of time-varying adherence to anti-hypertensive medication and
baseline comorbidities and demographics described in Section 2. We begin with
a discussion of models often used in this task that incorporate adherence as
non-dynamic information. We also address missing adherence measures for a
fraction of patients.

###  5.1 Non-dynamic models incorporating adherence

Adherence to medication is typically incorporated into outcome models in one
of two ways. The average adherence for the study period is either used
directly or is dichotomized into two groups, those below a certain threshold
and those above, indicating “poor” versus “good” adherence (Rose et al.,,
2011; Lee et al.,, 2006; Schroeder et al.,, 2004). Either of these approaches
includes adherence as a non-dynamic covariate in the model. Repeated outcome
measures are accompanied by patient-specific random effects.

We present these alternative models in the context of BP outcomes, a bivariate
measure. We label the outcomes $1$ for Systolic BP and $2$ for Diastolic BP.
These models take the form

$\displaystyle\begin{pmatrix}y_{it1}\\\
y_{it2}\end{pmatrix}=\begin{pmatrix}x_{i}\beta_{1}\\\
x_{i}\beta_{2}\end{pmatrix}+\begin{pmatrix}\delta_{i1}\\\
\delta_{i2}\end{pmatrix}+\begin{pmatrix}\bar{c}_{i}\gamma_{1}\\\
\bar{c}_{i}\gamma_{2}\end{pmatrix}+\begin{pmatrix}\varepsilon_{it1}\\\
\varepsilon_{it2}\end{pmatrix}$ (18)

with two different adherence measures, $\bar{c}_{i}$. We first consider
average adherence, $\bar{c}_{i}=T_{i}^{-1}\sum_{t=1}^{T_{i}}c_{it}$ over
$T_{i}$ days for patient $i$. In this case we would interpret the adherence
effect parameters $(\gamma_{1},\gamma_{2})$ as the differences in BP of being
fully adherent relative to being fully non-adherent, controlling for the
baseline covariates. We also consider a dichotomized summary
$\bar{c}_{i}=\mathds{1}\left\\{\left(T_{i}^{-1}\sum_{t=1}^{T_{i}}c_{it}\right)>2p-1\right\\}$,
an indicator of overall adherence. In this case $(\gamma_{1},\gamma_{2})$
would be interpreted as the differences in BP for those with “good”
($\bar{c}_{i}=1$) versus “poor” ($\bar{c}_{i}=0$) adherence. Several values of
$p$ were considered to assess the sensitivity to this choice, but we present
results for $p=0.8$ which is a conventional choice (Schroeder et al.,, 2004).
The model in Equation (18) includes patient-specific random effects
$\delta_{i1}\sim N(0,\sigma_{\delta 1}^{2})$ and $\delta_{i2}\sim
N(0,\sigma_{\delta 2}^{2})$.

We compare the fit of the above two non-dynamic models to our dynamic model
framework. The non-dynamic approach has a potential advantage of being more
robust to detailed model misspecification relative to our dynamic model,
particularly with the choice of the specific model for the evolution of the
health measures as a function of daily adherence. However, incorporating
average adherence may mask important time-varying effects of detailed
adherence. We explore this tradeoff in Section 5.3.

###  5.2 Missing adherence indicators

Not all of the patients in our cohort have completely observed adherence
indicators. Among the $503$ patients in our analyses, $70$ have at least one
day of missing adherence. Adherence could be missing due to MEMS cap
malfunctions, hospital inpatient stays in which the MEMS containers were not
used, or other causes. The majority ($50$) of the $70$ patients had only one
or two missing adherence values. In the most extreme case one patient was
missing $48$ out of $102$ adherence measures.

For patient $i$, let $\bm{c}_{i}^{obs}=\\{c_{it}\\}_{t\in\mathcal{T}_{obs}}$
be the set of observed adherence values and
$\bm{c}_{i}^{mis}=\\{c_{it}\\}_{t\in\mathcal{T}_{mis}}$ be the set of missing
adherence values. Letting $\eta_{i}$ represent the parameters describing a
potential adherence model for patient $i$, the posterior density in (15)
conditional on the observed adherence data only is

$\displaystyle
p\big{(}\theta|{\bm{y}}_{i\bm{\cdot}\bm{\cdot}},x_{i},\bm{c}_{i}^{obs}\big{)}$
$\displaystyle=\int{p\big{(}\theta,\bm{c}_{i}^{mis},\eta_{i}|{\bm{y}}_{i\bm{\cdot}\bm{\cdot}},x_{i},\bm{c}_{i}^{obs}\big{)}}\
d\bm{c}_{i}^{mis}\ d\eta_{i}$
$\displaystyle=\int{p\big{(}\theta|{\bm{y}}_{i\bm{\cdot}\bm{\cdot}},x_{i},\bm{c}_{i}^{obs},\bm{c}_{i}^{mis},\eta_{i}\big{)}p\big{(}\bm{c}_{i}^{mis}|{\bm{y}}_{i\bm{\cdot}\bm{\cdot}},x_{i},\bm{c}_{i}^{obs},\eta_{i}\big{)}p(\eta_{i}|{\bm{y}}_{i\bm{\cdot}\bm{\cdot}},x_{i},\bm{c}_{i}^{obs})}\
d\bm{c}_{i}^{mis}\ d\eta_{i}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\int{p\big{(}\theta|{\bm{y}}_{i\bm{\cdot}\bm{\cdot}},x_{i},\bm{c}_{i}^{obs},\bm{c}_{i}^{mis},\eta_{i}\big{)}p\big{(}\bm{c}_{i}^{mis}|\eta_{i}\big{)}p(\eta_{i}|\bm{c}_{i}^{obs})}\
d\bm{c}_{i}^{mis}\ d\eta_{i}.$ (19)

The last line in (5.2) labeled $(a)$ involves a set of modeling choices.
First, we assume that the adherence model parameters do not depend on other
covariates or the outcomes, and this is reflected by assuming
$p(\eta_{i}|{\bm{y}}_{i\bm{\cdot}\bm{\cdot}},x_{i},\bm{c}_{i}^{obs})=p(\eta_{i}|\bm{c}_{i}^{obs})$.
Other approaches, such as those described by Naranjo et al., (2013), provide a
framework for modeling missing time-varying covariates in DLMs with complex
models. However, our assumption is conservative but likely more robust to
model misspecification. Second,
$p\big{(}\bm{c}_{i}^{mis}|{\bm{y}}_{i\bm{\cdot}\bm{\cdot}},x_{i},\bm{c}_{i}^{obs},\eta_{i}\big{)}=p\big{(}\bm{c}_{i}^{mis}|\eta_{i}\big{)}$,
implies that the missing adherence values can be simulated directly from our
adherence model, represented by $\eta_{i}$, without regard to individual-level
characteristics. Again, this is a conservative choice because adding
information to the imputation model could help improve predictions if properly
modeled.

We simulate the missing adherence measures from Beta-Bernoulli distributions
that depend only on adherence measures for each patient separately. Letting
$\eta_{i}$ represent patient $i$’s overall adherence, the unobserved adherence
values are drawn from a Bernoulli distribution

$\\{c_{it}^{*}\\}_{t\in\mathcal{T}}|\eta_{i}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\text{Bernoulli}(\eta_{i})$

where the $c_{it}=2c_{it}^{*}-1$ take on values $\\{-1,1\\}$. Under a uniform
prior distribution on patient $i$’s adherence rate, the posterior distribution
for $\eta_{i}$ is given by

$\eta_{i}|\\{c_{it}\\}_{t\in\mathcal{T}_{obs}}\sim\mbox{Beta}(n_{i,1}+1,n_{i,-1}+1)$

when they were observed to be adherent $n_{i,1}$ days and non-adherent
$n_{i,-1}$ days.

We repeated the simulation process $20$ times and combined the posterior
samples for a proper Bayesian multiple-imputation analysis (Little and Rubin,,
2002) for our non-dynamic models. For our fully Bayesian model, the simulated
adherence values were incorporated into our posterior simulation analyses. We
found that in both cases inference of our non-dynamic parameters $\theta$ was
not sensitive to the missing adherence values.

###  5.3 Analysis of BP measures

Table 1 displays the number of BP measurements among our $503$ patients during
the study period. While a majority ($69.8\%$) had either one or two BP
readings, a substantial number of patients had $3$ or more. The number of BP
readings totalled $1152$, averaging $2.29$ readings per patient.

Number of BP Readings | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $7$ | $8+$
---|---|---|---|---|---|---|---|---
Number of Patients | $226$ | $125$ | $72$ | $28$ | $16$ | $12$ | $15$ | $9$
Percent of Patients | $44.90$ | $24.90$ | $14.30$ | $5.60$ | $3.20$ | $2.40$ | $3.00$ | $1.79$
Table 1:  Number of BP readings during the study period.

For the period of the study, the proportion of days that patients took their
medication varied widely. The proportions ranged from $8.9\%$ to $100\%$, with
a median of $95.1\%$ and a mean of $88.9\%$. These adherence rates on average
were high, with many patients being fully adherent throughout the study
period. Only $19.1\%$ of patients had below $80\%$ adherence. The high degree
of adherence is consistent with recruiting patients for the study who were
continual users of anti-hypertensive medication. On average, adherence was
recorded over $98$ days per patient (minimum and maximum of $21$ and $395$
days, respectively, with an inter-quartile range of $14$), with $75\%$ of
patients followed between $84$ and $112$ days.

Baseline summaries of the non-dynamic covariates appear in Table 2.

| Mean (Std Dev)
---|---
Age (y) | $60.3(11.1)$
DBP at Enrollment | $79.7(11.5)$
SBP at Enrollment | $133.5(19.6)$
| Percent
Female | 67.5
African-American | 54.9
Income below $\$20,000$ | 44.3
Obese | 60.6
Cerebral Vascular Disease | 5.6
Congestive Heart Failure | 3.8
Renal Insufficiency | 6.8
Coronary Artery Disease | 14.9
Diabetes | 36.2
Hyperlypidemia | 57.1
Peripheral Vascular Disease | 6.2
Table 2: Baseline socio-demographic and health characteristics of the patients
in the study cohort.

A majority of the cohort consisted of women, more patients of black race
(African or Caribbean descent) than white, and a large fraction of low-income
patients. The cohort also consisted of mostly obese patients, and had a
moderately high comorbidity burden. Based on 503 patients with BP readings
within 14 days of enrollment, the cohort on average had relatively well-
controlled hypertension at baseline, as all patients were prescribed anti-
hypertensive medication (though adherent to varying extents). The cohort
consisted of $26.2\%$ having DBP$>$80 mm Hg and SBP$>$130 mm Hg), and $7.8\%$
having DBP$>$90 mm Hg and SBP$>$140 mm Hg).

A prior distribution was assumed for the alternative models specified in (18).
For $k=1,2$, indicating systolic and diastolic BP measures separately, we
assumed the following.

$\displaystyle\rho_{k}$ $\displaystyle\sim U(-1,1)$ $\displaystyle\phi_{k}$
$\displaystyle\sim N(0,25)$ $\displaystyle\gamma_{k}$ $\displaystyle\sim
N(0,25)$ $\displaystyle\sigma_{\varepsilon k}$ $\displaystyle\sim
U(0,30),\rho_{\varepsilon}\sim U(-1,1)$ $\displaystyle\sigma_{\nu k}$
$\displaystyle\sim U(0,10),\sigma_{0k}\sim U(0,30)$ $\displaystyle\beta_{11}$
$\displaystyle\sim N(120,400),\beta_{12}\sim N(80,400)$
$\displaystyle\beta_{jk}$ $\displaystyle\sim N(0,400),j=2,...,p$
$\displaystyle\sigma_{Uk}$ $\displaystyle\sim U(0,30)$

The prior components were selected to be vague but proper. The intercepts
$\beta_{1k}$ had distributions centered near the typical systolic and
diastolic BPs, but had variances that were sufficiently large to acknowledge
the uncertainty in the effects. We assumed uniform prior components with
compact support for the standard deviation parameters, as recommended by
Gelman et al., (2006). The correlation and autocorrelation parameters were
assumed to have uniform prior components as in the dynamic model
parametrization. Convergence of the MCMC simulated values was inspected with
trace plots of multiple chains, and using the Gelman-Rubin convergence
statistic (Gelman and Rubin,, 1992).

Table 3 presents posterior means and 90% central posterior intervals for the
non-dynamic covariate effects using the DLM.

Variable | Systolic | ($90\%$ CI) | Diastolic | ($90\%$ CI)
---|---|---|---|---
Intercept | $132.97$ | $(129.46,136.51)^{\star\dagger\ddagger}$ | $85.68$ | $(83.65,87.72)^{\star\dagger\ddagger}$
Sex (male) | $-1.1$ | $(-3.59,1.38)$ | $0.53$ | $(-0.88,1.94)$
Age (group 1) | $-0.01$ | $(-3.3,3.27)$ | $-1.5$ | $(-3.4,0.4)$
Age (group 2) | $1.87$ | $(-1.45,5.19)$ | $-3.54$ | $(-5.45,-1.63)^{\star\dagger\ddagger}$
Age (group 3) | $5.22$ | $(1.57,8.89)^{\star\dagger\ddagger}$ | $-6.73$ | $(-8.83,-4.62)^{\star\dagger\ddagger}$
White | $-3.31$ | $(-5.61,-1.01)^{\star\dagger\ddagger}$ | $-1.58$ | $(-2.88,-0.26)^{\star}$
Obese | $3.07$ | $(0.76,5.38)^{\star\dagger\ddagger}$ | $1.68$ | $(0.35,3.01)^{\star\dagger\ddagger}$
Nicotine dependence | $-0.9$ | $(-5.22,3.42)$ | $1.3$ | $(-1.17,3.76)$
Hyperlipidemia | $-1.44$ | $(-3.71,0.83)$ | $-1.34$ | $(-2.64,-0.04)^{\star}$
Diabetes | $1.44$ | $(-0.96,3.85)$ | $-2.98$ | $(-4.34,-1.61)^{\star\dagger\ddagger}$
Peripheral vascular disease | $-1.28$ | $(-5.75,3.17)$ | $-2.6$ | $(-5.13,-0.06)^{\star\dagger}$
Renal insufficiency | $-0.67$ | $(-4.88,3.52)$ | $-2.73$ | $(-5.14,-0.33)^{\star\dagger\ddagger}$
Benign prostatic hypertrophy | $2.82$ | $(-3.66,9.3)$ | $-1.59$ | $(-5.34,2.15)$
Coronary artery disease | $-1.66$ | $(-4.77,1.45)$ | $-2.76$ | $(-4.54,-0.98)^{\star\dagger\ddagger}$
Congestive heart failure | $-0.78$ | $(-6.05,4.46)$ | $0.59$ | $(-2.44,3.59)$
Cerebral vascular disease | $2.45$ | $(-2.19,7.06)$ | $-0.17$ | $(-2.82,2.48)$
Table 3: Summaries of covariate effects for the bivariate dynamic linear
model. ⋆Effect with 90% posterior interval not containing 0 in DLM,
†Significant effect at the 0.1 level in average adherence model, ‡Significant
effect at the 0.1 level in dichotomized adherence model.

The point estimates and intervals are reported for the DLM only; the effects
in the alternative models that were significant at the $0.1$ level are
indicated with a $\dagger$ (for the average adherence model) or a $\ddagger$
(for the dichotomized adherence model). Effects with 90% central posterior
intervals not containing 0 marked with an asterisk ($\star$). Based on the
model fits, the estimated covariate effects tend to be similar across all
models with the point estimates tending to agree in magnitude and sign. Even
though the DLM covariate effects tended to have narrower intervals on average
($3\%$ reduction), the significance of the findings tended to agree as well.
In particular, the effect of race (white versus non-white) was significantly
negative, indicating that whites tended to have lower blood pressure
controlling for all other variables and time-varying adherence. Patients who
were obese at the beginning of the study tended to have significantly higher
DBP and SBP. These findings are consistent with the results of previous
studies (Kressin et al.,, 2010; Rose et al.,, 2011) in their significance and
direction of the effects. Both of these, except for the effect of being white
on mean diastolic BP, agreed with the alternative models in terms of
significance and directionality of the effect.

Table 4 reports inferences for the standard error and correlations of
$(\bm{\varepsilon}_{it})$ for all three models.

Dynamic linear model
---
Variable | Systolic | ($90\%$ CI) | Diastolic | ($90\%$ CI)
Standard Errors: | $13.37$ | $(12.46,14.25)$ | $8.11$ | $(7.68,8.54)$
Correlation: | $0.62$ | $(0.57,0.67)$ |  |
Average Adherence Model
Variable | Systolic | ($90\%$ CI) | Diastolic | ($90\%$ CI)
Standard Errors: | 14.53 | (13.87, 15.21) | 8.34 | (7.97, 8.72)
Correlation: | 0.58 | (0.54, 0.62) |  |
Dichotomized Adherence Model
Variable | Systolic | ($90\%$ CI) | Diastolic | ($90\%$ CI)
Standard Errors: | 14.5 | (13.86, 15.18) | 8.34 | (7.96, 8.72)
Correlation: | 0.58 | (0.54, 0.62) |  |
Table 4: Measurement error estimates for three models.

Comparing the sampling standard deviation estimates across models provides an
indication of the gains in modeling the adherence effects as time-varying. The
standard error estimates for the alternative adherence models tend to be
slightly larger than those given by the DLM, which is consistent with previous
work (Rose et al.,, 2011). The time-varying adherence explicitly captured in
the DLM may account for the extra variation in the outcomes of the alternative
adherence models through a reduction in the estimated measurement error
variance.

Table 5 contains the adherence effects for our models.

State Space Model
---
Variable | Systolic | ($90\%$ CI) | Diastolic | ($90\%$ CI)
Adherence effect: | -0.48 | (-0.84, -0.2)⋆ | -0.24 | (-0.43, -0.09)⋆
Asymptotic Adherence effect: | -3.87 | (-5.98, -1.83)⋆ | -3.15 | (-4.38, -1.94)⋆
Alternative Adherence Model: Average Adherence
Variable | Systolic | ($90\%$ CI) | Diastolic | ($90\%$ CI)
Adherence effect: | -9.24 | (-16.03, -2.22)⋆ | -9.46 | (-13.55, -5.52)⋆
Alternative Adherence Model: Dichotomized Adherence ($p=0.8$)
Variable | Systolic | ($90\%$ CI) | Diastolic | ($90\%$ CI)
Adherence effect: | -5.49 | (-8.26, -2.75)⋆ | -3.78 | (-5.27, -2.2)⋆
Table 5: Adherence effects for the models considered

The adherence effects across the three models are not directly comparable,
given the different approaches to incorporate adherence. The average adherence
model indicates that the difference between those who were fully adherent and
those who were fully non-adherent, controlling for other covariates, is about
$-9.2$ and $-9.5$ for systolic and diastolic, respectively. This implies, for
example, that a 10% additive increase in adherence corresponds to a $0.92$ and
$0.95$ reduction in systolic and diastolic blood pressure, respectively.
However, this effect size assumes that the relationship between average
adherence and blood pressure is linear and holds throughout the entire range
of adherence. Given the limited range of average adherence observed in the
data ($90\%$ of patients have average adherence above $70.5\%$), interpreting
this effect beyond this range is not recommended because it involves
extrapolating beyond the data. The dichotomized adherence model shows similar
results. In particular this approach involves comparing those with relatively
good overall adherence (above 80% adherent) to everyone else. Based on the
dichotomized adherence model, the benefit of being in the former group is
indicated by a lower blood pressure of $-5.49$ and $-3.78$ on average for
systolic and diastolic blood pressure.

The DLM gives similar results. The effect of taking medication can be inferred
on a daily basis. The results of the model fit suggest a small but significant
reduction of blood-pressure from taking the medication on a daily basis,
$-0.24$ mm Hg and $-0.12$ mm Hg on average for systolic and diastolic BP,
respectively. These estimates imply that, accounting for the correlation
estimate, a patient who is adherent over consecutive days would experience a
long-term reduction in systolic BP by $-3.9$ mm Hg, and a long-term reduction
in diastolic BP by $-3.15$ mm Hg. The magnitudes of an increase in long-term
BP for continued non-adherence is the same but in the opposite direction.
Overall, the adherence effects tend to agree in terms of the significance and
direction for the different models. However, the DLM provides a clearer
interpretation of these parameters that is consistent with the time-varying
nature of the data.

Another benefit of the DLM is our ability to infer the latent BP for
unobserved days using the procedure discussed in Section 4.2. Figure 3
contains an example of posterior draws of the mean process
$\bm{\mu}_{it\bm{\cdot}}=\bm{\beta}\bm{x}_{i}+\bm{\alpha}_{it\bm{\cdot}}$ for
the patient presented in Figure 1.

Figure 3: 90% central posterior intervals for the mean blood pressure
$\\{\mu_{it\bm{\cdot}}\\}$ for one study subject. The solid lines are the
estimated mean DBP and SBP, and the dashed lines are the point-wise posterior
intervals. Horizontal dotted lines are drawn at the baseline-estimated mean
DBP and SBP.

The solid line indicates the posterior mean of the latent process while the
dashed lines indicate a $90\%$ credible envelope across time. That is, for
each time $t$, a $90\%$ credible interval of
$\bm{\alpha}_{it1},\bm{\alpha}_{it2}|\left\\{\bm{y}_{i\bm{\cdot}\bm{\cdot}}\right\\}_{i=1}^{n}$
is shown. The estimated mean DBP and SBP processes vary gradually over time,
with shifts in direction indicated by medication adherence variation.

##  6 Discussion

In this article we propose a multivariate DLM that can be used for modeling
time-varying outcome measures as a function of detailed medication adherence
or other time-varying covariates. While DLMs are of common use, our particular
setting benefits from a two-stage computational approach made possible by
factorizing the posterior density into the dynamic and non-dynamic parameters.
Typical analyses in this setting ignore the time-varying structure of
medication adherence and instead examine this relationship between adherence
and outcomes via correlations with non-dynamic adherence measures (e.g., time-
averaged adherence). Our framework explicitly provides a measure of daily
impact of medication-taking and meaningful bounds for mean BP while
controlling for baseline covariates.

Our modeling approach is sufficiently flexible to permit a wide range of
assumptions distinct from those we included in our hypertension application.
For example, we assumed an AR(1) modeling structure for the mean outcome
process, and this assumption can be justified based on the pharmacokinetic
properties of the medications. In other settings, alternative mean processes
can be considered, including higher order AR processes, growth curves, and so
on. Our baseline covariates were modeled linearly, but our framework permits
non-linear inclusion of covariate information for both baseline non-dynamic
covariates, as well as the dynamic predictors in the mean process component of
our model. A crucial assumption of our framework is that the outcome
distribution is multivariate normal, as this assumption allows for the
marginalization strategy that leads to an efficient computational procedure.
However, with some non-normal outcome distributions, various strategies can be
employed that can potentially take advantage of the marginalization idea. For
example, non-normal outcome densities can be approximated by normal
distributions after which the marginalization can be performed; then a
Metropolis-Hastings algorithm would be incorporated into the posterior
sampling procedure that would account for the non-normality of the original
sampling distribution. Such a procedure is likely to be far more efficient
than sampling the dynamic parameters directly as part of model fitting.

The proposed modeling framework also permits extensions that are
straightforward to incorporate. For example, multivariate outcomes that are
recorded at staggered intervals can be accommodated by marginalizing the
distribution in (12) appropriately. We may also believe that the effect of
adherence to medication varies from person to person. We can then include
patient-specific adherence effects with a hierarchical prior on these
parameters to share information across individuals. Our analyses did not
acknowledge differences among anti-hypertensive medications, but our framework
permits distinguishing differential medication effects in a natural manner.
The effects of different medications, perhaps grouped by relevant
characteristics such as whether the medication is short-acting or long-
lasting, or by medication type (e.g., for hypertension, diuretics, ACE
inhibitors, etc.), can be included as separate time-varying effects in the
dynamic component of our model. Dosages and dosing frequencies can serve as
covariates for the effects of particular medications.

The framework we developed can help establish answers to questions of interest
to clinicians and medical researchers that have been difficult to assess
through simpler models. In particular, our model can establish the effects of
different socio-demographic or health factors on health outcomes that control
for detailed time-varying adherence to medication by recognizing that
medication-taking behaviors can change over time. From our framework, we can
also estimate the daily improvement in being adherent to one’s medication,
controlling for socio-demographic and health characteristics, but also the
likely long-range achievable mean outcomes. We can also use the model to
forecast health outcomes as a function of specified patterns of adherence,
potentially serving as a tool for medical decision-making by clinicians. Our
approach provides a robust framework for understanding the impacts of poor
medication adherence as clinicians and patients work together to improve their
medication treatment.

## References

  * Amery et al., (1985) Amery, A., Brixko, P., Clement, D., De Schaepdryver, A., Fagard, R., Forte, J., Henry, J. F., Leonetti, G., O’Malley, K., and Strasser, T. (1985).  Mortality and morbidity results from the European Working Party on High Blood Pressure in the Elderly trial.  The Lancet, 325(8442):1349–1354.
  * Andrieu et al., (2010) Andrieu, C., Doucet, A., and Holenstein, R. (2010).  Particle Markov chain Monte Carlo methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72(3):269–342.
  * Benjamin et al., (2018) Benjamin, E. J., Virani, S. S., Callaway, C. W., Chamberlain, A. M., Chang, A. R., Cheng, S., Chiuve, S. E., Cushman, M., Delling, F. N., Deo, R., et al. (2018).  Heart disease and stroke statistics—2018 update: a report from the american heart association.  Circulation, 137(12):e67–e492.
  * Bhattacharya and Wilson, (2018) Bhattacharya, A. and Wilson, S. P. (2018).  Sequential bayesian inference for static parameters in dynamic state space models.  Computational Statistics & Data Analysis, 127:187 – 203.
  * Carpenter et al., (2017) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017).  Stan: A Probabilistic Programming Language.  Journal of Statistical Software, 76(1):1–32.
  * Choo et al., (1999) Choo, P. W., Rand, C. S., Inui, T. S., Lee, M.-L. T., Cain, E., Cordeiro-Breault, M., Canning, C., and Platt, R. (1999).  Validation of patient reports, automated pharmacy records, and pill counts with electronic monitoring of adherence to antihypertensive therapy.  Medical care, 37(9):846–857.
  * Chopin et al., (2012) Chopin, N., Jacob, P. E., and Papaspiliopoulos, O. (2012).  SMC 2: an efficient algorithm for sequential analysis of state space models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 75(3):397–426.
  * Durbin and Koopman, (2001) Durbin, J. and Koopman, S. J. (2001).  Time series analysis by state space methods.  Oxford University Press, Oxford; New York.
  * Gelman et al., (2006) Gelman, A. et al. (2006).  Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper).  Bayesian analysis, 1(3):515–534.
  * Gelman and Rubin, (1992) Gelman, A. and Rubin, D. B. (1992).  Inference from iterative simulation using multiple sequences.  Statist. Sci., 7(4):457–472.
  * Homan and Gelman, (2014) Homan, M. D. and Gelman, A. (2014).  The No-U-turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.  J. Mach. Learn. Res., 15(1):1593–1623.
  * Kressin et al., (2010) Kressin, N. R., Orner, M. B., Manze, M., Glickman, M. E., and Berlowitz, D. (2010).  Understanding contributors to racial disparities in blood pressure control.  Circulation: Cardiovascular Quality and Outcomes, 3(2):173–180.
  * Lee et al., (2006) Lee, J. K., Grace, K. A., and Taylor, A. J. (2006).  Effect of a pharmacy care program on medication adherence and persistence, blood pressure, and low-density lipoprotein cholesterol: A randomized controlled trial.  JAMA, 296(21):2563–2571.
  * Little and Rubin, (2002) Little, R. J. and Rubin, D. B. (2002).  Statistical analysis with missing data. John Wiley & Sons.  New York.
  * Morisky et al., (1986) Morisky, D. E., Green, L. W., and Levine, D. M. (1986).  Concurrent and predictive validity of a self-reported measure of medication adherence.  Medical care, 24(1):67–74.
  * Murray, (2013) Murray, L. M. (2013).  Bayesian State-Space Modelling on High-Performance Hardware Using LibBi.  ArXiv e-prints.
  * Naranjo et al., (2013) Naranjo, A., Trindade, A. A., and Casella, G. (2013).  Extending the State-Space Model to Accommodate Missing Values in Responses and Covariates.  Journal of the American Statistical Association, 108(501):202–216.
  * Osterberg and Blaschke, (2005) Osterberg, L. and Blaschke, T. (2005).  Adherence to medication.  New England Journal of Medicine, 353(5):487–497.
  * Petris et al., (2009) Petris, G., Petrone, S., and Campagnoli, P. (2009).  Dynamic linear models.  In Dynamic Linear Models with R, pages 31–84. Springer, New York, NY, New York, NY.
  * Probstfield, (1991) Probstfield, J. L. (1991).  Prevention of stroke by antihypertensive drug-treatment in older persons with isolated systolic hypertension-final results of the Systolic Hypertension in the Elderly Program (SHEP).  Journal of the American Medical Association, 265(24):3255–3264.
  * Rose et al., (2011) Rose, A. J., Glickman, M. E., D’Amore, M. M., Orner, M. B., Berlowitz, D., and Kressin, N. R. (2011).  Effects of daily adherence to antihypertensive medication on blood pressure control.  The Journal of Clinical Hypertension, 13(6):416–421.
  * Schroeder et al., (2004) Schroeder, K., Fahey, T., and Ebrahim, S. (2004).  How can we improve adherence to blood pressure–lowering medication in ambulatory care? Systematic review of randomized controlled trials.  Archives of Internal Medicine, 164(7):722–732.
  * West and Harrison, (1997) West, M. and Harrison, J. (1997).  Bayesian forecasting and dynamic models (2nd ed.).  Springer-Verlag New York, Inc., New York, NY, USA.