In survival data analysis, the stratified Cox model becomes a popular option when the proportional hazards assumption of the conventional Cox model does not hold for certain covariates. For a stratified Cox model, when the observed survival times contain only right censoring, the method of maximum partial likelihood can still be implemented. However, if survival times include interval-censored observations, the method of maximum partial likelihood is not viable, and the partial likelihood approach cannot be applied. Furthermore, partial likelihood analysis does not supply a smooth estimate of the baseline hazard. In this paper, we consider the stratified Cox model under partly interval-censored survival times. We present a penalized likelihood method for estimating the model parameters, including the baseline hazards. Penalty functions are used to produce smoothed baseline hazards estimates, and also to relax the requirement on optimal number and location of the knots used in the baseline hazards estimates. We also derive a large sample normality result for the estimates, which can be used to make inferences on quantities of interest, such as survival probabilities, without relying on computing-intensive resampling methods.
Bounded lifetimes on (0,1), such as mortality proportions, are prevalent in biomedical and public health research; however, mean-based models can overlook tail behavior and heterogeneity. We introduce the unit power Maxwell distribution (UPMD) defined on (0,1) and develop a quantile regression framework that captures covariate effects across the outcomes. The UPMD is generated based on the log power transformation of the Maxwell-Boltzmann distribution, producing a tractable quantile and flexible shapes with skewness and heavy tails. This work examines several key distributional properties, including the quantile function, order statistics, residual life, and Rényi entropy. It also investigates inference via maximum likelihood, maximum product spacing, least squares, and Bayesian estimation. Extensive simulation analysis is performed to evaluate the efficiency of various estimation methods using bias, mean squared error, and interval probability coverage (CPs). For empirical evaluation, we analyze COVID-19 mortality proportions, where UPMD-based quantile regression delivers improved fit in upper and lower tails and better calibration relative to alternatives, as assessed by information criteria, formal goodness of fit tests, and residual diagnostics such as randomized quantile residuals and Cox-Snell checks. The results indicate that the UPMD provides a practical, interpretable tool for bounded biomedical outcomes.
While randomized controlled trials remain the gold standard for estimating causal effects in medical research, they are not always practical. Consequently, observational data become a necessary alternative, though it introduces confounding bias due to the lack of randomization. Propensity score modeling mitigates this bias but is highly sensitive to covariate selection. To address this critical limitation in survival analysis, we propose the outcome-adaptive Lasso (OAL)-Cox model-a novel variable selection framework for causal inference for right-censored data. The proposed method integrates the OAL with the Cox proportional hazards model by constructing penalty weights from coefficients estimated through the Cox partial likelihood. This design enables OAL-Cox to effectively identify outcome-related confounders while excluding irrelevant covariates, including instrumental and spurious variables, even under censoring and correlated covariate structures. Simulation studies demonstrate that OAL-Cox achieves stable variable selection across different censoring rates and correlation settings. Compared with the L-Cox and AL-Cox methods considered in this study, OAL-Cox more effectively reduces the inclusion of instrumental and spurious covariates while retaining outcome-related covariates, and it generally yields more stable restricted mean survival time-based causal effect estimates in the examined settings. We further illustrate the practical applicability of the proposed method by applying it to ovarian cancer recurrence time data to evaluate the effect of surgical sophistication on recurrence time.
In estimating the average treatment effect (ATE), the plug-in estimator with the efficient influence function and the targeted maximum likelihood estimator (TMLE) are semiparametrically efficient. However, the estimators suffer from the fact that too small/large a propensity score (PS) in the denominators can make the estimators unstable. This is usually overcome by trimming or truncating the data, that is, discarding observations with the PS too small/large or replacing the PS with some lower or upper threshold, which amounts to using a non-smooth weight to move the target parameter away from the ATE. This paper proposes a smooth weight, quadratic in PS instead, which leads to targeting the "overlap-weight (OW) treatment effect" that has several advantages over the ATE. We compare various estimators through simulation and empirical studies, including a trimmed efficient influence-function estimator, TMLE, and their OW versions. We find that the OW versions perform better than the trimmed ones in terms of the absolute bias and root mean squared error, when the ATE is a constant so that all estimators share the same estimand. We also find that "cross-fit" analogous to cross-validation improves the asymptotic variance estimators.
Diseases often involve multiple dimensions of interrelated impairments. Although significant advances have been made in joint models to simultaneously assess these processes in relation to clinical endpoints, they often fail to evaluate how these dimensions influence each other. We propose an original joint modeling framework to describe the temporal relationships between the processes involved in Alzheimer's disease and related dementias (ADRD) progression, and assess their association with ADRD diagnosis and death. The longitudinal submodel is a dynamic model that combines a structural multivariate mixed model based on differential equations to explain the instantaneous change over time of each latent dimension according to the others, and observation models that can accommodate ordinal, binary, and continuous (whether Gaussian or non-Gaussian) biomarkers. The association of the biomarkers with ADRD diagnosis and death are described via a shared random-effect joint modeling approach. The estimation procedure, carried out within the maximum likelihood framework is made available in the DynNet R package. The methodology is validated in a simulation study and is applied in a population-based French cohort study to disentangle the temporal relationships between three major drivers of ADRD natural history, depression, cognition, and functional dependency in link with the two major clinical events in ADRD progression: ADRD diagnosis and death. The methodology and application are designed to help understand the complex interplay between biomarkers over time.
A large proportion of clinical trials do not meet their recruitment targets. Trials using time-to-event endpoints come along with the additional complexity that the amount of available information depends on the number of observed events, which is random for a specified point in time. This means that even if a trial may have recruited the planned number of patients in its recruitment interval, the amount of available information at a fixed calendar time is random. The determination of the sample size underlying a clinical trial is a very important aspect when designing a trial. When adaptations of the pre-specified sample size are a desired design option to address the uncertainty of underlying assumptions, blinded interim analyses are recommended by the authorities. However, blinding cannot be assured when there is a high interest in a potential early trial stop or a high uncertainty about the underlying effect size. In those situations, (adaptive) group-sequential trial designs may be used. In this work, we propose a flexible, hybrid approach. Based on the chances of an early trial stop, either a blinded or an unblinded interim analysis is conducted such that the advantages of the two commonly applied approaches are combined. We focus on trials where the time to an event is the primary outcome. These studies usually come along with a long study duration, where the necessity of design adaptations may be especially appealing.
Bayesian conditional transformation models (BCTMs) address the direct estimation of the conditional distribution function of a random variable Y $Y$ conditional on a set of explanatory variables X ${\rm variables}\ \bm{X}$ . The BCTMs infer the conditional distribution by applying a transformation function of Y $Y$ given X = x $\bm{X} = \bm{x}$ towards a baseline distribution free of parameters to be estimated. The benefit of these models is that the explanatory variables X = x $\bm{X} = \bm{x}$ impact the whole conditional distribution of Y $Y$ given X = x $\bm{X} = \bm{x}$ instead of only the mean, variance, kurtosis, or skewness. The transformation functions are an essential part of the model, and they range from loss-complex and low-parameterized functions to complex relationships between explanatory variables and response variables represented by nonlinear functions. The general construction of the BCTM class explores monotonic B-splines for parameterizing the transformation function. Smoothness and regularization are accomplished through an adequate prior distribution for the parameters. We proposed a new estimation procedure for the BCTM based on the integrated nested Laplace approximation, which is tested through a simulation study. Also, two longitudinal studies using real data are considered. The first application is a cardiovascular study and compares our proposed algorithm, named integrated Laplace with Bayesian conditional transformation models (ILBCTM), with the original Markov chain Monte Carlo-based algorithm for BCTM. We obtained similar results with a shorter computational time. The second application considers the ILBCTM in a study of the mortality rate of bronchial and lung cancer in Brazil.
Beyersmann et al. propose a functional interpretation of hazards, viewing them as evolving quantities describing the entire event process rather than as pointwise causal contrasts. In this commentary, we elaborate on the implications of this view for causal inference in modern clinical trials with survival outcomes. We emphasize how censoring, competing events, and multistate structures shape not only identifiability but also the definition and transportability of hazard-based estimands. We highlight that, even within a functional framework, censoring mechanisms may implicitly determine the statistical estimand through time-dependent weighting, with direct implications for generalizability across studies and populations. We further discuss how these issues are amplified in competing-risks and multistate settings, where causal interpretation requires careful consideration of intercurrent events and selection induced by post-randomization state occupancy.
Meta-analysis of diagnostic test accuracy studies aggregates information from multiple studies on sensitivity and specificity. Classical approaches select a single pair of sensitivity and specificity per study (single threshold methods, STM), ignoring additional information if studies report results on multiple diagnostic thresholds. Recently, models have been proposed that consider all available information and enable inference on all diagnostic thresholds (multiple threshold methods, MTM). We compare five STM and six MTM to each other in a simulation study, evaluating their performance in various situations. Covering a broad range of real-life settings, we vary eight parameter dimensions in the data-generation mechanisms, including continuous or ordinal outcome type of an index test, and different numbers of diagnostic thresholds available per study. While model performances are comparable regarding bias, empirical coverage, and convergence, we observe a logit GLMM of the MTM type to perform best in many situations. Model performances depend strongest on the outcome type, while the number of thresholds only has a minor impact. We thus find the main advantage of using MTM by getting threshold-dependent estimates of sensitivity and specificity. Additionally, we illustrate differences between model estimates in two real-data examples on diagnosing type 2 diabetes using the continuous biomarker HbA1c and screening for any anxiety disorder using the ordinal questionnaire HADS-A. The applications reveal variations in model estimates within and between STM and MTM, which can be reduced by adjusting for the bias in the simulation settings resembling the real-data situation most closely.
Most original articles published in the medical literature report the results of multiple statistical tests. In a few simple cases, there is agreement on whether to adjust for the number of performed tests. For many cases encountered in practice, however, this is less clear, and the recommendations in the literature are contradictory, along different dimensions, or otherwise confusing. This lack of clear guidance may impair the conduct and interpretation of analyses, and encourage questionable research practices, ultimately jeopardizing the credibility of medical research. In this article, we refine, illustrate, and discuss a unifying guiding principle to assist both statisticians and applied researchers in deciding whether to adjust for multiple testing and, if so, over which set of tests. The principle is that multiple testing should be adjusted for if and only if authors, when reporting and interpreting their findings, put more emphasis on results of one or several of the tests because of their small p-value(s). We relate this principle to previously proposed rules and show how it can guide and clarify the choice of adjustment strategies in three complex multiple testing settings.
A hierarchical 2 × 2 $2\times 2$ factorial design is a type of two-level trial design where the first intervention is randomized at the cluster level and the second intervention is randomized at the individual level. With a continuous outcome, the linear mixed model with a random intercept can be used for this design to estimate the treatment effects while accommodating the within-cluster correlations. Such a model often serves as the basis for the sample size and power calculation. However, recent evidence in cluster randomized trials has shown that the cluster-level intervention effect might differ across clusters, leading to extra variability in the outcomes. This paper extends the existing literature on designing hierarchical 2 × 2 $2\times 2$ factorial trials to address heterogeneous treatment effects across clusters. Under the generalized least squares framework, we consider models with or without an interaction between the two interventions, and derive sample size formulas for testing the controlled effects, the marginal effects, and the interaction effect between the two treatments. Simulation studies were conducted to verify our sample size formulas in finite samples. The context of hierarchical 2 × 2 $2\times 2$ factorial trial on suicide prevention is used for illustrating our methods.
This article describes the design of a neutral comparison study in the context of empirical studies where the interest is in learning the functional relationship between a continuous error-prone exposure variable and a binary outcome. The performance of combinations of measurement error correction methods and flexible regression modeling techniques was compared using a simulation study. The project involved four independent teams, one devoted to Data Generation and Evaluation and the other three to specific Methods of measurement error correction (regression-calibration and multiple imputation, simulation-extrapolation, and Bayesian method). The study was conducted in three successive stages. In Stage 1, the first team simulated five datasets differing only by the true exposure-outcome functional form and distribution of true exposure. Furthermore, the implementation of flexible modeling methods (B-splines, P-splines, and fractional polynomials) was standardized. The three Methods teams, blinded to the underlying data generation process, created the codes to implement their methods and provided their results to the first team, who evaluated them. These codes were then used by this team in the next stages of the project. In Stage 2, the team simulated 150 additional datasets where other design parameters varied while using the same five exposure-outcome functions. Stage 3 consisted of simulating independent replications of each of the 150 scenarios considered in Stage 2 to quantify the sampling variance of the estimates. This work emphasizes the relevance of neutral comparison studies to fairly evaluate statistical methods aimed at addressing a complex analytical challenge and demonstrates their feasibility through a large collaborative project.
Hierarchical Composite Endpoints (HCEs), as analyzed with Generalized Pairwise Comparisons (GPC) statistics, are general methods of constructing endpoints in clinical trials across various therapeutic areas to establish the efficacy of novel treatments. Although GPC statistics do not require distributional assumptions for estimation, ignoring the underlying distributions can hinder the interpretation of treatment effects. We provide a formal definition of HCEs in a special case based on the "most-important outcome" principle over a fixed timeframe of evaluation. This can be a limitation, but it offers important advantages. In this case, although HCEs involve multiple outcomes, they result in univariate distributions, allowing for the application of the Brunner-Konietschke formula for the estimation of Mann-Whitney effect (MWE) variance. The Condorcet paradox in clinical trials can cause treatment effects to be nonaccumulative. We discuss classes of the stochastically ordered component distributions of the univariate HCE as a sufficient condition to avoid this paradox, emphasizing that HCEs should preferably be defined with outcomes where a consistent treatment benefit is expected. Maraca plots can be used to diagnose violations of conditions that could lead to a Condorcet paradox. We also provide a rank-based method for estimating the MWE when a threshold is used for pairwise comparison between patients in two groups and discuss the implications of using thresholds on the statistical power for the success odds test. Our formal definition of HCEs under stochastically ordered distributions provides a strong theoretical foundation for the use of HCEs in clinical trials, including when a threshold is applied to pairwise comparisons.
Human papillomavirus (HPV) is a well-established prognostic factor in head and neck (HN) cancer, with HPV-positive patients exhibiting markedly better survival outcomes compared to their HPV-negative counterparts. While advances in (cancer) genomics have been pivotal to precision medicine, existing gene screening methods for identifying molecular markers to predict survival often fail to account for HPV status. This oversight can result in missing important genes, whose effects are confounded or overshadowed by HPV, thereby limiting the biological interpretability and clinical utility of identified markers. To address these limitations, we propose a novel conditional screening method for ultrahigh-dimensional right-censored survival data that adjusts for HPV status. This approach identifies prognostic genes with independent associations with survival while also capturing HPV-specific interactions and synergistic effects. The proposed method employs a two-stage, model-free framework that combines nonparametric statistics for initial screening with a unified false discovery rate (FDR) control procedure to refine feature selection. Simulation studies demonstrate its advantages over existing alternatives. Application of the conditional screening framework to HN cancer data from The Cancer Genome Atlas revealed a set of robust prognostic genes, uncovering new insights into the molecular pathways driving survival outcomes across HPV subgroups.
In this paper, we present a Monte Carlo method for estimating a nonlinear function of the mean of a multivariate normal distribution. Building on this method, we propose a parametric estimation procedure for unimodal regression models, assuming that the response variable follows a gamma distribution while some covariates are contaminated with normal measurement error. Compared to existing approaches, the proposed method accommodates multivariate covariates and features a tractable bias-corrected likelihood function, enabling faster computation and more accurate estimation when the data distribution is correctly specified. To enhance the applicability of the proposed method, we also explore various model adequacy diagnostic tools and evaluate its robustness against distributional misspecifications. Notably, we introduce a goodness-of-fit test based on a unique characterization of the gamma distribution, designed to assess the validity of the distributional assumption for the response variable. Numerical studies and real-world data applications are conducted to evaluate the finite-sample performance of the proposed methods.
Although sample size calculation for open-cohort longitudinal cluster randomized trials (LCRTs) under a fixed design framework was developed by Kasza et al., unifying the closed-cohort and repeated cross-sectional sampling provided in Hooper et al. when a churn rate is constant, there has been no prior efforts in developing optimal open-cohort LCRTs that maximizes the design efficiency. This work assumes a prespecified number of periods T $T$ and a constant number of replaced individuals at each period in open-cohort LCRTs. We propose algorithms for deriving optimal sample size under a cost-efficiency framework and arrive at the local optimal design (LOD) for fixed correlation parameters and MaxiMin optimal design for addressing uncertainty in correlation parameters. When correlation parameters are known, as the number of replaced individuals increases, for open-cohort PA-LCRTs, the optimal cluster-period size generally decreases and then increases whereas the optimal number of clusters and power under LOD first increase and then decrease. In contrast, for CRXO trials and standard SW-CRTs, the optimal cluster-period size and churn rate under LOD increase whereas the optimal number of clusters and power under LOD decrease. When correlation parameters are unknown, but the parameter space is available, with a small number of replaced individuals, there is no difference in optimal designs between PA-LCRTs and CRXO trials. The number of replaced individuals also has less impact on the optimal cluster-period size than optimal number of clusters. We demonstrate our new optimal design methods using the context of two real-world LCRTs.
Unmeasured confounding remains a fundamental challenge permeating contemporary causal inference, substantially impeding the valid estimation of treatment effects. A rigorous identification strategy is presented for estimating average causal effects with unmeasured confounding by exploiting the information contained in primary and secondary outcomes. In contrast to existing literature, our approach is intended to construct the proxy confounder for inverse probability weighting-type estimation. Formal identification results and the asymptotic distribution theory for the proposed estimator are established. Through extensive simulation studies, it is demonstrated that the method achieves marked reduction in confounding bias and offers refinements to causal effect estimation. In practical applications, by integrating secondary outcomes that characterize cognitive aspects, we successfully supplemented information not captured by the covariates, enabling us to draw significant inferences regarding the effects of maternal delivery mode on child's test scores. This approach provides a promising methodology for data sets with multiple secondary outcomes.
Seamless phase II/III design aims to integrate a phase II trial for treatment selection and a phase III confirmatory trial. It offers valuable flexibility through mid-trial modifications, potentially optimizing resource utilization and reducing patient burden. In the first stage, multiple experimental treatments are evaluated to select promising ones, which enter the second stage for the confirmatory analysis against the control using data from both stages. Proper statistical methods are needed for the final analysis to control the overall Type I error rate at a prespecified level, regardless of the treatment selection rule used at the interim analysis. In this article, we provide a comprehensive evaluation of four classes of statistical approaches in the literature, which use different ways to integrate combination functions (or conditional error functions) and multiple testing methods (including Bonferroni, Simes, and Dunnett adjustments). Extensive simulation studies are performed to evaluate both the Type I error control and power. In addition, we illustrate the practical implementation of these approaches in real clinical trial settings.
Single gene mutations are increasingly being adopted as clinical biomarkers for the optimal application of various therapeutic areas (such as cancer and cardiovascular disease). A single nucleotide polymorphism (SNP), the most common type of genetic variation in human populations, can affect the abundance and function of gene products at the molecular level. In practice, a single-predictor model is usually built using the best single-SNP (or single-gene) biomarker with the most significant drug-SNP (or drug-gene) association in randomized clinical trials. Current statistical methods for the best single-SNP biomarker selection rely on a variety of ranking procedures. However, these existing approaches (i) cannot accurately distinguish drug-SNP interaction effects (i.e., predictive effects) from SNP main effects (i.e., prognostic effects); (ii) do not necessarily yield the true ordering or provide a confidence interval for ranking. In this paper, we propose a novel method called Multiple comparisons with the best using Marginal Means (3M). Specifically, 3M first calculates unbiased estimates of SNPs' predictive effects using marginal means (i.e., least squares means). To select the best predictive SNP (i.e., SNP with the largest unbiased drug-SNP interaction effect), 3M uses a nonparametric bootstrap method to construct constrained simultaneous confidence intervals under the framework of Multiple Comparisons with the Best. Simulation studies demonstrate that our proposed method is more reliable in distinguishing predictive effects from prognostic effects, and thus more powerful to identify the true best predictive SNP than existing methods. Finally, we applied our method to the IMPROVE-IT (IMProved Reduction of Outcomes: Vytroin Efficacy International Trial) pharmacogenomics genome-wide association study (GWAS) data to search for the best predictive SNP across the genome. One single best predictive SNP rs114462013 (STAG1) on chromosome 3 was detected, with the association previously identified in the literature.
