Analysis of N-of-1 trials using Bayesian distributed lag model with autocorrelated errors

Suppose we observe data from a patient on $n$ consecutive days. On day $t=1,\ldots,n$, let $X_{t}$ and $Y_{t}$ denote the binary treatment indicator and the outcome of interest, respectively. We consider a distributed lag autoregressive model for $Y$, described as follows:

for $t=p+1,...,n$, where the error term $\epsilon_{t}$ follows an autoregressive process,

$w_{t}$ is a white Gaussian noise with mean zero and unknown variance $\sigma^{2}>0$, and $L$ and $p$ are pre-specified. Note that for $t<L$, the maximum lag effect is of order $t-1$, and terms involving $X$ with non-positive subscript are not included in the model.

Model (1) is composed of two parts. First, for the structural component, the mean model is specified by lag coefficients $\boldsymbol{\beta}=(\beta_{0},...,\beta_{L})^{\prime}$ and control mean $\mu$. The immediate treatment effect is measured by $\beta_{0}$, and the carryover effect due to treatment on $l$ days ago is measured by $\beta_{l}$ for $l>0$. In the model, we assume the carryover effect beyond day $L$ is zero. As such, the total carryover treatment effect is captured as

Hence, the model naturally breaks down total treatment effect into $\beta_{0}$ and $\delta$.

Second, for the stochastic component, temporal dependency between errors is specified using an order-$p$ autoregressive error model with autoregression coefficient $\boldsymbol{\varphi}=(\varphi_{1},...,\varphi_{p})^{\prime}$. Let $B$ denote the backshift operator, that is, having $\Phi(B)=1-\varphi_{1}B-\varphi_{2}B^{2}-...-\varphi_{p}B^{p}$ so that the autoregression model for the error terms can be written as $\Phi(B)\epsilon_{t}=w_{t}$. It is often convenient to work with the transformed data $Y^{*}_{t}=\Phi(B)Y_{t}$ and $X_{t}^{*}=\Phi(B)X_{t}$ in the estimation steps. Thus, applying $\Phi(B)$ to both sides of model (1), we will rewrite the model

for $t=p+1,...,n$, where $\mu^{*}=\Phi(B)\mu$. To stack the data in vector form, we have

where $\boldsymbol{Y}^{*}=(Y^{*}_{p+1},...,Y^{*}_{n})^{\prime}$, $\boldsymbol{X}^{*}$ is a $(n-p)\times(L+1)$ matrix with $X^{*}_{k-l+p+1}$ being the $(k,l)$-th element of $\boldsymbol{X}^{*}$, $\boldsymbol{1}_{n-p}$ is a 1-vector of length $n-p$, and $\mathbf{I}_{n-p}$ is the identity matrix of dimension $n-p$. We denote $\tilde{\boldsymbol{\beta}}=(\mu,\boldsymbol{\beta}^{\prime})^{\prime}$ and $\tilde{\boldsymbol{X}^{*}}=(\Phi(B)\boldsymbol{1}_{n-p},\boldsymbol{X}^{*})$, so that $\tilde{\boldsymbol{X}^{*}}\tilde{\boldsymbol{\beta}}=\mu^{*}\boldsymbol{1}_{n-p}+\boldsymbol{X}^{*}\boldsymbol{\beta}$.

We consider normal prior distribution for $\tilde{\boldsymbol{\beta}}$, that is, having

where $\tilde{\mathbf{\Omega}}=diag(c_{0},\boldsymbol{\Omega})$ so that the prior variance of $\mu$ is $\sigma^{2}c_{0}^{-1}$ and the prior variance-covariance matrix of $\boldsymbol{\beta}$ is $\sigma^{2}\boldsymbol{\Omega}^{-1}$. We note that the prior variance depends on the variance $\sigma^{2}$ of the observations: such dependence renders a fused ridge penalized estimation procedure that is free of the variance parameters, resulting in computational stability; see Equation (7) below. We will postulate a non-informative prior on $\mu$ by setting $c_{0}$ to be a small number, and we will consider $\boldsymbol{\Omega}$ of the following form:

where the hyperparameters $\lambda_{l},\lambda_{l}^{*}>0$, for $l=0,\ldots,L,$ are constrained to increase over $l$. As a result of the monotonicity constraint, a lag coefficient $\beta_{l}$ at a greater lag $l$ is associated with a larger diagonal element in precision $\boldsymbol{\Omega}$, thus shrinking $\beta_{l}$ toward the prior mean (zero) to a greater extent. This effectively addresses collinearity of the lag coefficients without imposing strong parametric structure to $\boldsymbol{\beta}$. In addition, using the normal prior (5) with precision matrix (6), we can show that the maximum a posteriori probability estimate of $\boldsymbol{\tilde{\beta}}$ minimizes a fused ridge-type penalty:

where $\beta_{L+1}\triangleq 0$, thus giving insights on how the proposed prior constrains the lag coefficients: it regularizes not only the $\ell_{2}$-norm of the coefficients but also their successive differences, thereby enhancing local smoothness. The equivalence between the Bayesian inference and the fused ridge regularization (7) is proved in Appendix A.

There are many ways to specify $\lambda_{l}$ and $\lambda_{l}^{*}$ to meet the monotonicity constraints. In this article, we consider $\lambda_{l}=\exp\{\gamma_{1}(l+1)\}-1$ and $\lambda^{*}_{l}=\exp\{\gamma_{2}(l+1)\}-1$ for $\gamma_{1},\gamma_{2}>0$, so that $\gamma_{1}$ controls the rate at which the ridge penalty in (7) increases, and $\gamma_{2}$ controls the increasing rate of smoothness of the coefficient curve $\boldsymbol{\beta}$. Instead of treating these hyperparameters as fixed, we postulate a truncated standard exponential hyperprior on $(\gamma_{1},\gamma_{2})$, that is, having probability density function

where the support $S_{\boldsymbol{\gamma}}$ includes all pairs $(\gamma_{1},\gamma_{2})$ with which the precision matrix $\mathbf{\Omega}$ is positive definite. As such, the degree of ridge and smooth penalization can be determined according to the posterior distribution of the pair.

We put the Jeffreys prior for the error variance $\sigma^{2}$, that is, having density function

Note that any inverse-gamma prior for $\sigma^{2}$ would maintain conjugacy, and the Jeffreys prior can be regarded as an improper limit of inverse-gamma prior distribution.

For the autoregressive process, we consider a truncated normal prior for $\boldsymbol{\varphi}$ subject to the constraint that the error process is stationary. Specifically, we postulate

where ${S_{\boldsymbol{\varphi}}}(\boldsymbol{\varphi})$ denotes the support where all roots of the polynomial $\Phi(B)=1-\sum_{l=1}^{p}\varphi_{l}B^{l}$ are outside the unit circle. The process $\{\epsilon_{t}:t=1,2,\ldots\}$ is stationary when $\boldsymbol{\varphi}\in S_{\boldsymbol{\varphi}}(\boldsymbol{\varphi})$. [13] Note that the range of each $\varphi_{l}$ is $(-1,1)$; thus, a prior variance of 200 in (10) essentially amounts to a flat prior.

In this section, we evaluate the performance of the proposed BDLM-AR using simulation studies. At the end of each simulated trial, we fitted BDLM-AR with lag $L=7$ and AR(1), that is, having

where $\epsilon_{t}=\varphi\epsilon_{t-1}+w_{t}$ and $w_{t}\sim N(0,\sigma^{2})$. Posterior distributions were derived using the hybrid Metropolis Hastings/Gibbs algorithm described in the previous section with 50,000 iterations, a burn-in period of 25,000, and $a=0.2$ for sampling $\gamma$ in the MH step.

We compared BDLM-AR with some existing methods including the Bayesian distributed lag model (BDLagM),[10] Bayesian ridge DLM (BR-DLM) with a mean zero normal prior for $\tilde{\boldsymbol{\beta}}$, and a non-informative prior Bayesian DLM (NB-DLM) with a flat improper priors on each parameter in $\tilde{\boldsymbol{\beta}}$. These existing methods would use the same mean model (15) but assume independent errors without accounting for autocorrelation.

In addition, as a benchmark, we include the parametric Koyck’s DLM [7] which assumes the knowledge of the true autoregressive coefficients. For Bayesian models, we estimate the parameters using the posterior means and for Koyck model, we use the maximum likelihood estimates.

In each simulated N-of-1 trial, measurements were collected daily for 120 days, under one of two possible treatment sequences. In the first sequence, a participant would receive $x_{t}=1$ on the first 30 days and the last 30 days, and receive $x_{t}=0$ between days 31 and 90; that is,

In the second treatment sequence, a participant would switch treatments more frequently; specifically,

For each treatment sequence, the data were generated according to model (15) under five sets of lag coefficients (lag curves, LC):

1. LC1.

   Exponential decay curve: $\boldsymbol{\beta}=(5,2.5,1.25,0.625,0.3125,0,0,0)^{{}^{\prime}}$;

2. LC2.

   Exponential decay curve with oscillation: $\boldsymbol{\beta}=(5,2.5,-1.25,-0.625,0.3125,0,0,0)^{{}^{\prime}}$;

3. LC3.

   Slow absorption curve: $\boldsymbol{\beta}=(1.51,2.75,3.36,2.03,0.34,0,0,0)^{{}^{\prime}}$;

4. LC4.

   Slow absorption curve with oscillation: $\boldsymbol{\beta}=(1.51,2.75,-3.36,-2.03,0.34,0,0,0)^{{}^{\prime}}$;

5. LC5.

   No carryover effect: $\boldsymbol{\beta}=(10,0,0,0,0,0,0,0)^{{}^{\prime}}$.

The exponential decay curves (LC1 and LC2) specify coefficients that diminish in magnitude as lag lengthens. Specifically, the coefficients under LC1 decrease geometrically, which is aligned with the assumption of Koyck DLM. The slow absorption curves (LC3 and LC4) reflect scenarios where the carryover effect peaks at day 2 after treatment. LC5 is the null scenario where there is no carryover effect. The magnitudes of the coefficients were chosen in these scenarios so that $\|\boldsymbol{\beta}\|_{1}\approx 10$; in addition, the total carryover effects ($\delta$) are 4.69, 0.94, 8.48, -2.30 and 0 respectively for LC1-LC5. We consider $\sigma=10,20$ and $\varphi=0.5,0.2$ for the stochastic component in data generation. For each of these scenarios, the methods were evaluated using 100 simulated trials.

Figure 2 shows the bias and root mean squared error (RMSE) of the posterior means of individual lag coefficients. As expected, the biases of NB-DLM are relatively small; however, the method also has the largest RMSE uniformly because of the use of non-informative prior. The biases of the other methods are comparable, and are generally small compared to RMSE. While the $\ell_{2}$ penalty in BR-DLM on the lag coefficients reduces variability when compared to NB-DLM, the additional constraints on diminishing coefficients imposed by BDLagM and the proposed BDLM-AR further reduce RMSE for large lag $l$. Additionally, since the proposed BDLM-AR explicitly incorporates ridge-type regularization on the lag coefficients, it results in smaller RMSE for $\beta_{0}$ and the earlier lag coefficients (e.g. $\beta_{1}$). However, as a trade-off, the bias of BDLM-AR for early lag coefficients will be slightly inflated as compared to BDLagM, BR-DLM and NB-DLM, especially when true lag coefficient curve has frequent fluctuation. The benchmark method, Koyck DLM, performs best in the exponential decay case, where the true coefficient of autoregressive error ($\varphi$) is assumed to be known. The proposed BDLM-AR has very similar performance as Koyck DLM. Note that the coefficient of autoregressive error is estimated directly from the proposed BDLM-AR model, which is more practical in real application. In summary, the proposed BDLM-AR generally results in smallest RMSE for all lag coefficients.

To further examine the performance of each method in estimating the lag curve in aggregate, Figure 3 (top panel) plots the Euclidean distance between the vector of estimated lag coefficients and the vector of true lag coefficients. Under LC1 (exponential decay), the Koyck DLM has the best performance overall. This is not surprising because (i) the Koyck model mimics the coefficients under LC1 closely, and (ii) Koyck DLM assumes knowledge of the true autoregressive coefficients used in the simulation and hence it is not a method that can be implemented in practice. Thus, this comparison serves as a benchmark about the efficiency of the proposed BDLM-AR and other methods. The figure demonstrates that the proposed BDLM-AR produces smaller distance from the true lag curve than the other methods.

Table 1 gives the bias and RMSE in the estimation of the total effect ($\sum_{l=0}^{7}\beta_{l}$), the total carryover effect ($\delta=\sum_{l=i}^{7}\beta_{l}$) and the immediate effect ($\beta_{0}$) under different lag curves with $\sigma=10$ and $\varphi=0.5$ under treatment sequence $x_{t}^{(1)}$. Results for other values of $\varphi$, $\sigma$ and treatment sequence are similar and are provided in Figure A1 to A3 in the online Supporting Information. Overall, the proposed BDLM-AR yields consistently lower RMSE in estimating total effect, carryover effects, and immediate effects than other comparison methods, except for Koyck DLM. This is consistent with what we observe in Figures 2 and 3. We note that the advantages of BDLM-AR in terms of RMSE for the total carryover effect $(\delta)$ and immediate effect $(\beta_{0})$ are more pronounced than that for total effect ($\delta+\beta_{0}$). This is indeed the motivation that we set out to accomplish: to decompose the treatment effects and separate carryover effect from the immediate effect.

Figure 3 shows that the Euclidean distance between the vector of estimated lag coefficients and the truth under $x_{t}^{(2)}$ (bottom panel) is smaller than that under $x_{t}^{(1)}$ (top panel) suggesting frequently switching treatments will help improve the information content in N-of-1 trial data. This is in line with what we expect because collinearity of exposure lags will be lessened as treatments change frequently, while the total duration is held fixed. An implication in practice is that we should alternate intervention as frequent as it is feasible. The relative performance among methods is the same regardless of the treatment sequence, that is, the proposed BDLM-AR yields the shortest distance from the true lag coefficients $\boldsymbol{\beta}$.

In the previous subsections, BDLM-AR and other methods use a working mean model with $L=7$ and an AR(1) model for autoregressive errors. These working models correctly specify (or include the data generation model) in the previous simulation study. In this subsection, we investigate the robustness of BDLM-AR under model mis-specifications. Specifically, we will consider (1) the working mean model with $L=0,1,\ldots,7$; (2) the stochastic components that assume autoregressive error order of $p=0,1,7$. That is, we consider a total of 24 BDLM-AR models.

In data generation, we use LC1 as the true mean model, where $\beta_{l}>0$ for $l=0,1,2,3,4$, and we consider true scenarios for the errors:

1. E1.

   AR(1) with $\varphi=0.5$;

2. E2.

   Autoregressive model with $\varphi_{1}=0.5,\varphi_{2}=0,\varphi_{3}=0,\varphi_{4}=0.3,\varphi_{5}=0,\varphi_{6}=0.2$.

Note that, under the scenario LC1/E1, a working model with $L<4$ or $p=0$ under-specifies the true model. Likewise, under L1/E2, a working model with $L<4$ or $p=0,1$ under-specifies the true model.

Table 2 summarizes the RMSE of these 24 models under the two scenarios (LC1/E1 and LC/E2) with $\sigma=10$ under $x_{t}^{(1)}$. It can be seen that misspecified lag length has little influence on estimating total effect, total carryover effect and immediate effect, while under-specified error AR order will increase RMSE of parameters to a higher level than over-specified error AR order. Note that when choosing a small lag length value, we can hardly acquire estimation about the whole DL curve, as well as the information on the duration of carryover effect. Therefore, when the lag length is unknown, we suggest to fit data with a reasonable long lag length. For error autoregressive order, when the true orders are unknown, it is also suggested to fit a model with high autoregressive order.

Let $\pi(A|\cdot)$ denote conditional distribution of $A$ given all other variables in the model. The posterior distribution of $\tilde{\boldsymbol{\beta}}=(\mu,\boldsymbol{\beta}^{\prime})^{\prime}$ is

The MAP estimate of $\tilde{\boldsymbol{\beta}}$ is

Note that the model variance $\sigma^{2}$ is included in the prior distribution of $\tilde{\boldsymbol{\beta}}$ as a scaling parameter. In this way, the ridge and smoothness penalties are only controlled by $\lambda_{l}$ and $\lambda^{*}_{l}$, without depending on the model variance and prior covariance matrix of $\tilde{\boldsymbol{\beta}}$.