Multiple-use prediction and calibration for all future values play an important role in many areas, including health and medical research. Simultaneous tolerance bands (STBs) can be employed for these purposes. This article studies the construction of exact STBs for multiple regression over any given rectangular covariate regions and for polynomial regression over any given covariate intervals. We first fill the gap in constructing the exact STBs for multiple regression over rectangular regions. We then propose a general form of STBs and develop a method for computing the critical constants for all STBs considered, including existing ones in the literature. A new STB is also proposed for both multiple and polynomial regressions. This new STB and several existing forms are then compared under the average shift (AS) criterion. Numerical results indicate that the new STB performs best under the AS criterion and is therefore recommended. In addition, a computationally efficient algorithm is presented. Real-life examples are provided for illustration.
Identifying the optimal dose-schedule regimen in early-phase oncology trials is complicated by competing risks, such as disease progression (DP) and dose-limiting toxicity (DLT). Many existing dose-finding methods fail to adequately address these events or accommodate varying administration schedules. We propose CR-EffTox, a Bayesian adaptive phase I/II trial design that jointly models time-to-event DLT and DP using cause-specific hazard functions. Besides, the model proposed allows for dynamic information borrowing to account for associations among dose-schedule regimes. To guide regimen selection, a novel satisfaction score derived from cause-specific survival curves is introduced to quantify the benefit-risk trade-off. The operating characteristics of the method are evaluated through extensive simulations. The method generally outperforms methods that ignore competing risks or information borrowing, substantially improving correct selection probability and enhancing patient safety by reducing allocation to suboptimal regimens.
Biomarkers measured from biofluid samples and imaging scans are important to aid diagnosis of diseases and track their progression over time, especially in neurodegenerative diseases such as Alzheimer's disease (AD). Combining retrospectively obtained biomarker data across multiple studies can increase statistical power, but existing biomarker data may be generated using different assay platforms, scanner types, or processing protocols by different studies, which may significantly affect the measurements of biomarkers and hence render it necessary to harmonize the data across the studies. An optimal way of biomarker data harmonization is to re-analyze all the biofluid samples or imaging scans together on a single platform in a central reference lab, but this is often not practical because of the substantial cost involved as well as the limited amount of biofluid samples available. A more practical solution is to prospectively design a bridging study by re-measuring a subset of biofluid samples or imaging scans from the studies in a reference lab to evaluate how biomarker values may be harmonized across studies. An important design question for such a bridging study is the size of the retrospectively collected biofluid samples or imaging scans that will be re-measured. We aim to address this question by conceptualizing a latent but true biomarker that underlies the observed versions of the biomarker measured across retrospective studies and proposing methods to determine the sample size of a bridging study for estimating the biological correlation of the true biomarker with a standard and validated clinical outcome. We also report bridging data of several analytes from cerebrospinal fluid (CSF) in a multi-center biomarker study of AD and demonstrate that a small proportion of the CSF samples may be sufficient to design a future bridging study for estimating the correlations of CSF biomarkers with a cognitive and functional outcome.
In most multi-arm randomized controlled trials, all participants need to be eligible to be randomized to all the arms. This means some participants are excluded from participating in multi-arm studies. We propose to relax this by considering a 'differential randomization' approach that allows the inclusion of participants who are eligible for some but not all arms, as in some innovative trial approaches. We refer to a multi-arm design that employs differential randomization as an 'inclusive design' because it maximizes participant inclusion. Considering superiority comparisons for a normal endpoint, we evaluate the performance of some analysis methods for a three-arm inclusive design by simulation studies. Methods include pairwise or overall regression analysis, an averaging approach that pools summary statistics using prevalence rate, and a meta-analysis type framework. We compare the statistical power when the inclusive design employs different treatment allocation schemes. We find that using equal allocation ratio within each subpopulation leads to a higher disjunctive power than using equal allocation across arms, when the many-to-one comparisons are analysed using the corresponding pairwise data separately as opposed to using all data in a single regression analysis. We find that differential covariate-outcome relations among patients can affect the property of the inference.
Biomedical studies often collect mixed-type longitudinal data (e.g., clinical, lifestyle) to identify predictors of a binary health outcome. A common analytical challenge is to link these irregularly collected variables to the outcome, particularly when the number of potential predictors is large. While joint models (JM) can handle this complex data structure, a critical gap exists, as they lack a formal framework for variable selection, limiting their utility for identifying relevant predictors. This article fills that methodological gap by introducing two novel, structured Bayesian variable selection strategies. Our first approach (JM1) selects predictors for the binary outcome, while our more advanced approach (JM2) performs simultaneous two-level selection for both the outcome and the longitudinal trajectories. Our framework leverages shrinkage priors to handle high-dimensional predictors and interactions, preventing overfitting. To guide final variable selection, we extend false discovery rate-based rules to our complex, multi-part joint model, which also accommodates grouped categorical predictors. We apply our method to Women's Health Initiative and Life and Longevity After Cancer study data, identifying factors for post-treatment insomnia in breast cancer survivors. Our model identified a key predictor missed by conventional methods. This integrated approach provides a robust, interpretable, and computationally efficient framework, offering substantial gains in statistical power.
A simple device for balancing for a continuous covariate in clinical trials is to stratify by whether the covariate is above or below some target value, typically the predicted median. The object is to improve balance of the covariate and hence efficiency of the treatment estimate particularly if the trial is small, as may be the case if the disease in question is rare. This raises an issue as to which model should be used for modelling the effect of treatment on the outcome variable, Y. Should one fit, the stratum indicator, S, the continuous covariate, X, both or neither? When a covariate is added to a linear model there are three consequences for inference: (a) the mean square error effect, (b) the variance inflation factor and (c) second order precision. We consider that it is valuable to consider these three factors separately, even if, ultimately, it is their joint effect that matters. We present some simple theory, concentrating in particular on the variance inflation factor, that may be used to guide trialists in their choice of model. We also consider the case where the precise form of the relationship between the outcome and the covariate is not known. We conclude by recommending that the continuous covariate should always be in the model but that, depending on circumstances, there may be some justification in fitting the stratum indicator also.
Personalized treatment selection is often based on a single primary endpoint, even though complex diseases are typically monitored using multiple, heterogeneous outcomes. We develop a rank-based framework for individualized treatment selection that accommodates mixed collections of continuous, ordinal, binary, and right-censored survival outcomes. For each outcome and treatment, we summarize the conditional outcome distribution given patient covariates using clinically meaningful parameters, including conditional quantiles, survival probabilities, and tail probabilities at ordered thresholds, and assemble these summaries across treatment arms into vectors with a common structure. Each summary component induces a ranking of the treatments, and discrepancies between outcome-specific rankings and candidate ranking vectors are quantified using a Mahalanobis-type distance that incorporates user-specified outcome weights and an empirically estimated scaling matrix. The optimal treatment is defined as the arm with the most favorable overall rank. Outcome summaries are estimated using flexible tree-based ensemble methods, although other regression or machine-learning tools can be used without altering the decision rule. Monte Carlo simulations across a range of scenarios show reasonably high probabilities of correct treatment selection, with performance improving as the signal-to-noise ratio increases and only mildly affected by the choice of outcome weights. We illustrate the method using data from an AIDS clinical trial, deriving individualized antiretroviral regimen recommendations that balance CD4 response, CD4/CD8 ratio, and a composite time-to-event endpoint.
Health technology assessment frequently requires survival predictions well beyond the observed trial follow-up, yet single parametric models fitted to short horizons can accumulate long-term bias. Current guidance also encourages the principled use of external evidence, such as registry summaries or expert anchors, while maintaining fidelity to the trial data. We introduce an adaptive spline-weighted blended extrapolation on the cumulative-hazard scale that unifies these aims. The observation side is fitted on the trial window using a piecewise-exponential Cox-type model with a smooth prior-driven continuation beyond the administrative cutoff, implemented via INLA. The external side is an anchored Gompertz tail identified by a prespecified survival level at a clinically relevant time. To blend the two components, we compare their cumulative hazards, learn a monotone P-spline score over time, and pass it through a logistic link to obtain a data-driven, time-varying weight. A simple "temperature" scaling controls the slope of this weight and provides a practical guarantee of non-negative blended hazards on a chosen grid, preserving adaptivity while ensuring feasibility. Across Monte Carlo scenarios spanning multiple tail shapes and censoring levels, the method delivers consistently lower absolute survival error, smaller restricted mean survival time error, and improved stability compared with fixed-schedule blending and single-family parametric models. In a SEER registry study with three-year observation and ten-year extrapolation, the blended curve tracks the Kaplan-Meier estimates within follow-up and transitions smoothly toward the anchored tail, yielding small long-horizon errors across cancer sites and age strata. The framework is modular, interpretable, and easily extended to alternative tails and multiple anchors. An open-source implementation is available in the R package survblendr at https://github.com/haohaostats/survblendr.
While most methods for missing not at random (MNAR) data in regression models address MNAR outcomes assuming fully observed predictors, real-world observational health and longitudinal studies often violate this assumption. This paper compares several approaches for handling MNAR data in linear regression when missingness depends on both partially observed outcomes and predictors. Through extensive simulations, we evaluate complete-case analysis, multiple imputation assuming missing at random, maximum likelihood estimation via the Heckman selection model, uncertainty intervals, multiple imputation under the Heckman selection model, not-at-random fully conditional specification, imputation stacking, and random indicator imputation. None of the methods consistently produced unbiased estimates or nominal coverage across all scenarios. However, not-at-random fully conditional specification was straightforward to implement and yielded coverage close to the nominal level in most scenarios, provided that the sensitivity parameters were specified near their true values. Our results highlight the importance of sensitivity analyses exploring various full-data models and careful parameter specification when addressing MNAR in both outcomes and predictors. We illustrate such sensitivity analyses using Betula study data on the relationship between longitudinal memory change and grey matter volume in aging. The association remained significant across most considered MNAR scenarios, reinforcing existing evidence for this relationship.
Mediation analysis is a powerful tool for exploring the causal relationships between exposures and outcomes that are mediated by intermediate variables. In this paper, we propose a flexible Bayesian mediation analysis framework to accommodate zero-inflated mediators, compatible with a wide range of outcome distributions. This novel technique employs Bayesian models for both the mediator and the outcome, utilizing Markov Chain Monte Carlo algorithms for parameter estimation. While addressing the challenges posed by an excess of zeros, we further decompose the mediation effects into components influenced by either the probability of zero or the mean of the non-zero distribution in the mediator. An associated R package mediationBayes (https://github.com/jhcuibst/mediationBayes.git) has been developed to facilitate the application of this framework. Through comprehensive simulation studies, we demonstrate that our method outperforms alternatives in terms of the accuracy of point estimates, coverage probabilities, and the precision of the mediation effects decomposition. We further illustrate the practical applicability of our method by conducting an analysis on the REasons for Geographic And Racial Differences in Stroke Study to investigate the mediating influence of smoking pack-years on the association between educational levels and incident hypertension, where mediation effects are quantified on a risk ratio scale for binary outcomes.
Compartmental models based on ordinary differential equations quantifying the interactions between susceptible, infectious and recovered individuals within a population have played an important role in infectious disease modelling. The aim of the present paper is to explain the link between stochastic epidemic models based on the susceptible-infectious-recovered (SIR) model and methods from survival analysis. We illustrate how standard software for survival analysis in the statistical language R can be used to estimate pivotal parameters in the stochastic SIR model in the very much idealized situation where the epidemic is completely observed. Extensions incorporating interventions, age structure and heterogeneity are explored and illustrated.
Joint models for longitudinal and event-time data can be used to understand the co-development of multiple correlated non-fatal outcomes, such as those encountered when studying diabetic microvascular complications. However, such models were not developed for event-times characterized by a multistate process where the entry time into certain states is interval-censored, which is an important feature of interval cohort study designs. Our aim was to develop a joint model for this setting. Specifically, we formulated a shared random effects joint model with linear mixed-effects and proportional intensities progressive three-state Markov submodels, with interval-censoring of entry into the transient state and exact observation of entry into the absorbing state. Maximum likelihood estimation was used, with exploration of the functional forms and association structures. Bias and confidence interval coverage of the model parameters were assessed using simulation, and performance was compared to existing approaches that assume exact observation of all entry times, which were found to sometimes overestimate the magnitude of association between the longitudinal and multistate outcomes. In an application of the model, the trajectory of a routinely measured diabetic complication (retinopathy) was associated with progression through states of a more difficult to measure complication (neuropathy), suggesting potential efficiencies in screening and monitoring.
Individualized treatment regimes (ITRs) represent decision-making frameworks that tailor treatment assignments to individual patient characteristics. The value function of an ITR quantifies the expected outcome under a counterfactual scenario in which such a treatment rule is applied. However, estimating optimal ITRs for survival data remains a significant challenge when outcomes are right-censored and only a subset of patients has complete outcome information due to time and cost constraints. To overcome this challenge, we formulate the problem within a semi-supervised learning framework and adopt an induced missingness perspective to model partially observed survival outcomes. We propose an imputation-based semi-supervised approach that is robust and adaptable to various imputation models. Specifically, we employ a flexible single-index kernel smoothing imputation technique to effectively utilize unlabeled data in multidimensional covariate settings. The proposed estimators for the parameters indexing the optimal ITRs are shown to be consistent and asymptotically normal. Moreover, semi-supervised estimation enhances efficiency by reducing asymptotic variance relative to supervised estimation. Numerical experiments on both simulated and real datasets demonstrate the superior performance of our proposed semi-supervised approach.
Current status data are frequently encountered in many real life cross-sectional epidemiological, demographic, and medical studies, where each subject is examined only once, and the failure time of interest is never exactly observed but known to be either smaller or larger than the examination time for each subject by evaluating the failure status. Consequently, current status data are a mixture of left-censored and right-censored observations for the failure times of all subjects with or without covariates. In some real life studies, the test or diagnosis that is used to determine the failure status may be error-prone, and this leads to misclassified failure status for some or all subjects. The resulting data are referred to as misclassified current status data in the literature. In this paper, we study regression analysis of misclassified current status data and propose a novel estimation approach under the proportional odds model. Specifically, monotone splines are adopted to approximate the baseline odds function, and an efficient expectation-maximization algorithm is developed based on a data augmentation involving exponential and Poisson latent variables. An extension of the proposed method is also developed to account for the unknown test accuracy. The proposed method is shown to have excellent estimation performance in our simulation studies and is illustrated by an application to uterine fibroid data.
We develop a robust bias-corrected method of inference about the ratio of age-standardized rates (RASR) for comparing the age-standardized rate (ASR) between a subpopulation and the whole population. Unlike previous methods, the proposed approach does not rely on the proportional age-distribution (PAD) assumption, which is often unrealistic in many situations. Like an existing approach, the method corrects for bias resulting from sampling errors when using sample-based population estimates, instead of census-based populations, as the denominators for estimating ASRs. This broadens the applicability of the proposed method in studying cancer risk factors beyond the basic demographic characteristics. The robust bias-corrected estimator of the RASR and the associated variance estimator and confidence intervals are derived. We show empirically that the proposed RASR estimator performs significantly better than the existing estimator, which relies on the PAD assumption, especially when the latter assumption fails. Specifically, the proposed RASR estimator significantly reduces the bias without increasing the variance. On the other hand, when the PAD assumption holds, our RASR estimator performs similarly to the existing estimator. The proposed method has also shown highly desirable performance when at-risk population estimates used for calculating ASRs are subject to sampling errors. We also show empirically that the proposed variance estimator performs satisfactorily. A real-data application is discussed.
There is a growing interest in subject-specific predictions using neural networks, as large-scale biomedical data often exhibit dependency due to high-cardinality categorical features, which have been largely overlooked by traditional neural network frameworks. This article proposes a novel hierarchical likelihood learning framework that captures both nonlinear overall effects and subject-specific effects by incorporating gamma random effects into Poisson neural networks. The global maximizer of the proposed objective function yields maximum likelihood estimators for fixed parameters and best unbiased predictors for random effects. The proposed framework provides a robust end-to-end algorithm for clustered biomedical count data, in the sense that the corresponding estimating equations remain unbiased even when the random-effects distribution is misspecified. To enhance learning efficiency, we introduce an adjustment procedure for the random effects and variance component. Extensive simulation studies and real data analyses demonstrate the practical effectiveness of the proposed method for clustered biomedical count data. The proposed method achieves competitive predictive performance in terms of mean squared Pearson error and mean deviance across various random-effects distributions and real-world datasets.
In medicine, multiple continuous outcomes are often repeatedly measured on each subject over time to assess disease severity. Usually, it is of interest to investigate the association between those outcomes, which may be measured at different time points, resulting in unbalanced data. The multivariate linear mixed-effects model (MLMM) is a popular framework for this analysis. It considers the unbalanced nature of the data and accounts for the association of the outcomes via the random effects, often assuming a multivariate normal distribution. However, measuring and understanding the degree of connection between longitudinal outcomes remains challenging. We propose to enhance the MLMM by incorporating various interpretable association structures. Specifically, we consider that multiple longitudinal outcomes are related to the primary outcome through their current value, cumulative effect (total or partial), or both. Our research is motivated by Pompe disease, a rare, inheritable, progressive metabolic myopathy. Clinically, it is important to investigate how patient-reported outcome measures (primary outcomes) are associated with physical outcomes to determine whether improvements in physical outcomes are accompanied by improvements in health-related quality of life and other patient experiences. We found a positive association between them. The proposed models are fitted under the Bayesian framework using Hamiltonian Monte Carlo.
Recurrent health events often involve complex inter-relationships between longitudinal biomarkers and time-to-event outcomes, further complicated by sparse, irregular data collection and time-dependent correlations among events. Traditional statistical methods frequently struggle with these complexities, resulting in biased estimates and suboptimal modeling performance. To address these challenges, we propose the F unctional R egression with A utoregress I ve frai LTY (FRAILTY) method, a novel framework designed to jointly model longitudinal measurements and recurrent events, accommodating both scalar and functional covariates while capturing time-dependent correlations among events. The FRAILTY method employs a two-step estimation procedure. First, functional principal component analysis through conditional expectation (PACE) is applied to extract key temporal features from sparse and irregular longitudinal data. Second, the obtained scores are incorporated into a dynamic recurrent frailty model with an autoregressive structure to account for within-subject correlations across recurrent events. Simulation studies demonstrated that the FRAILTY method outperformed existing methods, such as those relying on B-spline basis functions and Bayesian joint modeling, by achieving lower integrated mean squared errors, higher concordance indices, and greater statistical power in detecting functional parameters. Its practical utility was further validated through applications to two datasets: the Systolic Blood Pressure Intervention Trial study and the Multicenter Collaboration to Study Treatment Outcomes in Nephrolithiasis Evaluation cohort.
In medical diagnostic studies, the area under the receiver operating characteristic curve (AUC) is a widely used metric that captures a continuous test's overall ability to discriminate between diseased and non-diseased individuals across all possible cutoffs. However, in practice, disease status is sometimes only partially verified, introducing verification bias that undermines the validity of AUC estimation. While numerous methods address bias correction for AUC estimation, approaches that directly construct confidence intervals for the AUC remain limited. This paper proposes two robust methods for constructing bias-corrected confidence intervals for the AUC under the missing-at-random assumption: one based on bootstrap resampling and the other on empirical likelihood. Both approaches accommodate missing disease verification by leveraging the bias-corrected ROC estimators introduced by Alonzo and Pepe. Extensive simulation studies and real-world data analyses demonstrate that our proposed methods yield valid and precise interval estimates for the AUC under various clinically relevant settings.
Maximum likelihood estimates of density function and regression coefficients in the proportional odds regression models are proposed and studied based on event-time data that are either completely or partly interval-censored. A smooth estimate of the survival function is then obtained. Theoretical results indicate that the proposed method enjoys an almost parametric n-consistency. Some simulation studies show that the proposed method not only gives density and smooth survival curve estimates but also outperforms the existing semiparametric method in terms of estimating both the regression coefficients and the survival curves for small and medium sample sizes. The proposed method is illustrated by an application to the HIV infection data.
