Research in evidence-based medicine is increasingly moving towards informing individualised clinical care. This has led to renewed interest in N-of-1 clinical trials as they provide the strongest level of evidence for individual decisions (Guyatt et al., 2000). In N-of-1 trials, patients undergo a series of treatment periods called cycles where each patient receives each treatment (active or placebo) sequentially in a randomised order. This allows each treatment to be trailed on each patient, enabling patients to act as their own control. A major benefit of such a design is that individual and population treatment effects can be estimated through hierarchical modelling approaches as proposed by Zucker et al. (1997). However, one potential drawback of such a design is that it remains fixed throughout the entire trial. This is potentially limiting as data collected throughout the trial are not used to target informative treatment allocations. For this purpose, we propose a novel adaptive design algorithm for the selection of treatments that target the estimation of both population and individual treatment effects. Such an approach will therefore provide more informative and efficient N-of-1 trials, ensuring the right treatment is selected for the right patient in personalised care.
N-of-1 trials are typically used to test new and/or competing treatments for chronic diseases that are relatively stable over time. In typical N-of-1 clinical trials, treatments do not permanently change the disease or condition, and therefore, once off treatment, the patient will return to their underlying stable state (after a sufficiently long wash-out period). Such features may seem restrictive in practice, however, N-of-1 trials have been applied widely to inform clinical care for, for example, arthritis, asthma, insomnia, attention deficit hyperactivity disorder, hypertension, sleep disturbance, and fatigue from cancer (Pope et al., 2004; Coxeter et al., 2003; Nikles et al., 2006; Alemayehu et al., 2017; Nikles et al., 2019; Mitchell et al., 2015). Most notably, N-of-1 trials will typically recruit far fewer patients than randomised controlled trials making them suitable for smaller patient cohorts such as those with rare diseases (Stunnenberg et al., 2018; Bellgard et al., 2019). Further, given individual treatment effects can be estimated, N-of-1 trials are suitable when there is significant variability in response to treatment as seen, for example, in response to the treatment of chronic pain (Germini et al., 2017).
To extend the N-of-1 trial design, we consider a Bayesian adaptive design framework as sequential learning through time fits naturally within this framework, and is a framework where uncertainty about, for example, the preferred treatment for a given patient is handled most rigorously. A variety of different adaptive design algorithms have been proposed in the literature (see Ryan et al. (2016) for a recent review). Such algorithms often use Markov chain Monte Carlo (MCMC) samplers, importance sampling or sequential Monte Carlo (SMC) methods to update the prior information at each iteration of the adaptive design process
for a recent review). Such algorithms often use Markov chain Monte Carlo (MCMC) samplers, importance sampling or sequential Monte Carlo (SMC) methods to update the prior information at each iteration of the adaptive design process(Palmer & Müller, 1998; Weir et al., 2007; McGree et al., 2012, 2016). However, such approaches are not appropriate for designing N-of-1 trials as they are either too computationally expensive and/or do not allow both population and individual treatment effects to be estimated. As such, new statistical methods are needed that allow N-of-1 trials to be efficiently and adaptively designed.
For adaptively designing N-of-1 trials, we consider a Laplace approximation to both the log-likelihood and the posterior distribution of the parameters (Lin & Breslow, 1996; Skaug & Fournier, 2006). As will be shown, this provides an accurate approximation to the posterior distribution of the population parameters and individual random effects, and can be derived efficiently, without relying on Monte Carlo methods like MCMC which are typically computationally burdensome. The standard Laplace approximation to the posterior distribution has been considered previously (Lewi et al., 2009; Long et al., 2015; Senarathne et al., 2019) but has been limited to independent data settings, and, such as, is not appropriate for N-of-1 trials. For the selection of treatments, we employ the Kullback-Leibler divergence where treatments are randomised based on the probability of being the preferred treatment for a given individual. Of note, our Laplace approximation allows the MAB design to be employed efficiently through the availability of estimated individual random effects. Secondly, we compare our designs with those based on a randomised N-of-1 trial design where treatments are randomly selected with equal probability.
but has been limited to independent data settings, and, such as, is not appropriate for N-of-1 trials. For the selection of treatments, we employ the Kullback-Leibler divergence(Kullback & Leibler, 1951) between the prior and the posterior distribution of the population and individual random effect parameters. Such an approach targets the collection of data which is expected to provide the most information about the parameters. To evaluate this approach, two alternative designs are considered. The first is based on the multi-armed bandit (MAB) design (Thompson, 1933; Scott, 2010; Villar et al., 2015)
where treatments are randomised based on the probability of being the preferred treatment for a given individual. Of note, our Laplace approximation allows the MAB design to be employed efficiently through the availability of estimated individual random effects. Secondly, we compare our designs with those based on a randomised N-of-1 trial design where treatments are randomly selected with equal probability.
The outline of the paper is as follows. In Section 2, our motivating N-of-1 trial is introduced, along with the statistical model used to analyse the data from this trial. Our Bayesian adaptive design framework is then defined in Section 3. In Section 4, we show how to efficiently approximate the log-likelihood and the posterior distribution of the parameters. To illustrate our methodologies, Section 5 focuses on two examples of aggregated N-of-1 trials. The paper concludes with a discussion of key findings and suggestions for future research.
2 Motivating Example
The majority () of advanced cancer patients experience fatigue, with such fatigue being related to cancer treatment or disease itself (Vogelzang et al., 1997; Lawrence et al., 2004). Cancer-related fatigue (CRF) often persists after the end of treatment and can last for months, or even years (Bower et al., 2006). As reported in previous studies, CRF is more severe and persistent than normal fatigue caused by lack of sleep or overexertion, and has a negative impact on work, social relationships, and daily activities (Curt et al., 2000; Bower, 2014). Despite this, and the fact that such fatigue has significant detrimental effects on the quality of life of cancer patients, in most cases it is under-treated as most patients consider fatigue a symptom to be endured (Vogelzang et al., 1997).
Methylphenidate (MPH), a psychostimulant, is a commonly prescribed medication for the treatment of CRF. However, studies of this medication have yielded conflicting evidence about it’s effectiveness in treating CRF with little to no effect observed by Minton et al. (2011) while Kerr et al. (2012) observed significant improvement in advanced prostate cancer patients. Collectively, the results of such studies indicated that the effectiveness of MPH on CRF varies depending on the condition of the patient and the type of the cancer. Such variation motivated the consideration of an N-of-1 clinical trial of MPH by Mitchell et al. (2015). The main goals of the trial were to estimate: (1) The population treatment effect of MPH on CRF in patients with advanced cancer; (2) Individual treatment effects; and (3) How variable individual treatment effects are within a population of advanced cancer patients. In this paper, we consider this study as motivation for our methodological developments that will enable adaptive N-of-1 trials to be efficiently designed.
2.1 Data collection and selection
Mitchell et al. (2015) conducted a series of N-of-1 trials on 43 patients over three cycles. In each cycle, patients were assigned to both MPH (treatment) and placebo in a randomised order. To ensure patients (and clinicians) were blinded throughout the study, both the treatment and placebo were administered in capsules that were identical in appearance and taste. To measure CRF, the Functional Assessment of Chronic Illness Therapy-Fatigue (FACIT-F) subscale was used. This is a survey comprised of 13 questions where each response is measured on a five-point Likert scale. A total score (which is the primary outcome) is calculated as the sum across all responses, with higher scores indicating less fatigue.
As reported by Mitchell et al. (2015), there were 24 patients who completed the trial. Among these, only 22 patients completed all six treatment periods (i.e. yielded complete data from three cycles of two treatments), and these data were considered in this paper to estimate population and individual treatment effects, and the variability of treatment effects within and between patients.
2.2 Modelling aggregated N-of-1 trials
The analysis of aggregated N-of-1 clinical trial data can be undertaken within a mixed effects modelling framework. Such an approach not only allows population effects to be estimated (akin to that provided by randomised clinical trials) but also allows individual treatment effects to be estimated. Zucker et al. (2010) provide the foundations for this for N-of-1 clinical trials through a model specified as follows:
where denotes the observation of the patient at the cycle, are the population parameters, are the random effects for the patient, and and are the treatment allocation and the residual deviation for the observation of the patient at the cycle, respectively. Here, can take the value 1 or 0, depending on the treatment allocation (treatment=1, placebo=0). In this model, random effects and the residuals and , respectively.
After estimating the parameters in the above model, it can be used to assess population and individual treatment efficacy. Here, represents the population level average response when on placebo, and represents the population level difference between the average response on treatment compared to placebo. For estimating the individual treatment and placebo effects, the relevant individual random effect values are added to the population effects. As such, the average placebo and treatment effects for the patient can be expressed as and , respectively.
3 Adaptive design
Consider an adaptive N-of-1 clinical trial where the goal is to determine with respect to a placebo: (1) The effectiveness of a given treatment in the population; (2) The effectiveness of a given treatment for each individual; and (3) How variable the effectiveness of treatment is in the patient cohort. Without loss of generality, we assume that such an adaptive trial can be constructed by considering treatment selection for each cycle for each patient, iteratively. That is, initially treatments within the first cycle for the first patient will be determined. Then, once data have been collected from this cycle, treatments for the first cycle for the second patient will be determined, and so on. Of course, in practice, patients are recruited at different times, thus this ordering may not be exactly how a given trial would proceed. However, as will be seen, our proposed methodology does not depend on this adaptive structure, and is actually flexible enough to handle the large variety of design problems that may be encountered in real N-of-1 trials. If we suppose there are a total of treatments (including placebo), then the response represents the observation collected at the cycle for the patient given the design (treatment) for , and .
Then, the design can be expressed as follows:
where is the total number of patients in the study. Within an adaptive design framework, the Bayesian inference problem is to approximate the joint posterior distribution of the population parameters
is the total number of patients in the study. Within an adaptive design framework, the Bayesian inference problem is to approximate the joint posterior distribution of the population parametersand the random effects at each iteration of the adaptive design. That is, we are required to sample from or approximate the following:
where denotes the random effect variances, and denotes the conditional likelihood of observing data from design given the population parameters and the random effects , for , and . Further, is the distribution of the random effects conditional on the parameters , and is the prior distribution of the population parameters.
The Bayesian adaptive design problem can then be stated as selecting at each iteration based on a proposed utility function. Here, depending on the iteration, the design can be either , or . That is, in each iteration of this design algorithm, we select a treatment for either the same patient, the next patient or for the first patient in the study to start the next treatment cycle. For selecting which treatment to administer, a utility function is defined to reflect the aim of the study which we assume, based on the three goals stated at the start of Section 3, is parameter estimation. In general, we denote the utility function as , where is a supposed outcome obtained from running design . However, as , and
where the above expected utility is defined for continuous data as such data are observed in all examples considered in this paper. Extensions to other data types are straightforward.
When the utility function does not depend on the population parameters and the random effects , Equation (2) can be simplified to yield,
Thus, at each iteration of an adaptive design, one seeks to find , and this is termed the optimal design. Unfortunately, the above expression for the expected utility generally cannot be solved analytically, and thus needs to be approximated. The most common approach for this is Monte Carlo integration through the simulation of prior predictive data as follows:
The adaptive design process described above is outlined in Algorithm 1 where initially the prior information about the parameters is defined. Typically, in N-of-1 trials, this prior information will be uninformative or vague. Then, throughout the iterative process, the next optimal design point is found by maximising the expected utility (line 5), and the next data point is collected (line 6) based on this selected optimal design. The prior information about the population parameters and the random effects is then updated based on the information gained from the new data point (line 7). For the examples considered in this desktop study, data cannot actually be collected. In place of this, we assume data are generated from an underlying model with specified parameter values. For our motivating study, this underlying model is based on analysing data from Mitchell et al. (2015).
In considering the adaptive design process as outlined in Algorithm 1, there are two main challenges (at least computationally). The first is efficiently updating prior information as new data arrive (line 7). Employing an algorithm like MCMC would require re-running this at every iteration of the adaptive design. Given the high-dimensional nature of the posterior distribution (i.e. many random effects to estimate), this will quickly become computationally infeasible. Thus, alternative approaches are needed. The second difficulty is evaluating the expected utility function which requires sampling from or approximating a large number of posterior distributions, see Equation (4). Thus, efficient approaches are needed, and this motivates the development of the new methods proposed in this paper.
4 Efficient approximation to the joint posterior distribution
To efficiently approximate the posterior distribution of the parameters, a Laplace approximation is used. This approximation is formed by finding the posterior mode and evaluating the inverse of the negative Hessian matrix at this mode. These two terms then form the mean and variance-covariance matrix of a multivariate normal distribution, and it is this distribution which is used to approximate the posterior distribution. To describe how this Laplace approximation is formed in this paper, we first show how the posterior mode is located. For this purpose, we first define the likelihood for random effect models such as those defined in Equation (1). For such models, the likelihood function is defined by integrating out the random effects as follows:
where is the joint log-likelihood function for the population parameters and the random effects.
In general, there is no closed-form solution to the integral given in Equation (5), and thus evaluating the likelihood requires approximation. Following the work of Breslow & Clayton (1993), a Laplace approximation to the log-likelihood function can be formed as follows:
where and is the Hessian matrix evaluated at .
Based on the above approximation to the log-likelihood function, the posterior mode of the population parameters can be found as follows:
However, from the above formulation, it can be seen that is conditional on and is conditional on . Hence, to find the posterior mode for both and , an iterative process is used. That is, for an initial value of the population parameters, values of the random effects are estimated. These estimated random effects are then used to approximate the log-likelihood which is subsequently used to find the posterior mode for . This process then iterates until convergence, upon which and denote the mode of the posterior distribution. It is this mode that is taken as the mean of the multivariate normal approximation to the posterior distribution.
Once the posterior mode has been found, the Hessian matrix at this point can be evaluated and used to form the variance-covariance matrix of the multivariate normal distribution. For this, we note that the model specified in Equation (1) assumes the population parameters are independent of the random effects, hence this variance-covariance matrix will be block diagonal. The block corresponding to the random effects can be found via the Hessian matrix as denoted above (under Equation (6)), and the block corresponding to the population parameters can be found as follows:
The process for approximating the posterior distribution of the population parameters and the random effects is outlined in Algorithm 2. To initialise the algorithm, a value for the population parameters is randomly drawn from the prior distribution. Given these values , the mode of the random effects and the Hessian matrix at this mode are found (line 3). These two quantities are used to form a Laplace approximation to the log-likelihood function given in Equation (6). Next, this approximate log-likelihood function is used to approximate a density that is proportional to the log-posterior distribution of the population parameters. This density is then used to locate the posterior mode of the population parameters (line 4). This process iterates until convergence, where the algorithm outputs a (joint) multivariate normal approximation to the posterior distribution of the population and individual parameters (line 7).
5 Simulation studies
Here, two aggregated N-of-1 trials are considered to demonstrate the adaptive design approach proposed in Section 3 with the approximations given in Section 4. Since the main objective of these trials is to determine the preferred treatment at the population and individual patient level, the Kullback-Leibler divergence utility (Kullback & Leibler, 1951) was implemented for treatment selection. This utility function can be defined as follows:
where and denote the prior and the posterior distribution of the population parameters and the random effects, respectively.
When both prior and the posterior distribution follow a multivariate Normal distributions with mean vectors
follow a multivariate Normal distributions with mean vectorsand covariance matrices , respectively, then the KLD utility can be approximated as follows:
where is the dimension of the two multivariate Normal distributions.
Throughout the examples, one active treatment will be compared to placebo, and hence . Given this, treatment will be coded as an indicator variable with ‘1’ denoting active treatment and ‘0’ denoting placebo. For selecting the optimal design, as there are only two possibilities for design selection, the expected utility for each will be evaluated, and the choice which yields the largest value for the expected utility will be selected. For comparison, we benchmarked this optimal design approach against the MAB design and a randomised N-of-1 trial design.
A MAB design refers to a sequential experiment in which the goal is to determine the choice or ‘arm’ that yields the largest ‘reward’. In the context of a clinical trial, each arm represents a particular treatment (active or placebo), and reward refers to benefit from treatment. The adaptive design approach proposed in this paper provides a framework to select MAB designs based on individual patients, and this offers more flexibility than previous approaches which were based on population parameters only (Scott, 2010) . Here, benefit from treatment is based on the posterior predictive distribution for each patient. For example, if a higher response indicates a higher level of disease severity, then reward can be quantified via the following probability:
. Here, benefit from treatment is based on the posterior predictive distribution for each patient. For example, if a higher response indicates a higher level of disease severity, then reward can be quantified via the following probability:
where denotes the posterior mean for treatment . Within a Bayesian framework, the reward probability given in Equation (10) can be expressed as an expectation of an indicator function. Let if , and otherwise. Then,
This reward probability can be approximated by drawing a large number of samples from the joint posterior distribution of the population parameters and the random effects of the patient, and evaluating the following:
For the randomised N-of-1 design, the entire treatment sequence for each patient for all treatment cycles (Eg: ) can be obtained before starting the experiment. This treatment sequence is generated in such a way that the patient receives both treatment and placebo at each treatment cycle in a random order. As such, even though we use an adaptive design approach, randomised N-of-1 designs do not use the information collected throughout the trial to inform treatment selection. This is exactly how typical N-of-1 trials are run.
In each of the two examples considered in this paper, a simulation study was undertaken where experiments were sequentially simulated within our design framework. That is, in each cycle, the next data point for the patient was generated based on the selected optimal design and the assumed model (Equation (1)) with true parameter values (line 6 in Algorithm 1). Once the data have been generated, they will be used to update the prior information of the population and individual parameters. Designs for each patient within a given cycle will be found sequentially such that data points will be collected for a given patient, then data points will be collected for the next patient. Once all patients have completed the current cycle, a new cycle for the initial patient will begin. This will continue until all patients complete three cycles. For both examples, the prior distribution for the parameters was assumed to be vague, and follow Normal distributions as given in Table 1.
Lastly, as the results are subject to variability through the simulated data, all simulated studies were repeated 20 times to explore the range of outcomes that could be observed. After the designs were obtained and the corresponding data for each example was generated, the posterior distributions were re-evaluated using a standard MCMC approach. This re-evaluation step was undertaken to remove any potential bias from the Laplace approximation to the posterior distribution. We note, however, that the posterior distributions from MCMC and the Laplace approximation were very similar, see Figure S10 in the online supplementary material which shows a comparison of these posterior distributions for Example 1. All simulations were run using R 3.5.2, and R code to reproduce the results in this paper is available via the following GitHub repository, https://github.com/SenarathneSGJ/Adaptive_N-of-1_trials_design.
5.1 Example 1
In this example, we investigate the effect of an active treatment over placebo by running an adaptive N-of-1 trial with 20 patients. Here, the response variable is assumed to follow a Normal distribution with higher response values indicating a higher level of disease severity (e.g. pain as measured on a VAS scale). For each patient, six observations were collected over three treatment cycles, and the above defined three approaches for treatment selection were considered.
In this example, we investigate the effect of an active treatment over placebo by running an adaptive N-of-1 trial with 20 patients. Here, the response variable is assumed to follow a Normal distribution with higher response values indicating a higher level of disease severity (e.g. pain as measured on a VAS scale). For each patient, six observations were collected over three treatment cycles, and the above defined three approaches for treatment selection were considered.
Within this simulation study, four design scenarios were considered, each differing in terms of the parameter values assumed in the underlying generative model. Firstly, we considered a group of patients in which the between-subject variability was much smaller than the within-subject variability of the outcome, and where there was a significant difference between treatment and placebo at the population level. This scenario was considered as a baseline setting, where the remaining scenarios were defined by changing a parameter value in this baseline. Specifically, in Scenario 2, we increased the between-subject variability such that it equalled the within-subject variability of the outcome. Scenarios 3 and 4 were constructed by changing the population treatment effect such that there was a large and no difference (respectively) between treatment and placebo at the population level. These four scenarios are summarised in Table 2 by defining the parameter values used in each.
|Parameter||Scenario 1||Scenario 2||Scenario 3||Scenario 4|
For each design scenario, separate simulation studies were conducted. In each simulation study, a set of random effects were generated based on the assumed parameter values, and these were considered as the true effects for each patient in the study, and thus used to generate ‘real’ data throughout each N-of-1 trial. Furthermore, the same set of values were also used to calculate the true treatment effect for each patient, and hence to determine the best treatment (active or placebo) for each patient via the model described in Equation (1).
Results: Figures 1, 2 and 3 summarise the results from Scenario 1. Here, we first assessed the posterior precision of the population and random effect parameters when data were collected based on the optimal (KLD), MAB and randomised N-of-1 designs. For this comparison, we evaluated the log-determinant of the variance-covariance matrix of the joint posterior distribution after each cycle of the experiment. Figure 1 shows boxplots of the distribution of the log-determinant values of each intermediate posterior distribution (for each cycle) for all simulations. As can be seen, the posterior distributions obtained from optimal designs have lower log-determinant values compared to those obtained from MAB and randomised N-of-1 designs. This indicates that the posterior distributions obtained from the optimal design have higher precision compared to those obtained from either MAB or randomised N-of-1 design. As more data are collected, the MAB design also performs relatively well for estimation when compared to the randomised N-of-1 design.
Next, we assessed the designs in terms of identifying the best treatment allocation (active or placebo) for each patient in the study. For this purpose, we first calculated the true treatment effects based on the true parameter values as explained in Section 2.2. Then, based on the true treatment effects, the best treatment assignment for each patient was identified. Next, we evaluated the probability of identifying the best treatment for each patient based on each of the three designs. For this, we obtained a large number of samples from the joint posterior distribution of the population parameters and the random effects. Then, for each posterior sample, we calculated the individual mean response based on each treatment allocation. The required probability for each patient was then estimated via Equation (11) with . These probabilities were evaluated for all independent simulations, and the results were averaged to compare the performance of optimal, MAB and randomised N-of-1 designs.
Figure 2 shows the probability (with credible intervals) of identifying the best treatment allocation for each patient after each treatment cycle when the optimal, MAB and randomised N-of-1 designs were considered for Scenario 1. According to Figure 2, the optimal design has higher probability values (with less uncertainty) compared to the MAB and randomised N-of-1 designs. For patients who preferred placebo over the treatment (negative treatment difference), the median probability values were close to (not 1) for all designs in this scenario. The reason these probabilities are not closer to one is that, under the true model (and therefore what will most likely be inferred from the data), it is likely that for a given patient, the active treatment will be preferred. Hence, stronger evidence is needed to shift individual effects towards placebo when compared to patients who prefer the active treatment. This is particularly noticeable in this simulation study as the between subject variability is small compared to the within subject variability meaning patients who favour placebo are much less likely to occur than those who favour treatment. This can be seen by comparing these results to those from Scenario 2 where the between subject variability is large. As can be seen, probabilities for patients with similar treatment effect differences are generally closer to 1. Of note, this is not a feature of our approach to treatment selection, but rather a feature of all designs considered here including the randomised N-of-1 trial design.
Figure 3 compares the proportion of times the best treatment was administered for each patient in each cycle when the optimal, MAB and randomised N-of-1 designs were considered. Here, MAB design chose the best treatment for each patient a larger number of times than the other two design approaches. The optimal design also chose the best treatment for each patient a reasonable number of times except for the eleventh patient in the study. As the eleventh patient had the highest true treatment effect difference, the optimal design selected the placebo a large number of times. It is quite reasonable to observe such a difference as the optimal designs focus on learning about parameter values while MAB explores with a preference for the preferred treatment (based on currently available information). For this patient, it was more beneficial (in terms of learning about parameter values) to administer the placebo more often than not.
Similar to design Scenario 1, we compared the performance of the optimal, MAB and randomised N-of-1 designs under the remaining three design scenarios. For all of these three scenarios, the optimal designs were able to precisely estimate the joint posterior distribution of the population parameters and the random effects when compared to MAB and randomised N-of-1 designs (see Figures S1, S2 and S3 in the online supplementary material). As shown in Figures S4 and S6, when there was large uncertainty about the random effects (Scenario 2) or a small difference between the two treatments at the population level (Scenario 4), it was difficult to determine the best treatment assignment for patients where the difference between treatment and placebo was small. However, when there was a clear difference between the two treatments at the population level (Scenario 3), it was relatively straightforward to determine the best treatment allocation for each patient in the study (Figure S5). In all these scenarios, the optimal design performed better than MAB and randomised N-of-1 designs for determining the best treatment assignment for each patient. Figures S7, S8 and S9 in the online supplementary material compare the proportion of times the best treatment was given to each patient when treatments were assigned using the optimal, MAB and randomised N-of-1 designs under design Scenario 2, 3 and 4, respectively. As can be seen, for all scenarios, MAB design chose the best treatment for each patient more often than the other two design methods. However, this did not translate into providing more information about which treatment is better for each patient (as the optimal design performed best for this). Such results are expected based on how treatments are selected within each approach.
5.2 Motivating example
We now return to our motivating example of MPH trial introduced in Section 2. Based on the data collected from all 22 patients who completed all three cycles in Mitchell et al. (2015), population and individual parameter estimates were found by fitting the model described in Equation (1). A simulation study was then undertaken with treatment selection based on the optimal, MAB and randomised N-of-1 designs. This simulation study was run exactly as outlined for Example 1, including the prior distributions used (see Table 1 ). Here, the response variable was assumed to follow a log-normal distribution with higher values indicating less fatigue (higher recovery).
). Here, the response variable was assumed to follow a log-normal distribution with higher values indicating less fatigue (higher recovery).
Results: For the purpose of comparing the parameter estimation results, log-determinant values of the posterior variance-covariance matrix of the population parameters and the random effects were evaluated and plotted. As shown in Figure 4, the joint posterior distribution based on the optimal designs have smaller log-determinant values compared to those based on MAB and randomised N-of-1 designs. That is, the optimal design was relatively efficient for parameter estimation. Again, as more design points were collected, MAB design performed relatively well for estimation when compared to the randomised N-of-1 design.
In Figure 5, we compare the probability (with credible intervals) of identifying the best treatment assignment for each patient after each treatment cycle when the optimal, MAB and randomised N-of-1 designs were considered for data collection. Similar to the first example, these probabilities were approximated using Equation (11), but with a different indicator function. Here, the indicator function equals to 1 if the mean response value of the best treatment assignment for a given patient is higher than that of the remaining treatment assignment, and 0 otherwise. As can be seen, the optimal design performed relatively well for this experimental goal when compared to the other two design methods.
The proportion of times the best treatment was selected for each patient in each treatment cycle is shown in Figure 6. Again, MAB design selected the best treatment for each patient a larger number of times when compared to the other two design methods.
In this work, we have developed a Bayesian adaptive design approach to find optimal treatment allocations for N-of-1 trials. As was seen, the designs derived from our approach can be used to determine the best treatment assignment for each patient with a fewer number of treatment cycles and/or provide more certainty about this after three cycles (typical duration of an N-of-1 trial). The empirical evidence presented in this paper demonstrates that the proposed Bayesian adaptive framework can be used to efficiently estimate both population parameters and the random effects in realistically sized N-of-1 trials, and benefits of this approach over alternatives were demonstrated. As such, we propose this method could be adopted in future N-of-1 trials to determine appropriateness in real-world settings.
In the first example, four design scenarios were considered to investigate the performance of our adaptive design approach under different parameter settings. It was found that the optimal designs were preferred for estimating the population parameters and the random effects when compared to the MAB and randomised N-of-1 designs. Furthermore, in using the optimal designs, we were able to determine the best treatment assignment for each patient in a fewer number of treatment cycles. When there was considerable variability in the random effects and a small difference between the individual treatment effects, it was difficult to determine the best treatment assignment with our optimal design approach but this was also observed for the other two approaches considered in this paper. In all four scenarios, MAB designs chose the best treatment for each patient a larger number of times than the optimal and randomised N-of-1 designs but we note that this did not translate into more certainty about which treatment was best for each patient.
When we considered the motivating example for this research, benefits were seen in adopting our optimal design approach compared to the MAB and randomised N-of-1 designs in terms of estimating model parameters. Of note, this again translated into more certainty about the preferred treatment for each patient during and at the end of the study. Thus, it seems as though more information from the N-of-1 study of MPH could have been obtained if Bayesian adaptive design methods were implemented. Accordingly, we hope to explore the use of our methods in real N-of-1 trials into the future.
The two Laplace approximations proposed in this paper to form an approximation to the posterior distribution of the parameter is different to what has previously been proposed in the literature by, for example, Overstall et al. (2018). In such work, authors have proposed a single Laplace approximation formed by considering the conditional likelihood i.e. not integrating out the random effects. Given this, it is of interest to compare the two approaches to determine which appears to yield a better approximation to the posterior distribution. To investigate this, a separate simulation study was undertaken where posterior distributions obtained from both approaches were compared to that obtained from MCMC. For this, the data generating model defined in Example 1 with five patients was used to simulate 50 data sets, each based on a typical N-of-1 design with 3 cycles. For each simulated data set, MCMC and the two Laplace-based approaches were used to form an approximation to the posterior distribution of the parameters. The posterior mean and variance for all parameters were recorded for each simulation, and the distribution of these is summarised in Table 3. As can be seen, the posterior distributions obtained from MCMC and the Laplace approach proposed in this paper have similar mean values, while the approximation from Overstall et al. (2018) are noticeably different, particularly for the random effect parameters. Both Laplace approximations appear to underestimate the variance of the parameters but this is much more apparent when the approach from Overstall et al. (2018) is used. An example of the posterior distributions obtained under these three approaches is shown in Figure 7. The discussed advantages of our Laplace methods are highlighted in this plot.
Future development of our adaptive design approach could include extensions to other types of trials. Of note, cross-over, single case experimental designs and step-wedge designs can be viewed as special cases of the N-of-1 trial design. Thus, our approach could potentially be adopted within such settings. Further, it would be interesting to explore our methodologies for designing platform trials (Berry et al., 2015). Such trials generally consider more treatments when compared to N-of-1s, and the availability of different treatments can varying depending on the patient. Our approach to targeting information at the population and individual patient level could prove useful in, for example, quickly discounting ineffective treatments. We plan to explore this in future endeavours.
When designing N-of-1 trials, typically there is some information available about treatment effects from previously collected data. Such prior information can be used to determine an appropriate model for which to undertaken adaptive design. However, if not sufficient, then uncertainty at the model level should be incorporated into the selection of treatments. Of note, the identified Bayesian adaptive design approach can be extended to incorporate model uncertainty by forming a non-trivial prior on the model space. Further, in cases where it is desirable to (additionally) learn which model is most appropriate for the data, dual-purpose utility functions could be considered (Borth, 1975; McGree, 2017). In cases where different data types and multiple primary and secondary outcomes are considered, extensions to the work of Senarathne et al. (2019) may prove fruitful in forming a more accurate approximation to the posterior distribution. These are also areas of research which we hope to explore in the future.
SGJS was supported by QUTPRA scholarship from the Queensland University of Technology. Computational resources and services used in this work were provided by the HPC and Research Support Group, Queensland University of Technology, Brisbane, Australia.
- Alemayehu et al. (2017) Alemayehu, C., Mitchell, G., Aseffa, A., Clavarino, A., McGree, J. & Nikles, J. (2017), ‘A series of N-of-1 trials to assess the therapeutic interchangeability of two enalapril formulations in the treatment of hypertension in addis ababa, ethiopia: study protocol for a randomized controlled trial’, Trials 18(1). 470.
- Bellgard et al. (2019) Bellgard, M., Snelling, T. & McGree, J. (2019), ‘Rd-rap: beyond rare disease patient registries, devising a comprehensive data and analytic framework’, Orphanet Journal of Rare Disease 14(1). 176.
- Berry et al. (2015) Berry, S. M., Connor, J. T. & Lewis, R. J. (2015), ‘The platform trial: an efficient strategy for evaluating multiple treatments’, JAMA 313(16), 1619–1620.
- Borth (1975) Borth, D. M. (1975), ‘A total entropy criterion for the dual problem of model discrimination and parameter estimation’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 37, 77–87.
- Bower (2014) Bower, J. E. (2014), ‘Cancer-related fatigue: mechanisms, risk factors, and treatments’, Nature reviews. Clinical Oncology 11(10), 597–609.
- Bower et al. (2006) Bower, J. E., Ganz, P. A., Desmond, K. A., Bernaards, C., Rowland, J. H., Meyerowitz, B. E. & Belin, T. R. (2006), ‘Fatigue in long-term breast carcinoma survivors’, Cancer 106(4), 751–758.
Breslow & Clayton (1993)
Breslow, N. E. & Clayton, D. G. (1993), ‘Approximate inference in generalized linear mixed models’,Journal of the American Statistical Association 88(421), 9–25.
- Coxeter et al. (2003) Coxeter, P., Schluter, P., Eastwood, H., Nikles, C. & Glasziou, P. (2003), ‘Valerian does not appear to reduce symptoms for patients with chronic insomnia in general practice using a series of randomised N-of-1 trials’, Complementary Therapies in Medicine 11(4), 215 – 222.
- Curt et al. (2000) Curt, G. A., Breitbart, W., Cella, D., Groopman, J. E., Horning, S. J., Itri, L. M., Johnson, D. H., Miaskowski, C., Scherr, S. L., Portenoy, R. K. et al. (2000), ‘Impact of cancer-related fatigue on the lives of patients: new findings from the fatigue coalition’, The Oncologist 5(5), 353–360.
- Germini et al. (2017) Germini, F., Coerezza, A., Andreinetti, L., Nobili, A., Rossi, P. D., Mari, D., Guyatt, G. & Marcucci, M. (2017), ‘N-of-1 randomized trials of ultra-micronized palmitoylethanolamide in older patients with chronic pain’, Drugs & Aging 34(12), 941–952.
- Guyatt et al. (2000) Guyatt, G. H., Haynes, R. B., Jaeschke, R. Z., Cook, D. J., Green, L., Naylor, C. D., Wilson, M. C., Richardson, W. S. & for the Evidence-Based Medicine Working Group (2000), ‘Users’ Guides to the Medical Literature:XXV. Evidence-Based Medicine: Principles for Applying the Users’ Guides to Patient Care’, JAMA 284(10), 1290–1296.
- Kerr et al. (2012) Kerr, C. W., Drake, J., Milch, R. A., Brazeau, D. A., Skretny, J. A., Brazeau, G. A. & Donnelly, J. P. (2012), ‘Effects of methylphenidate on fatigue and depression: A randomized, double-blind, placebo-controlled trial’, Journal of Pain and Symptom Management 43(1), 68 – 77.
- Kullback & Leibler (1951) Kullback, S. & Leibler, R. A. (1951), ‘On information and sufficiency’, The Annals of Mathematical Statistics 22(1), 79–86.
- Lawrence et al. (2004) Lawrence, D. P., Kupelnick, B., Miller, K., Devine, D. & Lau, J. (2004), ‘Evidence report on the occurrence, assessment, and treatment of fatigue in cancer patients’, JNCI Monographs 2004(32), 40–50.
- Lewi et al. (2009) Lewi, J., Butera, R. & Paninski, L. (2009), ‘Sequential optimal design of neurophysiology experiments’, Neural Computation 21(3), 619–687.
- Lin & Breslow (1996) Lin, X. & Breslow, N. E. (1996), ‘Bias correction in generalized linear mixed models with multiple components of dispersion’, Journal of the American Statistical Association 91(435), 1007–1016.
- Long et al. (2015) Long, Q., Scavino, M., Tempone, R. & Wang, S. (2015), ‘A Laplace method for under-determined Bayesian optimal experimental designs’, Computer Methods in Applied Mechanics and Engineering 285, 849–876.
- McGree (2017) McGree, J. M. (2017), ‘Developments of the total entropy utility function for the dual purpose of model discrimination and parameter estimation in Bayesian design’, Computational Statistics and Data Analysis 113, 207–225.
- McGree et al. (2012) McGree, J. M., Drovandi, C. C., Thompson, M., Eccleston, J., Duffull, S., Mengersen, K., Pettitt, A. N. & Goggin, T. (2012), ‘Adaptive Bayesian compound designs for dose finding studies’, Journal of Statistical Planning and Inference 142(6), 1480–1492.
- McGree et al. (2016) McGree, J. M., Drovandi, C. C., White, G. & Pettitt, A. N. (2016), ‘A pseudo-marginal sequential Monte Carlo algorithm for random effects models in Bayesian sequential design’, Statistics and Computing 26(5), 1121–1136.
- Minton et al. (2011) Minton, O., Richardson, A., Sharpe, M., Hotopf, M. & Stone, P. C. (2011), ‘Psychostimulants for the management of cancer-related fatigue: a systematic review and meta-analysis’, Journal of Pain and Symptom Management 41(4), 761 – 767.
- Mitchell et al. (2015) Mitchell, G. K., Hardy, J. R., Nikles, C. J., Carmont, S. S., Senior, H. E., Schluter, P. J., Good, P. & Currow, D. C. (2015), ‘The effect of methylphenidate on fatigue in advanced cancer: an aggregated N-of-1 trial’, Journal of Pain and Symptom Management 50(3), 289 – 296.
- Nikles et al. (2006) Nikles, C. J., Mitchell, G. K., Del Mar, C. B., Clavarino, A. & McNairn, N. (2006), ‘An N-of-1 trial service in clinical practice: Testing the effectiveness of stimulants for attention-deficit/hyperactivity disorder’, Pediatrics 117(6), 2040–2046.
- Nikles et al. (2019) Nikles, J., O’Sullivan, J., Mitchell, G., Smith, S., McGree, J., Senior, H., Dissanyaka, N. & Ritchie, A. (2019), ‘Protocol: Using N-of-1 tests to identify responders to melatonin for sleep disturbance in parkinson’s disease’, Contemporary Clinical Trials Communications 15, 100397.
Overstall et al. (2018)
Overstall, A. M., McGree, J. M. & Drovandi, C. C. (2018), ‘An approach for finding fully Bayesian optimal designs using normal-based approximations to loss functions’,Statistics and Computing 28(2), 343–358.
- Palmer & Müller (1998) Palmer, J. L. & Müller, P. (1998), ‘Bayesian optimal design in population models for haematologic data’, Statistics in Medicine 17(14), 1613–1622.
- Pope et al. (2004) Pope, J. E., Prashker, M. & Anderson, J. (2004), ‘The efficacy and cost effectiveness of N-of-1 studies with diclofenac compared to standard treatment with nonsteroidal antiinflammatory drugs in osteoarthritis.’, The Journal of Rheumatology 31(1), 140–149.
- Ryan et al. (2016) Ryan, E. G., Drovandi, C. C., McGree, J. M. & Pettitt, A. N. (2016), ‘A review of modern computational algorithms for bayesian optimal design’, International Statistical Review 84(1), 128–154.
- Scott (2010) Scott, S. L. (2010), ‘A modern Bayesian look at the multi-armed bandit’, Applied Stochastic Models in Business Industry 26(6), 639–658.
- Senarathne et al. (2019) Senarathne, S. G. J., Drovandi, C. C. & McGree, J. M. (2019), ‘Bayesian sequential design for Copula models’, TEST .
- Skaug & Fournier (2006) Skaug, H. J. & Fournier, D. A. (2006), ‘Automatic approximation of the marginal likelihood in non-gaussian hierarchical models’, Computational Statistics & Data Analysis 51(2), 699 – 709.
- Stunnenberg et al. (2018) Stunnenberg, B. C., Raaphorst, J., Groenewoud, H. M., Statland, J. M., Griggs, R. C., Woertman, W., Stegeman, D. F., Timmermans, J., Trivedi, J., Matthews, E., Saris, C. G. J., Schouwenberg, B. J., Drost, G., van Engelen, B. G. M. & van der Wilt, G. J. (2018), ‘Effect of Mexiletine on Muscle Stiffness in Patients With Nondystrophic Myotonia Evaluated Using Aggregated N-of-1 Trials’, JAMA 320(22), 2344–2353.
- Thompson (1933) Thompson, W. R. (1933), ‘On the likelihood that one unknown probability exceeds another in view of the evidence of two samples’, Biometrika 25(3/4), 285–294.
- Villar et al. (2015) Villar, S. S., Bowden, J. & Wason, J. (2015), ‘Multi-armed bandit models for the optimal design of clinical trials: Benefits and challenges’, Statistical science : a review journal of the Institute of Mathematical Statistics 30(2), 199–215.
- Vogelzang et al. (1997) Vogelzang, N., Breitbart, W., Cella, D., Curt, G., Groopman, J., Horning, S., Itri, L., Johnson, D., Saherr, S. & Portenoy, R. (1997), ‘Patient, caregiver, and oncologist perceptions of cancer-related fatigue: results of a tripart assessment survey’, Seminars in Hematology 34(3), 4–12.
- Weir et al. (2007) Weir, C. J., Spiegelhalter, D. J. & Grieve, A. P. (2007), ‘Flexible design and efficient implementation of adaptive dose-finding studies’, Journal of Biopharmaceutical Statistics 17(6), 1033–1050.
- Zucker et al. (2010) Zucker, D. R., Ruthazer, R. & Schmid, C. H. (2010), ‘Individual (N-of-1) trials can be combined to give population comparative treatment effect estimates: methodologic considerations’, Journal of Clinical Epidemiology 63(12), 1312–1323.
- Zucker et al. (1997) Zucker, D. R., Schmid, C. H., McIntosh, M. W., D’Agostino, R. B., Selker, H. P. & Lau, J. (1997), ‘Combining single patient (N-of-1) trials to estimate population treatment effects and to evaluate individual patient responses to treatment’, Journal of Clinical Epidemiology 50(4), 401–410.
|23.98 (23.00, 25.03)||1.16 (0.64, 1.92)||23.97 (22.99, 25.03)||2.25 (1.04, 4.47)||23.98 (22.94, 25.13)||0.88 (0.37, 1.47)|
|-1.94 (-3.4, -0.29)||1.93 (1.14, 3.35)||-1.93 (-3.42, -0.29)||3.44 (1.81, 6.64)||-1.94 (-3.34, -0.32)||1.45 (0.73, 2.69)|
|1.04 (0.75, 1.25)||0.02 (0.02, 0.03)||1.09 (0.80, 1.30)||0.02 (0.02, 0.03)||1.07 (0.71, 1.40)||0.02 (0.02, 0.02)|
|0.42 (-0.24, 0.98)||0.32 (0.17, 0.49)||0.46 (-0.22, 1.16)||0.57 (0.29, 0.74)||-2.59 (-4.16, 0.82)||1.86 (0.17, 2.67)|
|0.58 (0.14, 1.18)||0.41 (0.17, 0.54)||0.56 (0.10, 1.34)||0.67 (0.33, 0.81)||-2.31 (-2.99, 1.05)||2.35 (0.14, 2.75)|
|b01||1.17 (-0.03, 3.17)||1.01 (0.46, 1.67)||1.15 (-0.01, 3.11)||2.84 (1.04, 5.63)||0.77 (0.00, 3.2)||0.59 (0.00, 2.01)|
|b02||0.04 (-1.12, 1.19)||1.00 (0.46, 1.66)||0.05 (-1.08, 1.18)||2.65 (0.98, 5.18)||-0.07 (-1.63, 0.81)||0.52 (0.00, 1.80)|
|b03||0.20 (-1.15, 1.18)||1.00 (0.46, 1.66)||0.20 (-1.21, 1.16)||2.62 (0.99, 5.24)||0.10 (-0.59, 1.20)||0.50 (0.00, 1.86)|
|b04||-1.65 (-3.65, -0.21)||1.01 (0.47, 1.70)||-1.62 (-3.67, -0.24)||3.06 (0.98, 5.64)||-0.91 (-3.39, 0.00)||0.64 (0.00, 2.23)|
|b05||0.25 (-0.72, 1.34)||1.00 (0.47, 1.66)||0.25 (-0.67, 1.27)||2.61 (1.00, 5.10)||0.12 (-0.19, 0.88)||0.50 (0.00, 1.81)|
|b11||1.19 (-0.51, 3.96)||1.73 (0.96, 3.03)||1.18 (-0.59, 3.89)||4.38 (1.93, 8.71)||0.53 (0.00, 3.96)||0.44 (0.00, 3.64)|
|b12||-0.82 (-3.21, 0.57)||1.73 (0.96, 2.93)||-0.82 (-3.20, 0.61)||4.18 (1.90, 8.34)||-0.33 (-3.01, 0.00)||0.40 (0.00, 3.00)|
|b13||0.04 (-1.78, 1.57)||1.74 (0.95, 3.10)||0.04 (-1.81, 1.58)||4.02 (1.84, 7.69)||0.01 (-0.49, 0.38)||0.42 (0.00, 3.09)|
|b14||-0.81 (-3.56, 0.73)||1.72 (0.98, 3.03)||-0.82 (-3.42, 0.82)||4.39 (1.87, 9.35)||-0.34 (-3.20, 0.00)||0.43 (0.00, 3.13)|
|b15||0.39 (-1.19, 2.07)||1.75 (0.98, 3.06)||0.38 (-1.27, 1.90)||4.00 (1.85, 8.00)||0.18 (0.00, 1.28)||0.40 (0.00, 3.01)|
Supplementary Material for “Bayesian Adaptive N-of-1 Trials for Estimating Population and Individual Treatment Effects”