As circulating tumor DNA (ctDNA) levels can reflect disease progression, achieving a comprehensive understanding of the temporal evolution of ctDNA is key to informing clinical decision making. However, temporal changes can exhibit complex non-linear patterns and differ substantially across patients. Additionally, patient characteristics and outcomes may impact temporal change. Thus, traditional statistical approaches may be inadequate in characterizing ctDNA evolution over time. In this proof-of-concept study, we propose utilizing a new approach using a hierarchical random effects cubic spline model, which is sufficiently flexible to capture complex temporal ctDNA patterns while supporting the integration of patient characteristics. To demonstrate the benefits of the approach, a retrospective cohort of non-small cell lung cancer patients who received anti-EGFR therapies was analyzed. Model results are presented graphically in the form of patient-level response patterns, where each combination of patient characteristics produces a unique pattern. Patients with various ages, levels of health status, as well as mortality status were contrasted, where results provide examples of how the model can further our conceptualization of ctDNA dynamics and demonstrates how results can be used in targeted, patient-centered, clinical decision-making.
Article
Longitudinal Assessment of Circulating Tumor DNA: A Proposed Statistical Framework
https://doi.org/10.21203/rs.3.rs-3788054/v1
This work is licensed under a CC BY 4.0 License
You are reading this latest preprint version
Plasma-based next generation sequencing for comprehensive genomic profiling of oncology patients with advanced solid tumors is recommended by numerous guidelines and expert consensus statements including those from the National Comprehensive Cancer Network, the American Society of Clinical Oncology, the European Society for Medical Oncology, and the International Association for the Study of Lung Cancer. Plasma-based ‘liquid biopsies’ have many clinical applications in advanced cancers including genotyping, identification of actionable biomarkers with targeted therapeutic options, early assessment for treatment efficacy, monitoring of response to treatment and identification of acquired resistance mechanisms. Liquid biopsies have also been applied in cancer screening and early detection, detection of molecular residual disease after curative intent and predicting patients at risk of recurrence.1,2 Additional studies suggest serial measurements of circulating tumor DNA (ctDNA) can provide insight into disease progression. For instance, Sanz-Garcia et al. discussed how temporal ctDNA can have widespread utility in disease management across cancer types, while McLaren et al. proposed serial ctDNA testing be used to inform radiotherapy dosing.3,4 In addition to illustrating how serial ctDNA testing can support disease management, these studies highlight the power of providing a custom-tailored patient-level approach to oncology, a topic central to goals of this manuscript.
Interpreting serial ctDNA results is currently challenging. ctDNA levels can vary substantially within and between patients, posing distinct obstacles when analyzing these types of data.4 Similarly, because temporal ctDNA patterns might be impacted by patient characteristics such as age, gender, etc. these factors should be accounted in analyses.5,6 To handle the complexities inherent to temporal ctDNA patterns while simultaneously considering how patient characteristics might be associated with these patterns, we propose the adoption of a hierarchical cubic spline random effects model (HCSREM)—an extension of a hierarchal linear mixed model. While the model has been modified to meet the challenges posed by serial ctDNA data, the foundations of the model have been widely discussed and applied.7–17. As these prior works point out, this family of models has many advantages in comparison to traditional longitudinal analyses such as use of change scores, calculating ratios between baseline measures and future time points, or response profile analysis. In summary, the main benefits of HCSREM are: (1) The ability to capture biomarker evolution more accurately, (2) patient characteristics can be directly incorporated, (3) time points between measures need not be equally spaced and missing data is less problematic, (4) sample size can be preserved, that is, instead reducing sample size by performing sub-group analyses, a covariate that captures different strata can be integrated, and (5) longitudinal projections of existing, and new patients (those that do not exist in the data) can be generated.
To the best of our knowledge, HCSREM or similar models have not been utilized in analyzing complex longitudinal genomic data. As such, the main purpose of this paper is to discuss how to apply this model to serial ctDNAand to demonstrate the benefits of this framework by way of a specific application. Consequently, our results are for a proof-of-concept study to showcase the methodology and are not intended as a formal analysis to explore a predetermined clinical hypothesis.
General Approach
In this study the details of the HCSREM are discussed by demonstrating the model’s utility through analysis of a retrospective real-world cohort of patients who were diagnosed with advanced non-small cell lung cancer (NSCLC). Although the proposed framework can be applied to an assortment of longitudinal biomarkers, the biomarker of interest is ctDNA level as measured by the maximum variant allele fraction of all somatic variants detected through liquid biopsy. Since a major advantage of the method is the ability to incorporate patient information, several relevant covariates were considered. Finally, to enhance interpretation, model results are displayed graphically in the form of estimated longitudinal projections, each based upon a patient’s set of distinct traits. Within this process, patient-level projections are directly compared where comparisons are enhanced by velocity plots which are defined subsequently.
Data Source and Patient Cohort
The cohort used to illustrate the utility of the methodology is based on observational data and was sourced from the GuardantINFORM anonymized clinical-genomic database, which includes structured commercial payer claims collected from inpatient and outpatient facilities in both academic and community settings. Due to the observational nature of the data, the methods in this study adhered to the Strengthening the Reporting of Observational studies in Epidemiology (STROBE) guidelines where applicable.18 The GuardantINFORM database is fully deidentified and complies with Sections 164.514 (a)–(b)1ii of the US Health Insurance Portability and Accountability Act (HIPAA) regarding the determination and documentation of statistically deidentified data. The generation of de-identified data sets by Guardant Health for research purposes was approved by the Advarra Institutional Review Board; patient identity protection was maintained throughout the study in a de-identified database. Retrospective analysis of de-identified data is approved by Advarra IRB number Pro00034566. As the data is de-identified, obtaining patient consent was waived by the approving ethics committee.
Patients selected for the cohort were diagnosed with advanced non-small cell lung cancer (NSCLC) and had at least three Guardant360 (G360) liquid biopsy tests in the US between June 1st, 2014 and June 30th, 2023. Only patients receiving targeted therapies for EGFR mutations were included, with the following therapies considered: osimertinib, afatinib, dacomitinib, erlotinib, gefitinib, and amivantamab. All patients were required to retain at least three blood samples while on a specific line of anti-EGFR therapy, or within 30 days prior to line of therapy initiation and 30 days post end line of therapy. Patients whose first G360 test on the line of therapy was more than 120 days after the start of line of therapy were excluded. For patients with multiple lines of therapy meeting these criteria, the earliest line of therapy was selected for study inclusion. Finally, patients with suspected germline mutations were removed from the cohort.
Response Variable and Study Covariates
The response variable, ctDNA measurements captured over time, is reported as a percentage. In instances where samples contained ctDNA levels below the assay’s limit of detection, values were replaced with ctDNA levels of 0.04%, the lowest value in the cohort and consistent with the limit of detection of the test. All covariates except mortality were captured at baseline where the baseline period is defined as six months prior to the index date, i.e., the date of the patient’s first G360 test. Baseline covariates include age (in years), line of anti-EGFR therapy, smoking status (yes/no), gender (female/male), and the Van Walraven Elixhauser Comorbidity (ELIX) score specific to lung cancer patients (expressed as a weighted measure across multiple common comorbidities19). As the cohort is based on real-world-data, it is not possible to directly align treatment start date with the patient’s first G360 as is achievable in a prospective study. Therefore, days between the first G360 test and start of treatment was added as a covariate to serve as a statistical control and was set to zero days in the analysis to mimic a post treatment scenario. Patient mortality captured as alive vs deceased within the study timeframe was also included.
Proposed Statistical Model
In this section we provide an overview of the mathematical details of the HCSREM. As previously stated, this model was selected as it is malleable enough to capture variable nonlinear trends, and it allows for the direct incorporation of patient characteristics in the form of covariates. In addition to these properties, this model can provide a unique corresponding temporal ctDNA pattern for each combination of covariate values. It is the ability to provide this type of patient-specific information that makes this methodology attractive in targeted oncology efforts.
The model is partitioned into first- and second-level equations, which create the hierarchical structure. The first-level equation assumes the form of a truncated cubic spline and captures how a particular patient’s ctDNA levels change over time (see Eq. (1)). At a high level this is achieved by creating a function that is split into various pieces that traverse the abscissa. Within each piece, a cubic polynomial is used to fit the data where the ends of the consecutive cubic polynomials are connected by knots. Knot location, as well as the number of knots can be strategically devised based on data inspection, though “automated” methods for determining knot quantity and placement exist. Ultimately, a cubic spline model combines the separate pieces to form a single uniform function to represent the data.20
$${Y}_{ij}={\pi }_{0i}+{\pi }_{1i}{t}_{ij}+{\pi }_{2i}{t}_{ij}^{2}+{\pi }_{3i}{t}_{ij}^{3}+\sum _{k=1}^{K}{\pi }_{(k+3)i}{\left({t}_{ij}-{ϵ}_{k}\right)}_{+}^{3}+{\epsilon }_{ij}$$
1
where\({\left(t-ϵ\right)}_{+}=\left\{\begin{array}{cc}0& if t\le ϵ\\ t-ϵ& if t>ϵ\end{array}\right\}\)
In Eq. (1), ctDNA measurements (or a transformation thereof) captured over time are represented by the \({Y}_{ij}{\prime }s\), where \(i\) is used to index patients and \(j\) indexes the measurement occasion. Timepoints captured within the patient are given by \({t}_{ij}\), \(ϵ\) is the value of the \({k}^{th}\) knot, the \({\pi }_{ri}{\prime }s\) are the \(r\) response parameters, where each i.e. \({\pi }_{0i}\), \({\pi }_{1i}\), \(\cdots , {\pi }_{(k+3)i}\) varies across patients i.e. the random effects, and, \({\epsilon }_{ij}\) is the error term and is assumed to be normally distributed with a mean of 0 and variance \({\sigma }^{2}\). The response parameters are especially important as they collectively govern the shape of each patient’s unique longitudinal ctDNA trajectory and serve to bridge the first and second-level equations.
The significance of the second-level equations is they contain information about individual patient characteristics and associate these characteristics with the response parameters themselves. The second-level equations are given below.
$${\pi }_{ri}={\beta }_{r0}+\sum _{c=1}^{{C}_{r}}{\beta }_{rc}{X}_{ci}+{e}_{ri}$$
2
where, \({X}_{ci}\) represents a desired patient characteristic, \({\beta }_{rc}\) captures the linear relationship between the response parameter and the patient characteristic, \({\beta }_{r0}\) is the intercept for each corresponding \({\pi }_{ri}\), and \({e}_{ri}\) represents a random component and is assumed to adhere to the following multivariate normal distribution:
When the model contains covariates, it is referred to as a conditional model, otherwise it is an unconditional model. The unconditional model provides results at the cohort level and the conditional model is responsible for producing patient-level results.
A final, yet important element of the proposed methodology comes in the form of velocity plots, which are useful when examining the direction and speed in which ctDNA levels change at a given point in time i.e. the instantaneous rate of change (IRC) is of interest. Each model generated patient trajectory has at its heart, a cubic spline. An advantageous property of cubic splines is they are twice differentiable, thus, the IRC at a given time point can be calculated.21 In the case of the adopted spline model, this amounts to taking the first derivative of Eq. (1) with respect to time resulting in:
$$\frac{{dY}_{ij}}{d{t}_{ij}}={\pi }_{1i}+2{\pi }_{2i}{t}_{ij}+3{\pi }_{3i}{t}_{ij}^{2}+\sum _{k=1}^{K}{3\pi }_{(k+3)i}{\left({t}_{ij}-{ϵ}_{k}\right)}_{+}^{2}$$
3
The value of the IRC is given by the slope of the line tangent to the patient trajectory, where positive values correspond to an increasing IRC, negative values to decrease, and IRC values of zero indicate either a peak or trough was reached, or that the trajectory is flat. The further the IRC value is from zero, the more extreme the rate of change is.
STATISTICAL ANALYSIS and RESULTS
Data were extracted using SAS software package 9.4 (SAS Institute, Cary, NC, USA) and all statistical analysis for the HCSREM was performed using R version 4.1.3. A total of 400 patients with advanced NSCLC were identified from the GuardantINFORM database as having at least three G360 tests. Seventy-three patients were excluded as their first test was more than 120 days after therapy initiation and five were excluded due to germline mutations. Of the remaining patients, 163 received anti-EGFR therapy with a total of 561 ctDNA longitudinal measurements, where these 163 patients defined the cohort used in the analysis. The average age of these patients was 62 years, 66% of them were females, average line of anti-EGFR therapy was 1 and the average time between G360 test and treatment initiation was 0 days (range − 115 days to 30 days) (Table 1).
|
Characteristics (Total N = 163) |
N/Mean |
%/Standard Deviation |
|---|---|---|
|
Age (years) |
61.18 |
10.88 |
|
Female |
108 |
66% |
|
ELIX score |
1.89 |
1.86 |
|
Current or prior smoker |
123 |
75% |
|
Line of anti-EGFR therapy |
1.44 |
0.99 |
|
Time between G360 test and treatment initiation (days) |
0.29 |
31.98 |
|
ctDNA (%)* |
5.66 |
10.59 |
|
Deceased at the end of study period |
55 |
33% |
| *ctDNA value was extracted from each test and summarized, thus includes multiple ctDNA values for each patient. | ||
[Insert Table 1]
To meet model assumptions, ctDNA levels, expressed as a percentage, were transformed into logits. Histograms in Fig. 1 show that logit transformation of ctDNA level alleviated the extreme skewness of the raw data. Likewise, the logit transformed spaghetti plot accentuates the complexity and variability of ctDNA level over time, both within and between patients.
[Insert Fig. 1]
In our example we concentrate on conducting an exploratory investigation of the data. To begin, an unconditional model was fit to the transformed data using knots set at 50, 125, 250, 500, 750, 1000, and 1250 days respectively. To ensure consistency, other knot orientations were explored, although different orientations did little to alter results. Results are presented graphically as spline model parameter estimates are difficult to interpret although parameter estimates and related output are provided in the supplemental information for reference.22 The graphical manifestation of the unconditional model, referred to as a response pattern, is presented in Fig. 2.
[Insert Fig. 2]
The response pattern (given in black) suggests ctDNA levels drop substantially between the first G360 test and 30 days, then rise rapidly until 150 days, at which point ctDNA levels dip slightly and rise again at around 300 days, although at a less extreme rate. Additionally, from 550 days to 1000 days ctDNA levels drop, and then rise again from 1000 to 1600 days. The corresponding 95% confidence band expand over time as the number of datapoints decreases. The flexibility built into the unconditional model revealed details hidden within the data that simpler models would not detect. Despite this, the unconditional model only estimates the response pattern for the cohort and does not account for the contingency that patients with different characteristics may exhibit different response patterns. To assess the impact of incorporating patient characteristics, a conditional model that incorporated all baseline covariates was fit to the data. As is typical in hierarchical models, all numerical covaries were centered about their respective means (Hofmann and Gavin (1998)). Figure 3 shows how baseline age and health status, as measured by the ELIX score, impact response patterns in female non-smokers receiving their first line of EGFR-TKI treatment. Results are separated by patients who are alive vs deceased. As data becomes sparse after 400 days, we examine the first 400 days only.
[Insert Fig. 3]
Examples presented above reveal that patients with different characteristics have different response patterns. In the top-left panel, the response curve for a 30 and 80-year-old with average ELIX scores are contrasted. These results suggest 80-year-old patients did not exhibit initial post treatment drop in ctDNA levels in comparison with 30-year-old patients who demonstrated a rapid decrease followed by a rapid increase. The top-middle panel indicates response patterns for patients with an average age and a maximum ELIX score of 13 appear to be quite different compared to the same patient with a minimum ELIX score of 0, implying patients with many comorbidities exhibited a delayed treatment response. In the top-right panel, response patterns are displayed for older patients with high comorbidity burden and younger patients who are otherwise healthy, illustrating how the age/health status combination amplifies the disparity in response patterns. Though not shown, beyond 400 days, a decreasing trend in ctDNA values is observed for patients who remained alive at the end of the study while the trend increases for patients who died before study end.
To focus on the response pattern’s behavior, velocity plots that display the IRC for a corresponding response pattern were generated (Fig. 4).
[Insert Fig. 4]
In general, information presented in a velocity plot can be gleaned from the response patterns themselves, but the differences in the response patterns are accentuated when examining them through an IRC lens. Thus, comparing velocity plots can provide additional clues as to where response patterns are similar and where they diverge based on the IRC value. Another advantage of utilizing velocity plots occurs when baseline values between response patterns are dissimilar and therefore differences between response patterns may be due to the fact that biomarker values were different at the onset. In these instances, using velocity plots to make comparisons may be more appropriate as the IRC is invariant to the biomarker’s baseline value.
Interpreting a velocity plot is relatively straightforward. To demonstrate, we focus on the far-left panels. In the first 100 days, velocity plots for 80-year-old alive and deceased patients (red curves) exhibited different patterns. For survivors, the IRC was initially positive but slowed to zero around 20 days (indicating a peak in the corresponding response curve as referenced by the dashed line), and then decreased, where the fastest rate of decrease (-0.026 logits per day) occurred around 43 days. Beyond 43 days the IRC continued to decrease and remained relatively flat past 100 days. In contrast, the velocity plot for 80-year-old deceased patients displays a nearly opposite pattern.
The purpose of this paper is to introduce a statistical methodology that can accommodate the analysis of complex longitudinal genomic data. In our example we applied the methodology to observational data and used it for exploratory purposes. However, this methodology is applicable in different data settings and can be used for hypothesis generation, statistical inference, and patient monitoring. In this section we discuss these applications while considering their limitations and conclude with discussing model modifications and sample size requirements.
In observational data settings, as was the case in our example, the proposed method is often leveraged as an exploratory tool. When data exploration is the aim, the 95% confidence bands fail to retain their traditional inferential meaning and instead are used as ‘guidelines’ in identifying differences in response patterns. Since thousands of response patterns are available (there are 58,464 possible response patterns in our example), exploring results can be daunting. To alleviate this issue, RShiny, Excel, or a similar platform can be used to create an interactive tool that generates a response pattern for each covariate combination, allowing for the interrogation of a large volume of results. During this process, caution should be taken in extrapolating beyond the data range or investigating nonsensical covariate combinations. Many observations may be made in the data exploration process. For instance, results may support predefined conjecture, or may reveal unanticipated findings prompting further investigation in a follow-up study designed to evaluate the finding. For instance, results presented in Fig. 3 may motivate a study that assess whether older patients with more comorbidities exhibit a delayed treatment response in comparison to the response of younger healthier patients.
This paper has focused on analyzing observational data, but this framework can be applied to representative cohorts as well. If statistical inference is the goal, since the potential exists to generate and compare numerous response patterns, the number of comparisons should be minimized, based on a priori hypotheses, and common considerations such as controlling for type-I error should be made. Hypotheses may include comparing response patterns between patients with pre-determined sets of covariate values (where other study covariates can be used as statistical controls) but can also include hypothesizing about the nature of the relationship between response pattern behavior and the covariate values themselves.
Patient monitoring is another application of our proposed framework. The general idea is each response pattern is a reasonable portrayal of a patient as described by his or her own unique set of characteristics—and—in this way—the same response pattern can serve as a reference for a new patient that shares these characteristics. Additionally, if survival status (deceased or not) is incorporated into the model, a reference response pattern for survivors and non-survivors can be created. Thus, if the response pattern of a new patient is consistent with that of a survivor, intervention is unnecessary, but if the response pattern mirrors that of a non-survivor, intervention may be required. Utilizing velocity plots to compare response patterns can further enhance this process-especially if baseline values between response patterns are dissimilar. To ensure reliable classification, such a monitoring system should undergo internal and external validation. Internal validation may be achieved by creating training and test datasets and then apply say k-fold cross-validation to assess classification accuracy. If an acceptable level of accuracy is achieved, external validation can be accomplished if new patients i.e. not involved in cross-validation are also classified with a high degree of accuracy.
Modifications can be made to improve model performance. For instance, the relationship between the response parameter and covariate(s) may be non-linear, and therefore, imposing a linear restriction as done in Eq. (2) will fail to correctly specify the model. In such cases, a scatterplot can be used to identify the correct relationship and appropriate adjustments can be made. Model efficiency (i.e. models with smaller standard errors) can also be improved via model reduction-especially when multicollinearity is present. Although model reduction can be an arduous task due to the complicated nature of these models, the reduction process can be guided by utilizing information criteria such as the Akaike information criteria (AIC) or Bayesian information criteria (BIC) which balance model fit with model simplicity. Finally, applicable when timepoints are fully captured and equally spaced (not the case in our example), is to account for the correlation structure of the repeated measures within the model. Possible correlation structures include the first order autoregressive, compound symmetry, and spatial power structures, where the AIC or BIC can facilitate the identification of an optimal structure. Finally, selecting an appropriate sample size is an important consideration when conducting an analysis. Minimum sample size recommendations for hierarchical models are around 100 patients with at least 3 measurements per patient, though these types of models have been fit with sample sizes as low as 22 individuals.23–25
Changes in ctDNA levels can fluctuate significantly over time from patient to patient, and the results can be difficult to interpret. In this paper we introduce a statistical framework that is flexible and therefore capable of capturing these complexities while accounting for a diverse set of patient traits. Furthermore, we apply the modeling to analyze an observational dataset consisting of patients with NSCLC who received anti-EGFR therapy. Analytic results are presented graphically and demonstrate the utility of the approach in acquiring a comprehensive understanding of how response patterns evolve and how different patient characteristics influence these evolutions. In our example we demonstrate how the method can be used as a “high powered” exploratory technique, although we have outlined many other applications. Regardless of the desired purpose, a major advantage of our proposed framework is the ability to generate patient-level results, where such results can add to our understanding of ctDNA dynamics and enhance our ability to integrate ctDNA into clinical decision-making.
SOFTWARE
R Core Team (2022). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
Acknowledgments
We thank Dr. Aaron Hardin, PhD for insights regarding ctDNA capture and the limit of detection of the G360 test.
Author Contributions
C.R.P. conceived of the presented model and computational framework, analyzed the data, presented the results, and authored the manuscript with critical input from all co-authors. J.L. created the cohort, aided in the development of the model and computational framework, provided analytic support, and authored the manuscript. C.W., L.D., L.B., & A.D. provided critical content knowledge expertise, guidance, and input, and authored the manuscript.
Data Availability Statement
The datasets generated during and/or analyzed during the current study are not publicly available and cannot be shared due to the use of a third-party healthcare claims database. Researchers interested in replicating our study or pursuing new research topics should contact Guardant Health (https://guardanthealth.com/products/biopharma-solutions/real-world-evidence/) directly.
Additional Information
Competing Interests: All manuscript authors are employed by Guardant Health. Each receive an annual salary, bonus, and stock options that are commensurate with the author’s job description, experience, and level of education.
- Rolfo, C. et al. Liquid Biopsy for Advanced NSCLC: A Consensus Statement From the International Association for the Study of Lung Cancer. Journal of Thoracic Oncology 16, 1647–1662 (2021).
- Pascual, J. et al. ESMO recommendations on the use of circulating tumour DNA assays for patients with cancer: a report from the ESMO Precision Medicine Working Group. Annals of Oncology 33, 750–768 (2022).
- McLaren, D. B. & Aitman, T. J. Redefining precision radiotherapy through liquid biopsy. Br J Cancer 129, 900–903 (2023).
- Sanz-Garcia, E., Zhao, E., Bratman, S. V. & Siu, L. L. Monitoring and adapting cancer treatment using circulating tumor DNA kinetics: Current research, opportunities, and challenges. Sci. Adv. 8, eabi8618 (2022).
- Ørntoft, M.-B. W. et al. Age-stratified reference intervals unlock the clinical potential of circulating cell-free DNA as a biomarker of poor outcome for healthy individuals and patients with colorectal cancer. Int J Cancer (2020) doi:10.1002/ijc.33434.
- Huang, R. S. P. et al. Circulating Cell-Free DNA Yield and Circulating-Tumor DNA Quantity from Liquid Biopsies of 12 139 Cancer Patients. Clinical Chemistry 67, 1554–1566 (2021).
- Raudenbush, S. W. & Bryk, A. S. Hierarchical Linear Models: Applications and Data Analysis Methods. (Sage Publ, Thousand Oaks, Calif., 2010).
- Singer, J. D. & Willett, J. B. Applied Longitudinal Data Analysis: Modeling Change and Event Occurrence. (Oxford University PressNew York, 2003). doi:10.1093/acprof:oso/9780195152968.001.0001.
- Hedeker, D. & Gibbons, R. D. Longitudinal Data Analysis. xx, 337 (Wiley-Interscience, Hoboken, NJ, US, 2006).
- Fitzmaurice, G. M., Laird, N. M. & Ware, J. H. Applied Longitudinal Analysis. (Wiley, Hoboken, 2011).
- Welham, S., Cullis, B., Gogel, B., Gilmour, A. & Thompson, R. Prediction in linear mixed models. Aust NZ J Stat 46, 325–347 (2004).
- Mackenzie, M. L., Donovan, C. R. & McArdle, B. H. Regression Spline Mixed Models: A Forestry Example. Journal of Agricultural, Biological, and Environmental Statistics 10, 394–410 (2005).
- Straube, J., Gorse, A.-D., PROOF Centre of Excellence Team, Huang, B. E. & Lê Cao, K.-A. A Linear Mixed Model Spline Framework for Analysing Time Course ‘Omics’ Data. PLoS ONE 10, e0134540 (2015).
- Pretz, C. R., Kozlowski, A. J., Chen, Y., Charlifue, S. & Heinemann, A. W. Trajectories of Life Satisfaction After Spinal Cord Injury. Archives of Physical Medicine and Rehabilitation 97, 1706–1713.e1 (2016).
- Grajeda, L. M. et al. Modelling subject-specific childhood growth using linear mixed-effect models with cubic regression splines. Emerg Themes Epidemiol 13, 1 (2016).
- Yu, Z. et al. Beyond t test and ANOVA: applications of mixed-effects models for more rigorous statistical analysis in neuroscience research. Neuron 110, 21–35 (2022).
- Janssen, J. M. et al. Longitudinal nonlinear mixed effects modeling of EGFR mutations in ctDNA as predictor of disease progression in treatment of EGFR -mutant non‐small cell lung cancer. Clinical Translational Sci 15, 1916–1925 (2022).
- von Elm, E. et al. The Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement: guidelines for reporting observational studies. J Clin Epidemiol 61, 344–349 (2008).
- van Walraven, C., Austin, P. C., Jennings, A., Quan, H. & Forster, A. J. A modification of the Elixhauser comorbidity measures into a point system for hospital death using administrative data. Med Care 47, 626–633 (2009).
- Gauthier, J., Wu, Q. V. & Gooley, T. A. Cubic splines to model relationships between continuous variables and outcomes: a guide for clinicians. Bone Marrow Transplant 55, 675–680 (2020).
- Zhou, S. & Wolfe, D. A. ON DERIVATIVE ESTIMATION IN SPLINE REGRESSION. Statistica Sinica 10, 93–108 (2000).
- Shepherd, B. E., Rebeiro, P. F., & Caribbean, Central and South America network for HIV epidemiology. Brief Report: Assessing and Interpreting the Association Between Continuous Covariates and Outcomes in Observational Studies of HIV Using Splines. J Acquir Immune Defic Syndr 74, e60–e63 (2017).
- Huttenlocher, J., Haight, W., Bryk, A., Seltzer, M. & Lyons, T. Early vocabulary growth: Relation to language input and gender. Developmental Psychology 27, 236–248 (1991).
- Fan, X. & Fan, X. Power of Latent Growth Modeling for Detecting Linear Growth: Number of Measurements and Comparison with Other Analytic Approaches. The Journal of Experimental Education 73, 121–139 (2005).
- Muthén, B. O. & Curran, P. J. General longitudinal modeling of individual differences in experimental designs: A latent variable framework for analysis and power estimation. Psychological Methods 2, 371–402 (1997).
Competing interest reported. All manuscript authors are employed by Guardant Health. Each receive an annual salary, bonus, and stock options that are commensurate with the author’s job description, experience, and level of education. However please note this is not a formal study but instead a proof-of-concept manuscript with the intent to showcase a novel statistical methodology.