Statistical Methods in Medical Research

Randomized clinical trials in critically ill populations often use mortality at the end of follow-up as the primary outcome. When mortality is rare, composite outcomes that capture disease severity over time are frequently employed to improve efficiency. Motivated by the Crystalloid Liberal or Vasopressor Early Resuscitation in Sepsis (CLOVERS) trial—a randomized study comparing two fluid resuscitation strategies in hospitalized patients with sepsis—we evaluate multiple modeling approaches for estimating covariate-adjusted treatment effects on absorbing outcomes such as mortality. While analyses of binary end-of-follow-up outcomes (alive vs. dead) are natural, efficiency can be gained by incorporating low-cost longitudinal and ordinal severity data routinely collected in trials and by imposing reasonable modeling restrictions (constraints). We develop methods that leverage these data within the marginalized model framework, focusing on first-order marginalized transition models, which allows direct estimation of outcome-marginal (population-averaged) treatment effects. We highlight a key property of the marginalized transition model in the presence of absorbing states: the marginal and conditional sub-models correspond to cumulative incidence and discrete hazard models, respectively. Through simulation studies and an application to the CLOVERS data, we compare the validity and efficiency of competing approaches, including standard survival analyses.

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.

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.

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.

The development of chronic disease is a long-term process involving multiple endpoints. Although multi-state Cox models can estimate state-specific survival risks over time, they are not well suited for comparing the effectiveness of treatment regimes. A discrete-time split-state framework has been proposed, which divides disease states into substates by conditioning on past history. As this framework is both “memoryless” and “memorable,” the transition rates can be synthesized into summary measures, multimorbidity-adjusted life year (MALY) and substate-specific life year (SSLY). Building on this framework, we propose to investigate the causal effects of static and dynamic treatment regimes over the disease course under the assumptions of constant baseline confounders and instantaneous effects of interventions on transition rates. We identify the optimal treatment regime as the one that maximizes MALY and use SSLY to elucidate the mechanisms of how treatments influence disease progression. In the application, we identified the optimal weight targets in the ARIC study by modeling the disease course in healthy, at-metabolic-risk, coronary heart disease, heart failure, and mortality states. The estimated MALY was 1.80 years higher (95% CI: 0.62, 2.78) under regime “Normal weight initially and change to overweight if age >65 y” compared to regime “Normal weight across all states.” SSLY decomposition indicates that this gain arises from increased life year in all substates except the healthy state. In summary, our method provides a framework to evaluate the health benefits of treatment regimes over the disease course and has the potential to improve the precision prevention of chronic diseases.

Despite extensive research, the use of response-adaptive randomization (RAR) in clinical trials has remained controversial. Korn and Freidlin’s 2011 article reignited this debate back then, prompting numerous responses, including one by Zhu, Rosenberger, and Hu that remained unpublished until now. This article features the original response by Zhu, Rosenberger, and Hu, providing a valuable opportunity to revisit the original arguments, examine subsequent developments in RAR methodology, and offer a more complete historical perspective on this enduring debate. The piece also includes a concluding section by one of the special issue’s co-editor that explores the nuances and complexities of RAR implementation in the context of contemporary clinical trial design.

The Cox proportional hazards regression model is a widely used and valuable tool for modeling survival time with predictors, however its performance can deteriorate in the presence of multicollinearity. It can lead to unreliable estimates from the maximum partial likelihood estimator. In this paper, we introduce enhanced shrinkage estimators based on the Kibria–Lukman approach to obtain more efficient coefficient estimates. Specifically, we develop linear shrinkage, Stein and its positive counterpart, pretest, and shrinkage pretest estimators that incorporate prior information about model coefficients. We derive their asymptotic bias and variance properties and evaluate their performance through extensive Monte Carlo simulations. The results demonstrate significant improvements, highlighting the practical advantages of these methods for applied researchers. We also illustrate the application of our proposed estimators using a lung cancer dataset.

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.

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.

Response-adaptive randomization (RAR) can increase participant benefit in clinical trials, but also complicates statistical analysis. The burn-in period—a non-adaptive initial stage—is commonly used to mitigate this disadvantage, yet guidance on its optimal duration is scarce. To address this critical gap, this paper introduces an exact evaluation approach to investigate how the burn-in length impacts statistical operating characteristics of two-arm binary Bayesian RAR (BRAR) designs. We show that (1) commonly used calibration and asymptotic tests show substantial type I error rate inflation for BRAR designs without a burn-in period, and increasing the total burn-in length to more than half the trial size reduces but does not fully mitigate type I error rate inflation, necessitating exact tests; (2) exact tests conditioning on total successes show the highest average and minimum power up to large burn-in lengths; (3) the burn-in length substantially influences power and participant benefit, which are often not maximized at the maximum or minimum possible burn-in length; (4) the test statistic influences the type I error rate and power; (5) estimation bias decreases quicker in the burn-in length for larger treatment effects and increases for larger trial sizes under the same burn-in length. Our approach is illustrated by re-designing the ARREST trial.

Response-adaptive randomization (RAR) dynamically adjusts allocation probabilities of sequentially enrolled patients based on accumulating response information. It has gained increasing attention in clinical trials for its ability to enhance statistical efficiency by accelerating estimation of treatment effects, improve ethical allocation by favoring superior interventions, and maintain robust parameter estimation (e.g. via variance or uncertainty reduction) while accounting for adaptive modifications. While both regulatory agencies and the broader clinical research community are making efforts to promote the application of RAR, its real-world implementation remains limited. In this article, we review the application of RAR in clinical practice from 2015 to 2024, identifying key challenges such as managing diverse patient outcome types, incorporating repeated measurements, and addressing missing data. We categorize various RAR methods and propose practical solutions to these challenges, providing insights and practical references for clinical trial practitioners. Additionally, we discuss the limitations of existing RAR methods and outline potential future research directions. Our review aims to bridge the gap between theory and practice, promoting the broader adoption of RAR in clinical trials while advancing the development of RAR theory and methodology.

Joint modelling of longitudinal and time-to-event data is a powerful tool for analysing complex medical data. Joint latent class models (JLCMs), in particular, provide additional prognostic insight by clustering subjects into latent subgroups within a mixture of experts framework. Traditionally, the mixing proportions are modelled using a multinomial logistic regression. Recent work in mixture regression has proposed a neural network alternative, implemented through the neural network expectation–maximization algorithm, which offers greater flexibility by capturing nonlinearities that may be missed by logistic regression. In this article, we extend the JLCM by applying the neural network to the modelling the mixing proportions captured in the class-membership submodel. The model estimated using an expectation–maximization-type algorithm. We call this the neural-network JLCM (NN-JLCM). We formally define the method, evaluate its performance through extensive simulations, and apply it to real data. The results show that NN-JLCM improves predictive accuracy and class separation compared with traditional JLCMs, while maintaining competitiveness in terms of overall model fit. The code for the article is available at https://github.com/tharris0924/nn-jlcm_p1.

Monitoring medical imaging data is increasingly important in healthcare because image-derived quantitative features provide valuable information for assessing disease progression, diagnostic quality, and process stability. However, such data are often feature-rich, distributionally complex, and limited in sample size, and existing methods do not always provide a unified way for retrospectively monitoring joint changes in mean and covariance without relying on strong parametric assumptions. Motivated by these challenges, this paper develops a Phase I nonparametric change-point monitoring framework for individual observations. The framework combines existing high-dimensional two-sample statistics for mean and covariance changes into two control chart schemes, whose control limits are determined through simulation. A post-signal diagnostic procedure is also presented to estimate the change-point location, identify whether the detected shift occurs in the mean, the covariance, or both, and screen potentially affected variables. An extensive simulation study considering various distributional models shows that the proposed charts preserve the nominal in-control false alarm probability, achieve favorable out-of-control detection power, and yield informative post-signal diagnostic results. Finally, a real-world colposcopy dataset is used to illustrate the practical applicability of the proposed framework for the retrospective monitoring of medical imaging data.

The localization receiver operating characteristic (LROC) and the free-response receiver operating characteristic (FROC) curves are popular methods for evaluating the performance of diagnostic tests concerning detecting and locating lesions. However, the existing LROC and FROC curves have some limitations in their definition and estimation methods. In this article, we extend the LROC curve to the free-response tasks where each case can be marked with multiple suspicious locations and propose a nonparametric estimation of the LROC curve to compensate for the limitations of existing LROC and FROC curves. We derive the asymptotic properties for our estimators and propose a bootstrap method for constructing confidence bands for the proposed LROC curve. Simulations are conducted to compare the proposed LROC curve with the existing AFROC curve and to examine the theoretical results. Finally, we apply our method to evaluate the performance of a diagnostic test investigating AI-assisted pulmonary nodule diagnosis.

Count data are common in medical research. When these data have more zeros than expected by the most used count distributions, it is common to employ a zero-inflated regression model. However, the interpretability of these models is much lower than the most used count regression models. We present a more interpretable regression model that estimates the mean event rate and models covariate-dependent dispersion directly. Additionally, the dispersion parameter can be interpreted as an index of clumping, a control parameter for overdispersion. We discuss inferential and diagnostic tools and perform a Monte Carlo simulation study to evaluate the performance of the maximum likelihood estimator. Finally, the usefulness of the proposed regression model is illustrated through an application on antenatal care visits.

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.

In this paper we propose and examine randomization and allocation procedures in master protocol trials where each subtrial is performed as a single-arm study, which is motivated mainly by rare diseases. The subtrial analysis is done by comparing the single-arm observations data to historical data or expert judgment and following prespecified decision rules, without using data from any other study in the master protocol trial. The objective of this paper is to examine relevant options for participant allocation procedures in a variety of patterns in which the subtrials enter the master protocol trial, such as umbrella, platform and perpetual patterns. Our simulation study shows that our proposed allocation procedures using the predictive probability of current study success may bring substantial efficiency gains in the master protocol trial, making a clever allocation of the in-trial participants to achieve study success declarations sooner, benefiting the out-of-trial population.

The opportunity to utilize multivariate functional data types for conducting classification tasks is emerging with the growing availability of imaging data. Inspired by the extensive data provided by the Alzheimer’s Disease Neuroimaging Initiative, we introduce a novel classifier tailored for multivariate functional data. Each observation in this framework may be associated with numerous functional processes, varying in dimensions, such as curves and images. Each predictor is a random element in an infinite-dimensional function space, and the number of functional predictors can potentially be much greater than the sample size. By adopting a sparse deep rectified linear unit network architecture and incorporating the LassoNet algorithm, the proposed functional Bayesian information criterion deep neural network performs feature selection and classification simultaneously, in contrast to existing functional data classifiers. This approach addresses the challenge of complex inter-correlation structures among multiple functional processes without requiring distributional assumptions. A simulation study and a real data application demonstrate its favorable performance.

In online multiple testing, an a priori unknown number of hypotheses are tested sequentially, that is, at each time point, a test decision for the current hypothesis has to be made using only the data available so far. Although many powerful test procedures have been developed for online error control in recent years, most of them are designed solely for independent or at most locally dependent test statistics. In this work, we provide a new framework for deriving online multiple test procedures that ensure asymptotical (with respect to the sample size) control of the familywise error rate, regardless of the dependence structure between test statistics. In this context, we give a few concrete examples of such test procedures and discuss their properties. Furthermore, we conduct a simulation study in which the type I error control of these test procedures is also confirmed for a finite sample size, and a gain in power is indicated.

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.