Estimating mono- and bi-phasic regression parameters using a mixture piecewise linear Bayesian hierarchical model

Rui Zhao; Paul Catalano; Victor G. DeGruttola; Franziska Michor

doi:10.1371/journal.pone.0180756

Abstract

The dynamics of tumor burden, secreted proteins or other biomarkers over time, is often used to evaluate the effectiveness of therapy and to predict outcomes for patients. Many methods have been proposed to investigate longitudinal trends to better characterize patients and to understand disease progression. However, most approaches assume a homogeneous patient population and a uniform response trajectory over time and across patients. Here, we present a mixture piecewise linear Bayesian hierarchical model, which takes into account both population heterogeneity and nonlinear relationships between biomarkers and time. Simulation results show that our method was able to classify subjects according to their patterns of treatment response with greater than 80% accuracy in the three scenarios tested. We then applied our model to a large randomized controlled phase III clinical trial of multiple myeloma patients. Analysis results suggest that the longitudinal tumor burden trajectories in multiple myeloma patients are heterogeneous and nonlinear, even among patients assigned to the same treatment cohort. In addition, between cohorts, there are distinct differences in terms of the regression parameters and the distributions among categories in the mixture. Those results imply that longitudinal data from clinical trials may harbor unobserved subgroups and nonlinear relationships; accounting for both may be important for analyzing longitudinal data.

Citation: Zhao R, Catalano P, DeGruttola VG, Michor F (2017) Estimating mono- and bi-phasic regression parameters using a mixture piecewise linear Bayesian hierarchical model. PLoS ONE 12(7): e0180756. https://doi.org/10.1371/journal.pone.0180756

Editor: Zhongxue Chen, Indiana University Bloomington, UNITED STATES

Received: February 27, 2016; Accepted: June 21, 2017; Published: July 19, 2017

Copyright: © 2017 Zhao et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: Simulation data can be re-generated using procedures described in the paper; Data from VISTA study may be requested via Takeda pharmaceutical. The clinical trial data used in this paper were published in another publication, http://clincancerres.aacrjournals.org/content/22/16/4206.long.

Funding: This work was supported by the Dana-Farber Cancer Institute Physical Sciences Oncology Center (U54CA143798) and NIH grant R37 51164. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Mixed-effects models are particularly useful in medical research because of their ability to handle imbalances in the number of observations across patients and to identify between-subject and within-subject sources of variability [1–3]. Progress to extend mixed-effects models to include heterogeneity in data has been made by incorporating a finite mixture into the model [4, 5]. This extension is particularly relevant to clinical research, since clinical data often contain unobserved categorical variables corresponding to, for example, “responders” or “non-responders” to a given treatment. Ignoring such mixtures may result in biases in estimates. Xu and Hedeker investigated this idea and found there is ample evidence of non-homogeneous responses in two large psychiatric clinical trials [6]. Ketchum et al. further extended the mixed-effects mixture models to allow for differences in the variance-covariance matrices [7]. These improvements enable the random-effects models to better characterize heterogeneity in data. A book chapter by Verbeke and Molenberghs provides an excellent summary of heterogeneous mixed models [8].

In addition to population heterogeneity, changes in functional relationships between response variables and explanatory variables, particularly with time, are ubiquitous in longitudinal studies: HIV-1 viral load [9, 10], hepatitis B/C viral load [11, 12] and BCR-ABL expression levels in chronic myeloid leukemia [13, 14]. In those examples, biomarkers exhibit nonlinear changes over time, many of which are bi-phasic in nature—that is, patients biomarkers have two distinct patterns over time rather than one uniform trajectory. An example of bi-phasic decline patterns is that in some chronic myeloid leukemia patients, the initial decline of BCR-ABL expression levels is much faster than later declines, whereas in other patients the decline is uniform over time [13, 14]. These observations are contrary to the assumptions of many models that parameters are invariant over time. One method for accounting for changes in the longitudinal relationships over time is provided by nonlinear mixed-effects models [15, 16]. Morrel et al. applied a piecewise nonlinear mixed-effects model to a prostate cancer data set and found that patients with local lesions and metastatic lesions have similar initial prostate-specific antigen (PSA) trajectories. However, they found that the rates of PSA increase in a later phase were larger in patients with metastatic lesions than patients with local lesions. Naumova et al. used a piecewise mixed-effect model to analyze a prospective study on the development of obesity in female adolescents [17]. Cudeck and Klebe, and Harring et al. applied similar ideas to psychology-related data sets [18, 19]. Those examples demonstrate the flexibility of nonlinear mixed-effects models in investigating changing functional relationships over time.

Both heterogeneity in patient populations and changes in the longitudinal relationship have been addressed separately in several publications [6, 7, 15–19]; however, only a few publications have tackled both problems simultaneously. Pauler and Laird introduced a general framework for finite mixtures of nonlinear hierarchical models; they applied their methods to investigate non-compliance in a HIV clinical trial [20]. In their application, the mixture consists of a constant mean model for the compliant patients and a piecewise linear model for the non-compliant patients. Recently, Lu and Huang extended the general framework proposed by Pauler and Laird [20] to incorporate skewness in the distributions of individual regression parameters, relaxing the normal assumption [21]; they applied their methods to analyze a HIV viral load data set [22]. Their underlying nonlinear mixed-effects model was formulated based on the model structure of an ordinary differential equation (ODE) model describing viral loads over time [23]. One caveat associated with those approaches is that their models require extensive prior biological knowledge in order to specify the nonlinear models before analyzing the data. Misspecification of the model may have detrimental effects on parameter estimation and patient classification. Particularly, if the differences between categories in the mixture are not well separated, specifying the model becomes an even more challenging problem.

To address this issue, we developed a piecewise linear random-effects mixture model that does not require any prior knowledge on the model structure to account for both heterogeneity and change in the longitudinal relationship over time. The only assumptions of this model are that the underlying data may contain a mixture of mono- and bi-phasic observations, and that the bi-phasic observations are piecewise linear; no further constraints on the intercepts and slopes are necessary. The primary purpose of this model is to detect unobserved subgroups in a patient population and nonlinear longitudinal relationships over time. This method is particularly useful for current clinical trials, which often include diagnostic and prognostic hypotheses. With periodical follow-up visits gathering biomarker data, parameters associated with heterogeneous changes in biomarker trends can be detected with this method. Given the usually limited number of follow-up measurements in clinical trials, the current implementation of this model focuses primarily on mono- vs. bi-phasic changes; however, our model can easily be generalized to include multi-phasic changes and multi-category mixtures. In addition, because of the piecewise linear nature of the proposed model, other clinically relevant covariates can also easily be included.

Materials and methods

We designed a mixture piecewise linear Bayesian hierarchical model to estimate regression parameters and to determine the posterior distributions of these parameters, while accounting for both population heterogeneity and changes in longitudinal relationships over time. We consider situations in which the patient population consists two latent subpopulations: mono-phasic and bi-phasic patients. The individual-level trajectories for mono-phasic and bi-phasic patients are shown in Eq (1). Throughout the text, we use subscripts S and B to denote mono- (i.e. single-) and bi-phasic patients, respectively. For the i^th patient with a total of M_i observations, the dependent variable y_ij, corresponding to the quantitative measure of disease burden, which may either follow a mono-phasic or a bi-phasic regression line, depending on the latent indicator variable η_i: (1)

Here, ε_ij denotes the independent error term, which follows a normal distribution centered at 0 with variance σ²; η_i denotes the phasic indicator for patient i, with 0 and 1 denoting mono- and bi-phasic patterns, respectively; k_i, a latent variable, denotes the number of observations belonging to the first phase for patient i, if the response of patient i is bi-phasic. For the individual regression parameters, we assume hierarchical normal distributions for s_i and b_i: (2) (3)

We first consider the artificial case in which we do not know if a patient follows the mono- or bi-phasic pattern, but if this patient follows a bi-phasic pattern, the associated bi-phasic design matrix is known. That is, for each patient regardless phasicity, the true mono- and bi-phasic design matrices are known; the only unknown quantity is the phasicity, η_i. Assuming the prior distributions P(σ²) = (σ²)⁻¹, P(λ) = Beta(1, 1) = 1, P(Σ_S) = |Σ_S|^−(2+1)/2 and P(Σ_B) = |Σ_B|^−(4+1)/2, where λ denotes the proportion of bi-phasic patients and η_i ∼ Ber(λ), the posterior distribution is: (4) Here η = (η₁, …, η_N) is the missing indicator variable for phasicity, and and denote the individual mono- and bi-phasic design matrices, respectively: (5) (6)

However, the bi-phasic design matrix for subject i, , is not known. The estimation of this bi-phasic change point is a well-known problem in statistics, mathematics, and computer science with many applications in other fields; several methods for addressing this question have been suggested [24, 25]. Here we employed a Bayesian formulation of the change point problem, suggested by Carlin et al. [26]. For a particular patient i with M_i observations, there are M_i − 1 possible change points. The cases in which the bi-phasic transition point occurs before the first observation or after the last observation are ignored, because in such cases mono- and bi-phasic subjects are not distinguishable. Thus, for each patient i, there are M_i − 1 possible design matrices, for example: (7) (8)

For each corresponding design matrix , the probability associated with the j-th bi-phasic design matrix is denoted by π_ij. Assuming all patients comply with their clinic visit schedules, such that t₁₁ = t₂₁ = … = t_N1, …, and t_1M₁ = t_2M₂ = … = t_{NM_N}, and assuming an uninformative Dirichlet prior, Dir(α₁ = α₂ = … = α_{M_i−1} = 1), the posterior function is: (9) where ξ_ij is the unobserved indicator for the j^th bi-phasic design matrix for subject i, such that and ξ_ij ∼ Multinomial(π_ij), and ξ_i = (ξ_i0, …, ξ_{iM_i}) and ξ = (ξ₁, …, ξ_N); and π_ij is the probability that the j-th bi-phasic design matrix is selected for the i-th patient. Note that the inclusion of the Dirichlet prior results in a constant scaling factor, and hence it is not included in Eq (9).

Given the complexity of the model, we first used the Expectation Maximization (EM) algorithm to search for the mode of the posterior distribution, which was then used as the starting value for the Gibbs sampler. To implement the EM algorithm and to obtain Empirical Bayes estimators, we utilized the procedures derived by Verbeke and Lesaffre [5] and Xu and Hedeker [6]. Following the notation used in Xu and Hedeker, the Empirical Bayes estimators for individual regression parameters and the covariance matrices are given by (10) for a given set of S, B, Σ_S, and Σ_B.

In the expectation step, the quantity (11) is calculated. ζ_ij denotes the probability of the j-th bi-phasic matrix for the i-th patient is selected to be bi-phasic design matrix, and ζ_i = (ζ_i1, ζ_i2…ζ_{iM_i}). Similarly, the expected value for η_i can be calculated as (12) where (13) (14) and where Eqs (13) and (14) are the posteriors for the mono- and bi-phase pieces. However, in practice, given the added model complexity of the bi-phasic model compared to the mono-phasic model, the bi-phasic model has a larger likelihood than the single phasic model, resulting in most patients being classified as bi-phasic. To compensate for the difference in model complexity, we instead of using Eqs (13) and (14) to differentiate each patient’s phasicity, we used the negative Bayesian Information Criteria (BIC) for the mono- and bi-phasic models, respectively: (15) (16) The additional terms, −2log(M_i) and −4log(M_i) are constant factors, which can be seen as prior odds for distinguishing between mono- and bi-phasic models for each patient. Because they are constants, they only result in a proportional change in the posterior function, Eq (9). Similar methods of using BIC to determine the posterior model probabilities have been implemented and discussed by Kass and Raftery [27]. BIC corrects for the improvement in fitting associated with increasing model complexity and BIC has been shown to be a consistent model selector due to its quickly increasing penalty as a function of the sample size [28–30]. We also investigated AIC as an alternative penalty function. However, as the sample size increases, the penalty becomes too weak, resulting in mono-phasic patients being misclassified as bi-phasic patients. In addition, the use of BIC for model selection is analogous to the use of DIC (Deviance Information Criterion) in Bayesian mixture model [31].

The maximization consists of the following steps: (17) where .

Updating the probability weight for π_ij, if all patients adhere to the visit schedule, is straightforward by averaging ζ_ij. However, in practice, patients often miss scheduled visits entirely or have unscheduled visits. Such departure from the trial design creates misalignments in patients’ observation intervals. For instance, two bi-phasic patients, i and i′, are identical except for their j + 1^th visits. Patient i’s j + 1^th visit is 1 week later than the scheduled time and patient i′ is on time. Because of this difference in visit time, ζ_ij and ζ_i′j can no longer be simply averaged to update π_ij. Instead, to account for misalignments, the transition probability π_ij associated with bi-phasic transition design matrix Q_ij needs to adapt to each patient’s actual visit time. The phasic transition density as a function of time is (18) , where T denotes the maximum follow-up time for all patients. The denominator is the weighted sum of time interval lengths in which phasic transitions occur; this denominator serves as the normalizing factor to ensure θ(t) integrates to 1. Eq (14) results in a numeric stepwise function specifying the transition density at a given t based on the E-step. The weight for each interval for each patient can be updated by integrating over its corresponding interval: (19)

The starting values for the EM algorithm are obtained using an ad hoc grid search procedure as outlined in S1 File.

The estimated parameter mode from the EM algorithm are used as the starting values for the Gibbs’ sampler to simulate the posterior distributions of the parameters from the model specified in Eq (9), as outlined below:

Calculate ζ_i for each patient, using Eq (11).
For each patient, draw a ξ_i vector from a multinomial distribution with a parameter vector ζ_i, and obtain the corresponding bi-phasic design matrix , such that ξ_ij = 1.
Calculate z_i based on the mono-phasic design matrix and the bi-phasic design matrix from step (2).
Draw η_i from a Bernoulli distribution with parameter z_i for each patient.
Update θ(t) using Eq (18), with ζ_ij replaced by ξ_ij and z_i replaced by η_i.
Draw a vector π_i from a Dirichlet distribution with a parameter vector , for each patient.
Draw λ from a Beta distribution with parameters (, ).
Sampling s_i and b_i: (20) where and and and are from Eq (12). The design matrix for the bi-phasic is drawn from step (2).
Sampling S and B from s_i and b_i: (21)
Sampling Σ_S and Σ_B: (22)
Sampling σ²: (23)

The proposed Gibbs’ sampler is similar to the methods used for variable selection via Gibbs sampling proposed by George and McCulloch [32], and, Carlin and Chib [33]. The only difference between our methods and the model by Carlin and Chib is the lack of pseudoprior for the transition probability between mono- and bi-phasic models. An important consequence of this is that, as noted in Carlin and Chib’s publication, “it is tempting to skip the generation of actual pseudoprior values … although seemingly reasonable, such an algorithm is clearly not a Gibbs sampler in the strict sense, since the nodes visited are determined by the current value in the realized Markov Chain.” However, in practice, as shown by our simulation studies, this heuristic Gibbs sampler performs well. Other methods such as reversible-jump MCMC may also be used to sample the posterior distributions [34]; however given the simplicity of the Gibbs sampler and its close relationship with the EM algorithm, we decided to use Gibbs’ sampler to implement the MCMC chain.

Results

Simulation results

We designed three simulation studies to test the model’s abilities to categorize patients and to estimate associated parameters, Fig 1. All three scenarios have the same population-level regression parameters, as shown in Table 1; the differences between the three scenarios lie in the covariance matrices specifying between-patient variability. The first simulation scenario assumes that there is no between-patient variability; for the second scenario, there is between-patient variability in the intercepts and slopes in both mono- and bi-phasic patients but no correlation among these parameters, i.e. all non-diagonal entries in Σ_s and Σ_b are zero. The third scenario assumes a correlation of 0.5 between the first and second slopes among the bi-phasic patients. In each scenario, we simulated N = 100 patients. The probability of being a bi-phasic patient is λ = 0.60. According to a hypothetical clinical protocol, patient data are collected every 21 days with 1 at baseline and 17 at follow-up visits for a total follow-up duration of 357 days. In this simulation study, the actual visit time may deviate within ± 5 days from the scheduled time. The true individual regression parameters are drawn from multivariate normal distributions with respective population parameters and covariance matrices, (S, Σ_S) or (B, Σ_B), depending on patients’ phasicity. For each simulated bi-phasic patient i, the phasic transition time occurs at .

Download:

Fig 1. Longitudinal trajectories for the simulated patients in the three scenarios.

Blue lines indicate mono-phasic patients’ trajectories, and red lines, bi-phasic patients’ trajectories. Vertical solid lines indicate the median time at which phasic transitions occur for the bi-phasic patients; vertical dashed lines indicate the 10^th% and 90^th% phasic transition time. All bi-phasic patients have the same phasic transition time in scenario one; hence, the dashed and solid lines coincide.

https://doi.org/10.1371/journal.pone.0180756.g001

Download:

Table 1. The true and the means of the estimated parameters in the three simulation scenarios.

https://doi.org/10.1371/journal.pone.0180756.t001

We first applied the EM algorithm to estimate the parameter values that maximize the marginal likelihood. The true and estimated parameters, excluding the covariance for the three scenarios, are shown in Table 1. In all three scenarios, the proposed model was able to provide parameter estimates that are close to the true parameter values, except for the bi-phasic proportion parameter λ, which is biased towards the mono-phasic model in scenarios 2 and 3. Estimating the covariance matrices for scenarios two and three is more challenging, as shown in Table 2. In particular, we found that the proposed model consistently over-estimates the variance term associated with the second intercept for the bi-phasic patients. Three possible causes for this over-estimation are 1) the bi-phasic design matrices must be estimated and misclassification of observations between the first and second phases may result in an enlarged variance term for the second intercept; 2) estimation of the second intercept requires projection back to time zero, and any uncertainty is magnified by this projection; and 3) the phasic transition time in our simulated data is distributed according the Gaussian ratio distribution, with heavy tails [35]; thus, bi-phasic patients with extreme transition times may not be classified correctly. In addition to parameter estimation, the proposed method performed well in classifying patients according to their phasicities, as shown in Table 3

Download:

Table 2. The true and the means of estimated covariance components in the three simulation scenarios.

https://doi.org/10.1371/journal.pone.0180756.t002

Download:

Table 3. Classification accuracy for the three scenarios.

https://doi.org/10.1371/journal.pone.0180756.t003

In addition to the three scenarios outlined above, we also performed sensitivity analyses to test the effects of variability in population-level intercepts and slopes, and , on the classification accuracy, Fig 2. For each of the three scenarios, we tested a grid of values for the bi-phasic first slope, B₁ (−0.45, …, −0.26), and the second slope, (−0.24, …, −0.05), centering around the mono-phasic slope S₁ = −0.25. For this sensitivity analysis, the second intercept for the bi-phasic patients was kept at values such that the population-level phasic transition times occurred in the middle of the time span of the trial (178 days). All other parameters, S₀, S₁, B₀, σ and λ, were kept at the values used in the previous three scenarios. The covariance matrices, if applicable, were also kept at the values used in the three scenarios. As expected, as the bi-phasic first and second slopes approached the value of the mono-phasic slope, the specificity diminished in all three scenarios such that more bi-phasic patients were misclassified as mono-phasic patients. Due to the strong penalty induced by the BIC correction in deciding on patients’ phasicities, the proposed model is biased toward the mono-phasic model. Sensitivity is close to 100% in all three scenarios; hence, the contour plots for sensitivity are not shown. In addition, we also investigated the effects of the numbers of observations per patient and the effects of the numbers of patients on the method’s ability to distinguish between mono- and bi-phasic patterns. As expected, as the number of observations per patient decreases, specificity decreases. Interestingly, the model is not very sensitive towards the total number of patients as shown in Fig 3.

Download:

Fig 2. Specificity as a function of true bi-phasic slopes.

True mono-phasic slope is kept at −0.25; bi-phasic first slopes vary between −0.45 and −0.26; bi-phasic second slopes vary between −0.24 and −0.05. Population-level mono-phasic slope and bi-phasic first intercepts are 90 and 91 respectively; the second slopes for bi-phasic patients are selected such that the population-level phase transition occurs at 178 days, which is in the middle of 357-day trial period. Each graph is generated based on the averages of 10 simulations. Sensitivity is omitted since it is at 100% for all given scenarios; please refer to Table 3.

https://doi.org/10.1371/journal.pone.0180756.g002

Download:

Fig 3. Specificity as a function of the numbers of observations per patient and the total numbers of patients per simulated data set.

Parameter values are identical to those used in the three scenarios. The numbers of observations per patient vary in the figure on the left; these observations are evenly distributed between day 0 and day 357. The total numbers of patients, N, vary in the figure on the right; the number of observations per patient is kept at 18.

https://doi.org/10.1371/journal.pone.0180756.g003

We also compared our model and its estimates with a standard mixture model package, Flexmix, a publicly available package in R [36, 37]. For each simulation, we ran the Flexmix package with and without providing the true design matrix, and true mixture identity as the initial values. The true design matrix groups data points from the same phase together for each subject; the true mixture identity provides the initial clustering of data points from the same phase across different patients together. The means of the parameters from 1,000 simulations are shown in Table 4. Because the Flexmix package is not designed to estimate the change point, the primary comparison of interest is to compare Flexmix’s ability to identify and estimate parameters associated with the three components of the mixture: single-phasic, bi-phasic first phase, and bi-phasic second phase. For scenario 1, in which is no between-subject variability, without providing both true design matrix and mixture identity, the mixture model was not able to estimate the parameters accurately, as compared to our model Table 1. Providing the true design matrix and mixture identity greatly improves the mixture model’s ability to estimate the parameters; however, this improvement is only limited to the case in which there is no between-subject variability. Once between-subject variability is introduced in scenarios 2 and 3, the mixture model was not able to estimate these parameters correctly.

Download:

Table 4. Comparison of estimates from a mixture model with and without a true design matrix and true cluster identity using the Flexmix package.

https://doi.org/10.1371/journal.pone.0180756.t004

In addition to obtaining the maximum likelihood parameter estimates, Gibbs’ sampling was implemented to obtain the posterior densities for the estimated parameters for the three scenarios. The key parameter of interest in this simulation study is the coverage probability for the proposed model. We simulated 1,000 independent data sets using identical parameter values for each scenario. The EM algorithm was first applied to maximize the likelihood; using the maximum likelihood parameter estimates as starting values, a Markov Chain Monte Carlo simulation was performed for each data set, with the number of iterations per simulation equal to 30,000. With samples generated from the posterior distributions, we constructed a 95% simultaneous rectangular credible region for each simulated data set, using the method outlined by Held [38, 39]. The coverage probability is calculated as the probabilities of the simultaneous credible regions covering all 8 parameters for scenario one and covering all 21 parameters for scenarios two and three, out of the 1,000 simulated data sets. The coverage probabilities are 80.1%, 71.5% and 65.9% for the three scenarios, respectively. The convergence of the Gibbs’ sampler is shown in S1 File.

We then performed detailed analyses to determine the parameter with the worst coverage probability in each scenario. For scenario one, the second intercept for the bi-phasic patients, , had the worst coverage probability of only 80.1%. Further analyses of scenarios two and three revealed that the variance component for the bi-phasic second intercept had the lowest coverage rate among all parameters; the 95% simultaneous credible regions for all 21 parameters were able to cover the covariance term for the second intercept only at rates of 82.5% and 80.1% for scenarios two and three, respectively. In this case, the covariance components obtained from the Gibbs’ sampler were consistently higher than the true covariance used in the simulated data. The same reasons used to explain the enlarged covariance structure for the EM results can be applied here.

Another parameter of particular interest is the correlation between the first and second slopes in scenario three. Focusing only on this parameter, our model was able to detect a correlation in 62.1% of the simulation runs; detection refers to 0 being excluded from the 95% credible region. Overall, the actual coverage probabilities from the proposed model are lower than the nominal probabilities. Model complexity appears to contribute to this poor coverage, as previous research has shown that even in the simple binomial case, the coverage probability rarely agrees with the nominal probability [40]. In addition, the parameter values, particularly the slopes, used in our simulation have considerable overlaps, which renders identifying patients’ phasicities difficult, and hence lowers the coverage probabilities.

Application results

To further demonstrate the utility of the proposed methods, we applied it to the M-protein data from the Velcade as Initial Standard Therapy in Multiple Myeloma: Assessment with Melphalan and Prednisone (VISTA) trial [41]. Briefly, the VISTA trial is a randomized, open-label phase III study, consisting of 682 patients with newly diagnosed, previously untreated, symptomatic, measurable multiple myeloma. In this study, patients were randomized to treatment with either melphalan and prednisone with (VMP cohort) or without (MP cohort) bortezomib (Velcade, Johnson & Johnson Pharmaceutical R&D and Millennium). Measurable disease was defined as the presence of quantifiable M-protein in serum or urine, or measurable soft-tissue or organ plasmacytomas. The longitudinal M-protein data from patients in the VISTA trial are shown in Fig 4.

Download:

Fig 4. Longitudinal trajectories for patients in the VISTA trial separated by treatment cohorts.

The mono-phasic (blue) and bi-phasic (red) lines indicate the population-mean trajectories based on the maximum likelihood estimates from the EM algorithm.

https://doi.org/10.1371/journal.pone.0180756.g004

The parameter estimates from our model revealed several interesting features associated with the M-protein dynamics, Table 5. First, the differences between the first and second slopes for the bi-phasic patients in both cohorts are striking. For the bi-phasic patients, the gradient of the first slope was lower than the second slopes in both cohorts, judging by the posterior credible regions. Second, more patients in the VMP cohort displayed bi-phasic trajectories than in the MP cohort. Third, the gradient of the bi-phasic first slope in the VMP cohort is lower than that of the bi-phasic first slope in the MP cohort. Fourth, the bi-phasic first intercepts are similar in both cohorts. Fifth, the long-term declines for the bi-phasic patients in both cohorts are similar. Sixth, in both cohorts, the intercepts for the mono-phasic patients tend to be substantially smaller than the first intercepts of the bi-phasic patients. Lastly, in the MP cohort, despite the large differences in the rates of initial declines, the long-term declines are very similar between the mono-phasic and the bi-phasic patients, as shown by the similarity in the estimates for S₁ and . Those observed differences in the M-protein dynamics between cohorts suggest that the tumor dynamics of multiple myeloma are highly complex.

Download:

Table 5. 95% simultaneous credible regions for the MP and VMP cohorts in the VISTA trial.

https://doi.org/10.1371/journal.pone.0180756.t005

Discussion

We have proposed a piecewise linear mixture random-effects model to investigate the extent of heterogeneity and time-varying functional relationships in longitudinal biomarker data. The combination of heterogeneity and a time-varying functional relationship is where the innovation in the proposed model lies. Our model assumes a simple yet robust piecewise linear functional form. The major advantage of this piecewise linear functional form over other more complex nonlinear functions is that the likelihood can be maximized analytically, using empirical Bayes estimators and standard expectation-maximization algorithms. No prior knowledge of the functional relationship other than the piecewise assumption is required; thus this method is particularly useful for initial exploratory analyses. In addition, the ease of interpretation of the parameter estimates is another advantage of the proposed model. Lastly, in the extreme case in which all patients are mono-phasic, the proposed model completely reduces to linear mixed-effects model.

One minor drawback of our approach is that for the bi-phasic patients, the proposed model produces a point of discontinuity between k_i and k_{i + 1} observations Eq (1). Nonlinear models, such as the broken-stick model, Bacon Watts model, and the polynomial model suggested by Matthews et al. offer potential solutions to this problem [42]; however, analytical solutions do not exist for those nonlinear functions. Another minor problem is that there is a small bias for the parameter denoting the proportion of bi-phasic patients, λ, in the EM algorithm and Gibbs’ sampler. Closer investigation reveals that this bias is not due to our proposed model; rather, it is an artifact of the data generation process for the simulation studies (see S1 File for detail). The simulated data for the bi-phasic patients are generated using a multivariate normal distribution. An unwanted implication of this generation mechanism is that the phasic transition time follows a Gaussian Ratio distribution, with heavy tails, such that for some bi-phasic patients the actual transition time may exceed the window of observation. This problem is particularly confounding in the case in which the first and second slopes are close in value to the slope of the true single-phasic slope parameter. As a result, such bi-phasic patients are indistinguishable from single-phasic patients. Thus, this small bias indicates that our proposed method was able to classify these bi-phasic patients “correctly” as single-phasic, based on the observed data. Fig A in S1 File shows a few such examples. The four figures in the bottom left corner of Fig A in S1 File are from true bi-phasic patients; however given the late phasic transitions and steep second slopes, our algorithm is likely to classify them as mono-phasic. This “misclassification” is due to the similarity of those patients’ trajectories to those mono-phasic patients rather than a systematic mistake in our algorithm. Lastly, the MCMC algorithm requires a large sample size to be implemented due to its model complexity—21 parameters in total. We recommend a sample size of at least 100 patients and with sufficient numbers of mono- and bi-phasic patients, greater than 40 each, to ensure that the number of patients is greater than the number of parameters.

When applying this algorithm to a data set from the VISTA trial, although from our analysis we have not found a significant correlation between phasicities and patients outcomes, the distinct mono- and bi-phasic trajectories may have significant medical implications warranting further investigation and validation. Those findings on the distinct treatment responses for patients randomized to the same treatment arm may help generate new hypotheses for improving patient prognosis and disease management.

From the prospective of clinical trial design, one interesting question is how to design a trial to maximize phasicity detection, if phasicity is important for patient management. Our linear framework may offer a simple approach to address these issues. In addition, our model can be extended beyond between-patient variability to include additional layers inside the hierarchy. For instance, patients with metastatic solid tumors or multiple tumors at different sites may demonstrate a large degree of similarity in terms of individual tumor trajectories, yet also exhibit kinetic differences depending on an individual tumor’s microenvironment. Modeling treatment responses in such scenarios would then require the incorporation of between-patient variability and within-patient/between-tumor variability into the model. Furthermore, our model can be extended to include multi-category and multi-phasic changes. The computational tractability problem associated with multi-phasic changes can be addressed from a practical point of view, such as restricting the number of minimum observations in each phase to be greater than 5 data points. These additional extensions can further enrich our proposed model.

Supporting information

S1 File. EM starting value search algorithm, additional EM examples, and assessment of Gibbs’ sampler convergence.

https://doi.org/10.1371/journal.pone.0180756.s001

(PDF)

Acknowledgments

This work was supported by the Dana-Farber Cancer Institute Physical Sciences Oncology Center (U54CA143798) and NIH grant R37 51164. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References

1. Laird NM, Ware JH. Random-effects models for longitudinal data. Biometrics. 1982; p. 963–974. pmid:7168798
- View Article
- PubMed/NCBI
- Google Scholar
2. Lindstrom MJ, Bates DM. Newton—Raphson and EM algorithms for linear mixed-effects models for repeated-measures data. Journal of the American Statistical Association. 1988;83(404):1014–1022.
- View Article
- Google Scholar
3. Fitzmaurice GM, Laird NM, Ware JH. Applied longitudinal analysis. vol. 998. John Wiley & Sons; 2012.
4. Belin TR, Rubin DB. The analysis of repeated-measures data on schizophrenic reaction times using mixture models. Statistics in Medicine. 1995;14(8):747–768. pmid:7644856
- View Article
- PubMed/NCBI
- Google Scholar
5. Verbeke G, Lesaffre E. A linear mixed-effects model with heterogeneity in the random-effects population. Journal of the American Statistical Association. 1996;91(433):217–221.
- View Article
- Google Scholar
6. Xu W, Hedeker D. A random-effects mixture model for classifying treatment response in longitudinal clinical trials. Journal of biopharmaceutical statistics. 2001;11(4):253–273. pmid:12018779
- View Article
- PubMed/NCBI
- Google Scholar
7. Jessica M, et al. A within-subject normal-mixture model with mixed-effects for analyzing heart rate variability. Journal of Biometrics & Biostatistics. 2012;.
- View Article
- Google Scholar
8. Verbeke G, Molenberghs G. Linear mixed models for longitudinal data. Springer; 2009.
9. Pakker NG, Notermans DW, De Boer RJ, Roos M, De Wolf F, Hill A, et al. Biphasic kinetics of peripheral blood T cells after triple combination therapy in HIV-1 infection: a composite of redistribution and proliferation. Nature medicine. 1998;4(2):208–214. pmid:9461195
- View Article
- PubMed/NCBI
- Google Scholar
10. Arnaout RA, Wodarz D, et al. HIV–1 dynamics revisited: biphasic decay by cytotoxic T lymphocyte killing? Proceedings of the Royal Society of London Series B: Biological Sciences. 2000;267(1450):1347–1354. pmid:10972131
- View Article
- PubMed/NCBI
- Google Scholar
11. Tsiang M, Rooney JF, Toole JJ, Gibbs CS. Biphasic clearance kinetics of hepatitis B virus from patients during adefovir dipivoxil therapy. Hepatology. 1999;29(6):1863–1869. pmid:10347131
- View Article
- PubMed/NCBI
- Google Scholar
12. Neumann AU, Lam NP, Dahari H, Gretch DR, Wiley TE, Layden TJ, et al. Hepatitis C viral dynamics in vivo and the antiviral efficacy of interferon-α therapy. Science. 1998;282(5386):103–107. pmid:9756471
- View Article
- PubMed/NCBI
- Google Scholar
13. Michor F, Hughes TP, Iwasa Y, Branford S, Shah NP, Sawyers CL, et al. Dynamics of chronic myeloid leukaemia. Nature. 2005;435(7046):1267–1270. pmid:15988530
- View Article
- PubMed/NCBI
- Google Scholar
14. Stein AM, Bottino D, Modur V, Branford S, Kaeda J, Goldman JM, et al. BCR–ABL Transcript Dynamics Support the Hypothesis That Leukemic Stem Cells Are Reduced during Imatinib Treatment. Clinical Cancer Research. 2011;17(21):6812–6821. pmid:21903771
- View Article
- PubMed/NCBI
- Google Scholar
15. Morrell CH, Pearson JD, Carter HB, Brant LJ. Estimating unknown transition times using a piecewise nonlinear mixed-effects model in men with prostate cancer. Journal of the American Statistical Association. 1995;90(429):45–53.
- View Article
- Google Scholar
16. Pinheiro JC, Bates DM. Mixed effects models in S and S-PLUS. Springer; 2000.
17. Naumova EN, Must A, Laird NM. Tutorial in biostatistics: evaluating the impact of ‘critical periods’ in longitudinal studies of growth using piecewise mixed effects models. International journal of epidemiology. 2001;30(6):1332–1341. pmid:11821342
- View Article
- PubMed/NCBI
- Google Scholar
18. Cudeck R, Klebe KJ. Multiphase mixed-effects models for repeated measures data. Psychological methods. 2002;7(1):41. pmid:11928890
- View Article
- PubMed/NCBI
- Google Scholar
19. Harring JR, Cudeck R, du Toit SH. Fitting partially nonlinear random coefficient models as SEMs. Multivariate Behavioral Research. 2006;41(4):579–596. pmid:26794919
- View Article
- PubMed/NCBI
- Google Scholar
20. Pauler DK, Laird NM. A mixture model for longitudinal data with application to assessment of noncompliance. Biometrics. 2000;56(2):464–472. pmid:10877305
- View Article
- PubMed/NCBI
- Google Scholar
21. Lu X, Huang, Y. Bayesian analysis of nonlinear mixed-effects mixture models for longitudinal data with heterogeneity and skewness. Statistics in medicine. 2014;.
22. Hammer SM, Vaida F, Bennett KK, Holohan MK, Sheiner L, Eron JJ, et al. Dual vs single protease inhibitor therapy following antiretroviral treatment failure: a randomized trial. Jama. 2002;288(2):169–180. pmid:12095381
- View Article
- PubMed/NCBI
- Google Scholar
23. Perelson AS, Essunger P, Cao Y, Vesanen M, Hurley A, Saksela K, et al. Decay characteristics of HIV-1-infected compartments during combination therapy. Expert Reviews in Molecular Medicine. 1997;.
24. Basseville M, Nikiforov IV, et al. Detection of abrupt changes: theory and application. vol. 104. Prentice Hall Englewood Cliffs; 1993.
25. Chen J, Gupta AK. Parametric statistical change point analysis: with applications to genetics, medicine, and finance. Springer; 2011.
26. Carlin BP, Gelfand AE, Smith AF. Hierarchical Bayesian analysis of changepoint problems. Applied statistics. 1992; p. 389–405.
27. Kass RE, Raftery AE. Bayes factors. Journal of the american statistical association. 1995;90(430):773–795.
- View Article
- Google Scholar
28. Nishii R, et al. Asymptotic properties of criteria for selection of variables in multiple regression. The Annals of Statistics. 1984;12(2):758–765.
- View Article
- Google Scholar
29. Hannan EJ, Quinn BG. The determination of the order of an autoregression. Journal of the Royal Statistical Society Series B (Methodological). 1979; p. 190–195.
- View Article
- Google Scholar
30. Shibata R. Selection of the order of an autoregressive model by Akaike’s information criterion. Biometrika. 1976;63(1):117–126.
- View Article
- Google Scholar
31. Celeux G, Forbes F, Robert CP, Titterington DM, et al. Deviance information criteria for missing data models. Bayesian analysis. 2006;1(4):651–673.
- View Article
- Google Scholar
32. George EI, McCulloch RE. Variable selection via Gibbs sampling. Journal of the American Statistical Association. 1993;88(423):881–889.
- View Article
- Google Scholar
33. Carlin BP, Chib S. Bayesian model choice via Markov chain Monte Carlo methods. Journal of the Royal Statistical Society Series B (Methodological). 1995; p. 473–484.
- View Article
- Google Scholar
34. Green PJ. Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika. 1995; p. 711–732.
- View Article
- Google Scholar
35. Hinkley DV. On the ratio of two correlated normal random variables. Biometrika. 1969;56(3):635–639.
- View Article
- Google Scholar
36. Leisch F. Flexmix: A general framework for finite mixture models and latent glass regression in R. Journal of Statistical Software. 2004;.
- View Article
- Google Scholar
37. R Core Team. R: A Language and Environment for Statistical Computing; 2014. Available from: http://www.R-project.org/.
38. Held L. Simultaneous posterior probability statements from Monte Carlo output. Journal of Computational and Graphical Statistics. 2004;13(1).
- View Article
- Google Scholar
39. Besag J, Green P, Higdon D, Mengersen K. Bayesian computation and stochastic systems. Statistical Science. 1995; p. 3–41.
- View Article
- Google Scholar
40. Agresti A, Coull BA. Approximate is better than “exact” for interval estimation of binomial proportions. The American Statistician. 1998;52(2):119–126.
- View Article
- Google Scholar
41. San Miguel JF, Schlag R, Khuageva NK, Dimopoulos MA, Shpilberg O, Kropff M, et al. Bortezomib plus melphalan and prednisone for initial treatment of multiple myeloma. New England Journal of Medicine. 2008;359(9):906–917. pmid:18753647
- View Article
- PubMed/NCBI
- Google Scholar
42. van den Hout A, Muniz-Terrera G, Matthews FE. Smooth random change point models. Statistics in medicine. 2011;30(6):599–610. pmid:21337356
- View Article
- PubMed/NCBI
- Google Scholar

[ref1] 1. Laird NM, Ware JH. Random-effects models for longitudinal data. Biometrics. 1982; p. 963–974. pmid:7168798
View Article
PubMed/NCBI
Google Scholar

[2] View Article

[3] PubMed/NCBI

[4] Google Scholar

[ref2] 2. Lindstrom MJ, Bates DM. Newton—Raphson and EM algorithms for linear mixed-effects models for repeated-measures data. Journal of the American Statistical Association. 1988;83(404):1014–1022.
View Article
Google Scholar

[6] View Article

[7] Google Scholar

[ref3] 3. Fitzmaurice GM, Laird NM, Ware JH. Applied longitudinal analysis. vol. 998. John Wiley & Sons; 2012.

[ref4] 4. Belin TR, Rubin DB. The analysis of repeated-measures data on schizophrenic reaction times using mixture models. Statistics in Medicine. 1995;14(8):747–768. pmid:7644856
View Article
PubMed/NCBI
Google Scholar

[10] View Article

[11] PubMed/NCBI

[12] Google Scholar

[ref5] 5. Verbeke G, Lesaffre E. A linear mixed-effects model with heterogeneity in the random-effects population. Journal of the American Statistical Association. 1996;91(433):217–221.
View Article
Google Scholar

[14] View Article

[15] Google Scholar

[ref6] 6. Xu W, Hedeker D. A random-effects mixture model for classifying treatment response in longitudinal clinical trials. Journal of biopharmaceutical statistics. 2001;11(4):253–273. pmid:12018779
View Article
PubMed/NCBI
Google Scholar

[17] View Article

[18] PubMed/NCBI

[19] Google Scholar

[ref7] 7. Jessica M, et al. A within-subject normal-mixture model with mixed-effects for analyzing heart rate variability. Journal of Biometrics & Biostatistics. 2012;.
View Article
Google Scholar

[21] View Article

[22] Google Scholar

[ref8] 8. Verbeke G, Molenberghs G. Linear mixed models for longitudinal data. Springer; 2009.

[ref9] 9. Pakker NG, Notermans DW, De Boer RJ, Roos M, De Wolf F, Hill A, et al. Biphasic kinetics of peripheral blood T cells after triple combination therapy in HIV-1 infection: a composite of redistribution and proliferation. Nature medicine. 1998;4(2):208–214. pmid:9461195
View Article
PubMed/NCBI
Google Scholar

[25] View Article

[26] PubMed/NCBI

[27] Google Scholar

[ref10] 10. Arnaout RA, Wodarz D, et al. HIV–1 dynamics revisited: biphasic decay by cytotoxic T lymphocyte killing? Proceedings of the Royal Society of London Series B: Biological Sciences. 2000;267(1450):1347–1354. pmid:10972131
View Article
PubMed/NCBI
Google Scholar

[29] View Article

[30] PubMed/NCBI

[31] Google Scholar

[ref11] 11. Tsiang M, Rooney JF, Toole JJ, Gibbs CS. Biphasic clearance kinetics of hepatitis B virus from patients during adefovir dipivoxil therapy. Hepatology. 1999;29(6):1863–1869. pmid:10347131
View Article
PubMed/NCBI
Google Scholar

[33] View Article

[34] PubMed/NCBI

[35] Google Scholar

[ref12] 12. Neumann AU, Lam NP, Dahari H, Gretch DR, Wiley TE, Layden TJ, et al. Hepatitis C viral dynamics in vivo and the antiviral efficacy of interferon-α therapy. Science. 1998;282(5386):103–107. pmid:9756471
View Article
PubMed/NCBI
Google Scholar

[37] View Article

[38] PubMed/NCBI

[39] Google Scholar

[ref13] 13. Michor F, Hughes TP, Iwasa Y, Branford S, Shah NP, Sawyers CL, et al. Dynamics of chronic myeloid leukaemia. Nature. 2005;435(7046):1267–1270. pmid:15988530
View Article
PubMed/NCBI
Google Scholar

[41] View Article

[42] PubMed/NCBI

[43] Google Scholar

[ref14] 14. Stein AM, Bottino D, Modur V, Branford S, Kaeda J, Goldman JM, et al. BCR–ABL Transcript Dynamics Support the Hypothesis That Leukemic Stem Cells Are Reduced during Imatinib Treatment. Clinical Cancer Research. 2011;17(21):6812–6821. pmid:21903771
View Article
PubMed/NCBI
Google Scholar

[45] View Article

[46] PubMed/NCBI

[47] Google Scholar

[ref15] 15. Morrell CH, Pearson JD, Carter HB, Brant LJ. Estimating unknown transition times using a piecewise nonlinear mixed-effects model in men with prostate cancer. Journal of the American Statistical Association. 1995;90(429):45–53.
View Article
Google Scholar

[49] View Article

[50] Google Scholar

[ref16] 16. Pinheiro JC, Bates DM. Mixed effects models in S and S-PLUS. Springer; 2000.

[ref17] 17. Naumova EN, Must A, Laird NM. Tutorial in biostatistics: evaluating the impact of ‘critical periods’ in longitudinal studies of growth using piecewise mixed effects models. International journal of epidemiology. 2001;30(6):1332–1341. pmid:11821342
View Article
PubMed/NCBI
Google Scholar

[53] View Article

[54] PubMed/NCBI

[55] Google Scholar

[ref18] 18. Cudeck R, Klebe KJ. Multiphase mixed-effects models for repeated measures data. Psychological methods. 2002;7(1):41. pmid:11928890
View Article
PubMed/NCBI
Google Scholar

[57] View Article

[58] PubMed/NCBI

[59] Google Scholar

[ref19] 19. Harring JR, Cudeck R, du Toit SH. Fitting partially nonlinear random coefficient models as SEMs. Multivariate Behavioral Research. 2006;41(4):579–596. pmid:26794919
View Article
PubMed/NCBI
Google Scholar

[61] View Article

[62] PubMed/NCBI

[63] Google Scholar

[ref20] 20. Pauler DK, Laird NM. A mixture model for longitudinal data with application to assessment of noncompliance. Biometrics. 2000;56(2):464–472. pmid:10877305
View Article
PubMed/NCBI
Google Scholar

[65] View Article

[66] PubMed/NCBI

[67] Google Scholar

[ref21] 21. Lu X, Huang, Y. Bayesian analysis of nonlinear mixed-effects mixture models for longitudinal data with heterogeneity and skewness. Statistics in medicine. 2014;.

[ref22] 22. Hammer SM, Vaida F, Bennett KK, Holohan MK, Sheiner L, Eron JJ, et al. Dual vs single protease inhibitor therapy following antiretroviral treatment failure: a randomized trial. Jama. 2002;288(2):169–180. pmid:12095381
View Article
PubMed/NCBI
Google Scholar

[70] View Article

[71] PubMed/NCBI

[72] Google Scholar

[ref23] 23. Perelson AS, Essunger P, Cao Y, Vesanen M, Hurley A, Saksela K, et al. Decay characteristics of HIV-1-infected compartments during combination therapy. Expert Reviews in Molecular Medicine. 1997;.

[ref24] 24. Basseville M, Nikiforov IV, et al. Detection of abrupt changes: theory and application. vol. 104. Prentice Hall Englewood Cliffs; 1993.

[ref25] 25. Chen J, Gupta AK. Parametric statistical change point analysis: with applications to genetics, medicine, and finance. Springer; 2011.

[ref26] 26. Carlin BP, Gelfand AE, Smith AF. Hierarchical Bayesian analysis of changepoint problems. Applied statistics. 1992; p. 389–405.

[ref27] 27. Kass RE, Raftery AE. Bayes factors. Journal of the american statistical association. 1995;90(430):773–795.
View Article
Google Scholar

[78] View Article

[79] Google Scholar

[ref28] 28. Nishii R, et al. Asymptotic properties of criteria for selection of variables in multiple regression. The Annals of Statistics. 1984;12(2):758–765.
View Article
Google Scholar

[81] View Article

[82] Google Scholar

[ref29] 29. Hannan EJ, Quinn BG. The determination of the order of an autoregression. Journal of the Royal Statistical Society Series B (Methodological). 1979; p. 190–195.
View Article
Google Scholar

[84] View Article

[85] Google Scholar

[ref30] 30. Shibata R. Selection of the order of an autoregressive model by Akaike’s information criterion. Biometrika. 1976;63(1):117–126.
View Article
Google Scholar

[87] View Article

[88] Google Scholar

[ref31] 31. Celeux G, Forbes F, Robert CP, Titterington DM, et al. Deviance information criteria for missing data models. Bayesian analysis. 2006;1(4):651–673.
View Article
Google Scholar

[90] View Article

[91] Google Scholar

[ref32] 32. George EI, McCulloch RE. Variable selection via Gibbs sampling. Journal of the American Statistical Association. 1993;88(423):881–889.
View Article
Google Scholar

[93] View Article

[94] Google Scholar

[ref33] 33. Carlin BP, Chib S. Bayesian model choice via Markov chain Monte Carlo methods. Journal of the Royal Statistical Society Series B (Methodological). 1995; p. 473–484.
View Article
Google Scholar

[96] View Article

[97] Google Scholar

[ref34] 34. Green PJ. Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika. 1995; p. 711–732.
View Article
Google Scholar

[99] View Article

[100] Google Scholar

[ref35] 35. Hinkley DV. On the ratio of two correlated normal random variables. Biometrika. 1969;56(3):635–639.
View Article
Google Scholar

[102] View Article

[103] Google Scholar

[ref36] 36. Leisch F. Flexmix: A general framework for finite mixture models and latent glass regression in R. Journal of Statistical Software. 2004;.
View Article
Google Scholar

[105] View Article

[106] Google Scholar

[ref37] 37. R Core Team. R: A Language and Environment for Statistical Computing; 2014. Available from: http://www.R-project.org/.

[ref38] 38. Held L. Simultaneous posterior probability statements from Monte Carlo output. Journal of Computational and Graphical Statistics. 2004;13(1).
View Article
Google Scholar

[109] View Article

[110] Google Scholar

[ref39] 39. Besag J, Green P, Higdon D, Mengersen K. Bayesian computation and stochastic systems. Statistical Science. 1995; p. 3–41.
View Article
Google Scholar

[112] View Article

[113] Google Scholar

[ref40] 40. Agresti A, Coull BA. Approximate is better than “exact” for interval estimation of binomial proportions. The American Statistician. 1998;52(2):119–126.
View Article
Google Scholar

[115] View Article

[116] Google Scholar

[ref41] 41. San Miguel JF, Schlag R, Khuageva NK, Dimopoulos MA, Shpilberg O, Kropff M, et al. Bortezomib plus melphalan and prednisone for initial treatment of multiple myeloma. New England Journal of Medicine. 2008;359(9):906–917. pmid:18753647
View Article
PubMed/NCBI
Google Scholar

[118] View Article

[119] PubMed/NCBI

[120] Google Scholar

[ref42] 42. van den Hout A, Muniz-Terrera G, Matthews FE. Smooth random change point models. Statistics in medicine. 2011;30(6):599–610. pmid:21337356
View Article
PubMed/NCBI
Google Scholar

[122] View Article

[123] PubMed/NCBI

[124] Google Scholar

Figures

Abstract

Introduction

Materials and methods

Results

Simulation results

Application results

Discussion

Supporting information

S1 File. EM starting value search algorithm, additional EM examples, and assessment of Gibbs’ sampler convergence.

Acknowledgments

References