Comparison of Bayesian Networks, G-estimation and linear models to estimate causal treatment effects in aggregated N-of-1 trials with carry-over effects

Background The aggregation of a series of N-of-1 trials presents an innovative and efficient study design, as an alternative to traditional randomized clinical trials. Challenges for the statistical analysis arise when there is carry-over or complex dependencies of the treatment effect of interest. Methods In this study, we evaluate and compare methods for the analysis of aggregated N-of-1 trials in different scenarios with carry-over and complex dependencies of treatment effects on covariates. For this, we simulate data of a series of N-of-1 trials for Chronic Nonspecific Low Back Pain based on assumed causal relationships parameterized by directed acyclic graphs. In addition to existing statistical methods such as regression models, Bayesian Networks, and G-estimation, we introduce a carry-over adjusted parametric model (COAPM). Results The results show that all evaluated existing models have a good performance when there is no carry-over and no treatment dependence. When there is carry-over, COAPM yields unbiased and more efficient estimates while all other methods show some bias in the estimation. When there is known treatment dependence, all approaches that are capable to model it yield unbiased estimates. Finally, the efficiency of all methods decreases slightly when there are missing values, and the bias in the estimates can also increase. Conclusions This study presents a systematic evaluation of existing and novel approaches for the statistical analysis of a series of N-of-1 trials. We derive practical recommendations which methods may be best in which scenarios. Supplementary Information The online version contains supplementary material available at 10.1186/s12874-023-02012-5.

In our data simulation model, we used a similar model as described in Percha et al. [25] for the stochastic time series component. Let T denote the treatment variable, where T k (t) = 1 indicates that treatment k is given at time point t and T k (t) = 0 otherwise. We assume that the wash-in τ k and wash-out γ k of the treatment effect can be described by the following ordinary differential equation: This equation captures the decay of the wash-in and wash-out. Together with the linear effect of the treatment on the outcome, w T k ,O , and by setting Z t T k = E t k in equation (3) in the main manuscript, the direct carry-over adjusted treatment effect of treatment k at timepoint t on the outcome is w T k ,O · E t k . In the simulation, we incorporated exponential decay in the modeling of the full treatment effect at discrete time points t. For that, we used equation ES1 as the difference between two time points: As a result, we have a recursive formula for generating the exponential decay. This also shows that E t k can be understood as an exponential decay treatment indicator. Figure FS1 shows the simulated treatment effect with the exponential decay.
Further, a baseline drift B depending on time t is introduced through a Wiener process: [25] for more details). With that, we can simulate positive or negative time trends, when we set µ ≥ 0 or µ ≤ 0, respectively. Figures FS1 and FS2 show an illustration of the simulated data. In Figure FS1, the treatment effects w T k ,O E k (blue and orange) and the overall treatment effect ∑ k w T k ,O E k (green) are plotted. The patient has 4 treatment periods with a length of 28 days each and was exposed to T 1 , T 2 , T 1 , T 2 . Figure FS2 illustrates the different components of the outcome: the baseline drift (blue), the underlying state (orange) and the observed outcome values (green). In this example, the mean of the noise in the baseline drift µ b was set to 0 so the baseline does not change over time. If µ b were different from 0, then it would be expected that the baseline drift increases or decreases over time. For instance, if µ b were set to -1, it would be expected that the outcome on day i decreases compared to the previous day E(O i ) = E(O i−1 ) − 1.

Supplementary Text 2. DETAILS ON SIMULATION STUDY ON NONSPECIFIC LOW BACK PAIN
In the first step, we identified relevant variables as shown in Table TS1 and defined the variable measurements. The treatment was simulated as a binary variable identifying whether the patient was getting treatment 1 or 2, and had a certain effect size, wash-in effect and wash-out effect depending on the simulation scenario. All patients in the cohort were exposed to each treatment in each treatment block. We simulated a bigger effect for back stretching (treatment 2) in each data simulation by setting the effect from treatment 1 on the outcome to -2 and from treatment 2 on the outcome to -4. The outcome variable nonspecific low back pain was assessed on discrete pain levels on a scale from 0 to 15, where 15 indicates highest pain and 0 indicates no pain, drawn from a normal distribution through the Wiener process. To ensure that the simulated variables were in the defined range and discrete, they were truncated and rounded, if applicable.
As there are several scales to quantify the outcome [33], the simulated outcome can be interpreted as the score from the low back pain rating scale [33] divided by four. In all simulated data sets, the noise variable was generated with the same variance of 1. The mean µ b of the baseline drift was set to 0 and the standard deviation θ b to 0.2 to have an effect on the outcome over time.
We assumed that all patients were between 20 and 60 years old (drawn uniformly) and had a balanced male-female distribution.
As a proxy for the activity of a patients, Steps per Day (Activity) were simulated from a normal distribution N(µ = 8000, θ = 500) and the effects from Education, Demographics, Work and Chronic Diseases were added. The Quality of Sleep was measured on a scale of 0 (good quality) to 10 (bad quality), which could come from a questionnaire. The parameter was simulated with a truncated normal distribution N(µ = 5, θ = 3). Level of Stress could come as well from a questionnaire where the patient quantify the stress on a scale from 0 (no stress) to 10 (high stress). For the simulation, a Poisson distribution with λ = 2 was used. See Table TS1 for more details on the parameter values. Table TS1. Overview of the parameters used in the data simulation for the treatment, outcome and covariates. U(lower, upper) denotes a uniform distribution with a lower and an upper bound, λ is the decision parameter for the step function explained in the main manuscript and N(µ, σ) refers to a Normal distribution with corresponding values for the mean and the standard deviation. The simulation order of the variables in the data set matches the order of the variables in the table.

Category
Variable Name Simulation Parameter Baseline

Supplementary Text 3. SIMULATION OF MISSING VALUES
For simulating realistic data sets, missing values were introduced by randomly setting selected observations to None. The probability of setting a data point to missing was increased linearly over time so that the patient was more likely to have missing values as the trial progresses. With that, missing values are not completely at random, but as we have set the treatment order randomly we do not expect systematic bias on the treatment effect difference between the two treatments.
Furthermore, random blocks with missing values were added. These blocks had a predefined length and were inserted at a random time point with equal probability over time. In summary, for generating data with missing values, two parameters were specified in the simulation. The first parameter describes the dropout rate, which is the ratio of missing values to all data points. The second parameter defines the length of a block of missing values.

Supplementary Text 4. COAPM MODEL
In our data generation model, we used the following equation: For simplicity, we set L = [0] and substitute it in the Z t i=o : In the next step, we extract the treatment variables from the sum, so that j is a set of covariates C: In order to estimate w T 1 ,O and w T 2 ,O , we used the COAPM model defined as: Assuming that β 0 = E(U), we estimate w T k ,O withβ k in the COAPM model. Also, we estimate the treatment effect difference withβ 1 −β 2 .

Supplementary Text 5. FULL RESULTS OF THE SIMULATION STUDY
For each scenario, 100 samples with 5, 10, 25, 50 and 100 patients were drawn and analyzed separately. Estimates of the mean (µ) and standard error of the mean (σ) treatment effect difference were obtained using the different models.  Table TS4. Scenario 3: Activity interaction. Overview of the estimates (of the mean,μ, and standard error of the mean,σ) for the different sample sizes for the models: carry-over adjusted parametric model (COAPM), unadjusted Bayesian Network (BN), Bayesian Network with time adjustment (BNt), G-Estimation -independence (GEST-i), G-Estimation -AR1 (GEST-AR1), Linear Model (LM), and Sample Mean (SM).