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.
In resource-limited or time-sensitive care settings, there is interest in assessing the impact of time to treatment (TTT) on mortality. Traditional Cox proportional hazards models, which specify the effect of TTT as an unrestricted term in the log hazard ratio, can produce counterintuitive results where the survival probability may not decrease monotonically with longer delays. Moreover, hazard ratios from such models quantify the effect of TTT conditional on surviving until treatment, rather than the effect of delayed treatment at baseline. We propose a class of bounded hazard ratio (BHR) Cox models that constrain the hazard ratio for TTT to attenuate towards the null with increasing treatment delay, such that hazard for death after treatment cannot exceed the hazard without treatment. Estimation can be performed using direct optimization of the partial log-likelihood or with an iterative linearized estimation procedure for large sample sizes. From BHR models, the estimated hazard ratio curve describes how treatment benefit diminishes with delay. Additionally, we propose a survival probability difference that provides an absolute measure for comparing survival under different treatment delays. We evaluate model performance in simulations and apply the method to examine treatment delays for colon cancer with data from the National Cancer Database.
Longitudinal zero-inflated count data frequently arise in various fields such as medicine and social sciences. Standard hurdle models separate zero and positive counts, but fail to directly infer marginal means, limiting their ability to assess overall covariate effects. This paper introduces overall marginalized hurdle random effects models (OMHREMs) for zero-inflated count data, which extends the traditional hurdle model by directly modeling the marginal mean while considering random effects to account for heterogeneity. OMHREMs enable population-average effects of covariates like odds ratio, providing how covariates influence the overall mean in zero-inflated count data. Through simulation studies, we evaluate the performance of OMHREMs. Furthermore, we apply our approach to systemic lupus erythematosus data to compare its effectiveness against existing models.
Interactions and correlations among features are essential in biology, as well as in other fields. This article introduces a novel approach for linear interaction models characterized by complex correlation structures. By integrating local linear approximation and Laplacian smoothing penalty with l1 or l1 and l2 penalties, our methods effectively estimate and predict highly correlated interaction models. Theoretical analysis confirms that both methods converge to an oracle solution within two iterations, demonstrating a rapid convergence rate. In simulation studies, our proposed methods outperform existing techniques in terms of prediction accuracy, estimation precision, and variable selection. When applied to protein microarray data for Alzheimer's disease analysis, they reveal substantial main and interaction effects with notably lower prediction errors. This highlights the potential of our methods as powerful tools for analyzing linear interaction models with intricate correlations, applicable across a wide range of biological research and other fields.
