Biostatistics

Estimating the correlation between subject-specific rates of change is a common objective in longitudinal studies. A natural approach is to fit ordinary least squares regressions per subject and correlate the resulting slope estimates directly. Under complete data, this slopes-first estimator is algebraically equivalent to a contrast-based plug-in estimator derived from the joint variance-covariance matrix of the two signals across time. Under independent dropout, however, the slopes-first approach implicitly applies heterogeneous weights across missing-data patterns, introducing bias whenever the cross-covariance between the two signals is not constant over time. We formalise this result and propose an unbiased alternative that estimates the joint covariance matrix before applying the contrast. We further extend our results to include covariates as well. Simulations confirm that the slopes-first estimator is biased under non-constant cross-covariance while the proposed alternative remains unbiased. Further, we show in a real-world example that the two estimators can lead to different scientific conclusions and recommend that researchers working with longitudinal data inspect the missing data pattern to make an informed decision about the estimator of choice.

Longitudinal non-Gaussian data are common in biomedical and microbiome studies, where responses may not be normally distributed. However, many existing dynamic graphical models are developed for Gaussian responses, while many non-Gaussian graphical models focus on cross-sectional dependence and do not explicitly estimate temporal transition effects. We propose a latent Gaussian dynamic graphical model for longitudinal non-Gaussian responses, including binary, count, and compositional outcomes. The proposed framework jointly estimates the contemporaneous conditional dependence network and the lag-one temporal dependence structure. We further develop a penalized Monte Carlo expectation maximization procedure for estimation. Simulation studies across multiple graph structures and outcome types show that the proposed method can recover both dependence structures, with more stable recovery for the lag-one temporal dependence structure and more structure-dependent recovery for the contemporaneous conditional dependence network. We apply the proposed method to longitudinal data from the Integrative Human Microbiome Project (HMP2) using a multinomial observation model. The results reveal persistent lag-one temporal dependence and contemporaneous conditional dependence among dominant microbial taxa, demonstrating the utility of the proposed framework for dynamic network analysis of longitudinal non-Gaussian biomedical data.

The rapid growth of technologies such as wearable devices and high-throughput sequencing has generated unprecedented volumes of longitudinal data, often with multiple correlated features. This trend highlights the need for robust clustering methods tailored to multi-feature longitudinal datasets. Although previous reviews have benchmarked clustering methods in this context, most evaluations relied heavily on simulated and overly simplified datasets, which substantially limit their interpretability and applicability to real-world data. Missing data is also ubiquitous in longitudinal data, yet understudied in the context of clustering. To address these gaps, we systematically evaluate clustering methods across three distinct types of real-world longitudinal datasets: (1) a high-dimensional, low-timepoint metagenomic study, (2) a low-dimensional, high-timepoint multi-site repeated-measure study, and (3) a high-dimensional, high-timepoint repeated-measure study. We further investigate how missing data and different imputation strategies influence clustering performance. In addition, we introduce a novel resampling framework based on the empirical cumulative distribution function (eCDF) and a copula. This framework enables us to simulate multi-feature longitudinal data that preserves the correlation structure of real datasets, while also allowing for controlled variations such as sample size. Our simulation framework naturally accommodates missing data, thereby providing a flexible tool for benchmarking and comparing clustering methods in more realistic longitudinal settings.

Understanding the causal pathways linking the gut microbiome to downstream biomarkers and clinical outcomes is central to microbiome research. When the microbiome acts as the exposure, the high dimensionality and compositional structure of the data introduce unique methodological challenges for mediation analysis. To address these challenges, we propose a latent variable mediation framework that captures variation in microbial composition through a microbial balance, defined as the log-ratio between two unknown subsets of taxa. This balance serves as a latent scalar exposure and simplifies the estimation of the overall indirect effect at the community level, while simultaneously identifying specific taxa that contribute to the overall direct and indirect effects. We describe the model's estimation and inference, and illustrate the method using data from the Multiethnic Cohort--Adiposity Phenotype Study, where we examine the mediation pathway from the gut microbiome through lipopolysaccharide-binding protein, a marker of metabolic endotoxemia, to percent liver fat.

Somatic mutations in sequenced tumors are often detected with imperfect sensitivity. They can be missed because tumor heterogeneity results in some mutations being absent in a sample at hand or present at very low frequencies which, combined with technical limitations like insufficient sequencing depth, reduces their detectability. It becomes apparent with multi-regional sequencing that some mutations confined to specific areas can be overlooked in single-region analyses. Adjusting for such sensitivity becomes critical when mutation frequencies between groups of subjects are compared, such as responders and non-responders to a given therapy. Otherwise, the estimate of the mutation prevalence, or of the difference in prevalences between the groups, can be biased, leading to inflated type I error. Here, we propose a weighted approach to adjust for imperfect sensitivity of the mutation detection and compared it to the unweighted approach. We assume that a surrogate of mutation sensitivity is observed only among patients with called mutation, and we derive its unconditional mean. Such adjustment leads to less biased estimates of mutational frequencies and better control of the type 1 error.

Personalized medicine aims to tailor treatments to individual patients, but its implementation is limited by scarce clinical resources and heterogeneous treatment effects. Electronic health records now offer longitudinal data that can support data-driven therapy selection. We focus on SGLT2 inhibitors, medications with broad metabolic and cardiovascular benefits, whose clinical impact varies significantly based on individual patient profiles. We propose a counterfactual Structural Causal Model (SCM) to estimate individual treatment effects (ITE) of SGLT2 inhibitors from observational data. The model assumes counterfactual stability and uses a Gumbel-Max sampling mechanism to capture stochastic treatment assignment while enabling differentiable inference. Optimal treatments are defined by maximizing the counterfactual value of future diagnostic outcomes. We evaluate the Counterfactual Gumbel-Max SCM using a doubly robust value estimator and show that, across diagnostic categories, it outperforms real-world prescriptions in predicted treatment value, highlighting subpopulations where SGLT2 inhibitors yield superior outcomes. By integrating causal inference with electronic health data, this framework provides individualized counterfactual reasoning about SGLT2 treatment effects and can support clinicians in identifying optimal therapies for type 2 diabetes and comorbid conditions.

Sample-specific molecular networks, which capture how an individual's gene regulatory structure may deviate from the population as a whole, are important tools for understanding heterogeneity in health and disease. However, most existing methods either do not accommodate uncertainty quantification or cannot handle missing data, non-Gaussian modalities, and longitudinal structure common in modern multi-omic studies. For example, copy number alterations in cancer genomic studies are non-Gaussian and often only partially observed, and their regulatory effects on gene expression may vary substantially across patients, variation that a single, population-averaged network cannot capture. To address these challenges, we propose PRISM (Personalized Relational Inference via Sample-specific Modeling), an empirical Bayes framework for estimating sample-specific covariance matrices from cross-sectional and longitudinal multi-omic data. PRISM models each individual sample's covariance matrix as drawn from an inverse-Wishart prior whose parameters are estimated from the data and allowed to vary across subpopulations or evolve over time, yielding a closed-form posterior and edge-level uncertainty quantification, ensuring scalability to high-dimensional omics data. In simulation experiments, PRISM achieves point-estimate accuracy comparable to or better than existing methods, substantially improved uncertainty calibration under non-Gaussian marginals, and lower computational cost, while maintaining robust coverage under moderate missingness. We apply PRISM to a yeast cell-cycle time course dataset, recovering time-varying biologically meaningful regulatory relationships independently documented in literature, and to paired copy number and gene expression data from ovarian cancer patients, identifying copy-number-expression interactions associated with survival, including TFs HNF4G and PAX8 which have independent evidence of clinical relevance in this cancer.

Estimating population-level dependence between brain regions is challenging when the number of regions greatly exceeds the number of individuals. We introduce an information-borrowing generalized Bayesian method for high-dimensional covariance estimation from region-level brain measurements. A divide-and-conquer parameterization separates marginal variances from correlation calibrations: an empirical \emph{pilot} correlation matrix--by default, obtained by standardizing a sample covariance matrix (SCM)--provides a reference structure that is adaptively calibrated through region-specific parameters shared across pairwise correlations, reducing the number of free parameters from $\mathcal{O}(d^2)$ to $\mathcal{O}(d)$ ($d$ is the number of regions) while guaranteeing a positive-definite covariance through the inclusion of a ridge term. Horseshoe priors on the calibration parameters enable selective borrowing and shrinkage of pilot information. Under suitable conditions on pilot accuracy, we establish theoretical guarantees for covariance estimation while retaining scalability in ultra-high-dimensional settings via a low-rank likelihood representation and a Woodbury-accelerated Hamiltonian Monte Carlo algorithm. Simulations showed lower estimation error than the SCM, bagging and iteratively refined principal orthogonal complement thresholding (POETRY). Applications to a semi-synthetic cortical-thickness and ADHD resting-state fMRI data revealed significant group differences in the default-mode--visual network and somatomotor network--limbic network blocks after Benjamini--Hochberg correction.

Selection bias is a common concern in causal inference. Two mechanisms can give rise to selection bias: collider bias or stratifying on an effect modifier of the exposure-outcome causal effect. Research into the impact of selection bias on instrumental variable (IV) analyses has focused primarily on collider bias. Here, we show that IV estimates can be biased if selection is an effect modifier (or caused by an effect modifier) of the exposure-outcome effect. In a simulation study, we investigate the size and direction of the bias when estimating effects on the mean difference or odds ratio scale. Bias can also occur due to heterogeneity in the instrument-exposure relationship, although its magnitude is typically smaller. Finally, we discuss methods to deal with this bias in practice, including inverse probability weighting (IPW).

For the purpose of differential expression (DE) analysis in single-cell RNA-sequencing (scRNA-Seq), phenotype differences between samples are often tested within specific cell types. Cell types are regularly imputed by clustering the same gene expression data which is later used for phenotype testing. This creates the potential for a "double-dipping" or post-selection inference problem resulting in inflated rates of false discoveries. While this selection bias is known to inflate significance in cell-type marker identification, its effect on sample-level phenotype testing, e.g. in patient cohorts, has never been explored despite the growing preponderance of this type of analysis. To address this, we perform an extensive simulation study and demonstrate that naive clustering on uncorrected embeddings can severely inflate the False Discovery Rate (FDR) in the presence of strong phenotypic differences. However, we further show that applying batch-correction methods to remove phenotypic effects prior to clustering resolves the FDR inflation with no obvious loss of power. Finally, we provide measures of phenotypic imbalance that can be applied to real datasets which closely track the false discovery proportion and thus can be used to as part of data exploration to gauge the risk of post-selection inflation of p-values in a particular dataset.

Voxel-level prediction of treatment response from longitudinal fluorodeoxyglucose (FDG) positron emission tomography/computed tomography (PET/CT) can enable spatially adaptive dose planning in advanced non-small cell lung cancer (NSCLC), but point predictions alone lack the uncertainty quantification needed for safe clinical decision-making. Standard conformal prediction (CP) provides distribution-free coverage but yields spatially uniform intervals and assumes exchangeability, an assumption violated by intra-tumoral spatial correlation. We propose a multiscale CP framework built on a re-optimized VoxelForecast generalized least squares model that explicitly captures voxel-level covariance. Variogram sensitivity analysis identified the Stable model as significantly superior to the previously used Matern family. The framework combines hierarchical nested leave-one-patient-out calibration, which restores exchangeability for valid voxel-level coverage, with residual-variance CP (VarCP), which learns a heteroscedastic, location-dependent scale to produce adaptive intervals while preserving finite-sample marginal coverage guarantees. On two prospective NSCLC trial cohorts, FLARE-RT (N=25 locally advanced, chemoradiotherapy; 11,100 voxels) and BRIGHT (N=19 metastatic, chemoimmunotherapy; 26,980 voxels), multiscale feature integration reduced voxel-level prediction error by 35% and 31%, respectively. VarCP maintained nominal coverage while producing 7-18% (FLARE-RT) and 4-8% (BRIGHT) narrower voxel-level intervals than a hierarchically calibrated standard CP baseline (paired Wilcoxon signed-rank, significant at most miscoverage levels), with comparable lesion-level gains (7-18% and 6-10%) via median-based aggregation of voxel residuals. VarCP delivers tighter yet mathematically guaranteed spatial uncertainty maps, supporting both voxel-level dose adjustment and whole-lesion treatment planning in advanced NSCLC.

Many biological quantities in omics are inherently periodic or directional, including circadian phase, cell-cycle position, and cellular orientation. Treating such quantities as ordinary linear variables can introduce artificial discontinuities and obscure biological interpretation. Circular methods are well developed across directional statistics, bioinformatics, and spatial statistics, but practical guidance for omics applications remains limited. We present key considerations for circular data analysis in single-cell and spatial omics, including the provenance and observational unit of the circular quantity, the inferential or predictive objective, the distinction between tissue composition and residual spatial dependence, and the choice of appropriate models. We describe relevant approaches for circular representation, inference, clustering, regression, machine learning, and spatial analysis, together with representative R software. Their practical use is illustrated through spatial transcriptomic case studies of derived circadian peak phase in a mouse model of Alzheimer's disease and latent cell-cycle position in human ovarian cancer. These examples demonstrate how the definition of the angle, its observational unit, and tissue heterogeneity shape appropriate analysis and interpretation, including when spatial modeling is not warranted. We conclude with practical recommendations and directions for future methodological development.

Multiview data analysis has become increasingly important for integrating heterogeneous information from multiple data sources. Cooperative learning provides a flexible framework for multiview integration by bridging early and late fusion strategies. However, its use of squared loss makes it sensitive to outliers. In this study, we propose robust cooperative learning based on the soft-IPOD framework. The proposed method introduces an observation-specific outlier absorption vector into the cooperative learning objective, thereby reducing the influence of outlying observations while preserving information sharing across views. The resulting optimization problem can be formulated as a weighted Lasso-type problem using an augmented design matrix and can be naturally extended to more than two data views. Numerical simulations and a real-data application to labor onset prediction demonstrate the predictive performance of the proposed method in the presence of outliers.

Single cell RNA sequencing (scRNAseq) has revolutionized cellular biology over the past decade. At its most basic, a scRNAseq dataset is a cell by gene matrix with entries giving the transcription count of a particular gene in a particular cell. A common workflow applies log-normalization and z-scoring to the count matrix and then embeds the cells in a lower dimensional space using principal component analysis (PCA). A large number of computational tools use the PCA embedding as a starting point, but statistical models for the distribution of cells in the PCA space are lacking. As a result, many tools depend on ad-hoc assumptions that likely lead to signal loss. In this work, we address this gap using the metacell model, a statistical model of count matrices introduced by previous authors. We show that application of the workflow to count matrices generated under a metacell model yields PCA embeddings distributed as a Gaussian mixture model (GMM). To our knowledge, this is the first result connecting a statistical model of counts to the distribution of cells in PCA space. Our derivation depends on methods from random matrix theory (RMT) and we extend and establish connections between RMT and scRNAseq count matrices. We then use the predicted GMM to test the metacell model on 11 datasets. The predictions are accurate for counts generated by the metacell model but underestimate the variance of real cells for most of our datasets. We show that this error largely arises from gene correlations within metacells, which the metacell model does not capture. Some of our datasets are fit relatively well by the metacell model, suggesting that the fit of the metacell model to true datasets depends on sufficiently dense sampling. Finally, we show that the strength of noise in the PCA space varies significantly across metacells, suggesting that statistical modeling of PCA embeddings can improve existing downstream workflows. Our results provide a framework for deriving and testing statistical models of PCA embeddings from statistical models of counts and identify a limitation of the metacell description for scRNAseq datasets.

Researchers must often rely on observational data to draw conclusions about the causal relationship between an exposure and an outcome. This is common in environmental epidemiology research, where the assignment to an exposure can be unethical or impossible. Consider a setting in which we seek to estimate excess mortality during heat waves, but are unable to randomize individuals to moderate or extreme heat. We propose to follow a causal framework for non-randomized time series that uses what we call historical time series matching to construct hypothetical randomized experiments, in which background covariates are balanced. We applied this framework and assessed the following causal question: Are daily mortality counts among elderly Americans higher during hot versus warm summer time periods? We addressed this causal question by drawing frequentist, Fisherian, and Bayesian inferences for five American cities. We estimated the average treatment effect of heat on three-day mortality with 95% Fisherian intervals to be 0.5 [-3.6, 4.4] deaths in Pittsburgh, 3.7 [0.7, 6.6] deaths in Seattle, 12.5 [6.1, 18.9] deaths in Chicago, 8.5 [1.4, 15.6] deaths in Los Angeles, and 20.3 [11.2, 29.5] deaths in New York.

Transcriptome-wide association studies (TWAS) provide an effective framework for identifying genes associated with complex traits. Population-scale single-cell transcriptomic data enable genetically regulated expression (GReX) prediction and TWAS analyses at cell-type resolution, but the predictive performance of existing single-cell TWAS methods remains limited. Here, we develop MOD-scTWAS, a module-based method that jointly models GReX for genes within co-expression modules to borrow information across genes. Starting from a generative model for single-cell gene expression, MOD-scTWAS accounts for the heteroscedasticity and cross-gene correlation of individual-level pseudobulk expression in joint GReX prediction. In cross-validation analyses of the OneK1K dataset, MOD-scTWAS achieved higher mean GReX prediction accuracy than scTWAS across all 14 cell types and increased the number of imputable genes. When applied to TWAS analyses of UK Biobank quantitative hematological traits, MOD-scTWAS identified more significant cell type-gene-trait associations than scTWAS. These results demonstrate the potential of leveraging gene co-expression through joint modeling to improve cell-type-specific GReX prediction and TWAS discovery.

Background: Induced pluripotent stem cell (iPSC)-derived kidney organoids hold great promise for disease modeling, drug screening, and regenerative medicine. However, since these organoids are synthetic constructs, it is essential to evaluate how faithfully their constituent cell populations recapitulate the corresponding cell types in the native kidney. Single-cell transcriptomic label-transfer tools, which assign cell-type labels by integrating cells with well-annotated reference data and classifying them, are increasingly used for this purpose. Yet the reliability of the resulting labels - particularly for immature or partially differentiated cells that might have no clear in vivo counterpart - has not been systematically addressed. Methodology: To quantify the reliability of cell-type labels assigned by label-transfer algorithms, we introduce two complementary robustness measures for each cell type - the classification stability score, which assesses label stability under simulated Poisson noise, and the prediction specificity score, which assesses how well the predicted cell-type label probabilities match the assigned labels. We apply this methodology to scRNA-seq datasets from the human fetal kidney and four iPSC-derived kidney organoid protocols, using a mouse fetal kidney dataset as a reference. Results: We find that label-transfer robustness scores are generally high in most human fetal kidney cell types, except for cells of immature epithelial structures. In the organoid-derived cell subpopulations, robustness scores are high for the cap mesenchyme but vary for epithelial lineages according to each protocol's specific differentiation strategy. These observations are consistent across both classification stability and prediction specificity scores, which are generally well correlated with each other. We further validate our approach through a simulated cell-removal experiment, confirming that the proposed measures correctly identify unreliable classifications. Significance: These robustness measures provide a generalizable framework for assessing the reliability of cell-type label transfer in single-cell transcriptomic datasets, and can be extended to other biomedical settings that rely on automated single-cell classification.

Electronic health records provide a source of real-world information for disease risk prediction, but many machine learning models learn implicit associations that are difficult to inspect or use for clinical reasoning. Causal Bayesian networks (CBNs) offer a form of causal-interpretable machine learning by representing conditional dependencies, putative directional relations, and probabilistic evidence propagation in a directed acyclic graph. However, conventional CBN structure learning becomes unstable and computationally expensive when applied to large-scale EHR data. To address this limitation, we propose UPEBNL, a scalable framework for causal-interpretable cancer prescreening based on parallel CBN learning. UPEBNL integrates adaptive data slicing, quality-aware structure aggregation, and global DAG construction to learn stable and interpretable dependency structures from large-scale observational data. We evaluated UPEBNL in high-dimensional and multi-million-sample simulations and applied it to EHR-based prescreening for esophageal and colorectal cancer. In simulations, UPEBNL improved structural recovery accuracy by nearly 40% and achieved up to a 221.28-fold speedup over conventional CBN learning strategies. For cancer risk prediction, the learned CBNs provided interpretable evidence paths and achieved validation AUCs of 0.8171 for esophageal cancer and 0.784 for colorectal cancer. Calibration and decision curve analyses further supported the reliability and clinical utility of the models. These findings suggest that scalable CBN learning can support interpretable cancer prescreening from large-scale EHR data.

In hypertension screening, previously diagnosed cases provide reliable positive labels, whereas individuals without a prior diagnosis comprise both undiagnosed hypertensive and non-hypertensive individuals, creating a canonical positive-unlabeled (PU) learning setting. We propose PU-Boost, a two-stage reliable-sample reconstruction framework that uses only readily obtainable demographic, behavioral, medical-history, and anthropometric variables. During training, only previously diagnosed hypertension cases are treated as labeled positives (P), while all remaining individuals are treated as unlabeled (U); no additional reference labels are supplied to the two-stage reconstruction procedure. Elastic Net first provides a global ranking of P versus U, after which a random forest refines the uncertain region and highly ambiguous observations are rejected. XGBoost is then trained on the reconstructed class-balanced pseudo-labeled sample. In the independent test set, PU-Boost achieved a sensitivity of 0.817 and a balanced accuracy of 0.682, the highest observed values among the compared models. It produced 729 false negatives, 41.1% fewer than standard XGBoost, although specificity decreased to 0.547. SHAP analysis indicated that established hypertension-related features, including age, body mass index, waist circumference, and family history, contributed strongly to model output. PU-Boost therefore provides a statistically motivated framework for reducing missed cases when screening for undiagnosed hypertension under incomplete diagnostic labeling.

Understanding when and how pathogens spread from person to person is fundamental to characterizing transmission, interpreting surveillance data, and designing effective control measures. Yet key epidemiological delays -- including the latent period, incubation period, generation interval, and serial interval -- are difficult to estimate because key events such as infection, infectiousness, symptom onset, and transmission are rarely observed directly. Current approaches typically estimate these delays using flexible statistical distributions rather than explicitly linking them through the underlying processes that generate them, making it difficult to interpret the underlying biology and to ensure consistency across related delay estimates. Here, we develop a stochastic event-time framework built around a dynamic within-host pathogen-growth model that treats epidemiological delays as coupled outcomes of a shared biological process. By representing these key events as linked stochastic events, the framework produces internally consistent delay distributions while inferring parameters with mechanistic interpretations conditional on the assumed within-host model that govern pathogen growth, clearance, infectiousness, and symptom onset. The framework also accommodates censored and partially observed infection and symptom-onset times through a unified observation model. Across simulation studies, the framework accurately recovered generation-interval characteristics across sample sizes and disease time scales. In applications to SARS-CoV-2 and mpox transmission-pair data, the framework showed improved predictive fit relative to conventional lognormal delay models in several analyses, although the magnitude of improvement varied across datasets and observation models; it estimated shorter generation intervals for Omicron than Delta and revealed alternative within-host parameters compatible with the mpox transmission data. This mechanistic framework provides a biologically interpretable alternative to independently fitted delay distributions, yields internally consistent estimates of epidemiological delays, and establishes a foundation for integrating biological and epidemiological data to better understand, predict, and control pathogen transmission.