Biometrics

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.

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.

In competing risks settings, covariate effects and group comparisons are usually assessed one event at a time - through log-rank or Cox tests on the cause-specific hazards, or Gray's test or Fine-Gray regression on a cumulative incidence function (CIF). This can obscure a clinically important quantity: the ratio between the event of interest and the competing event, since groups may differ little on the individual events yet differ sharply in their ratio. The generalized competing event (GCE) framework makes this ratio the object of inference; on the cause-specific scale the hazard ratio omega+(t) = lambda_1(t)/lambda_2(t) is estimated efficiently from a single stacked (Lunn-McNeil) model. We extend the framework to two scales that describe realized incidence. The subdistribution hazard ratio omega-tilde+(t) = lambda-tilde_1(t)/lambda-tilde_2(t) is estimated by a stacked, risk-set-weighted extension of the Lunn-McNeil construction; the cumulative-incidence ratio rho(t) = F_1(t)/F_2(t) - the odds that a subject's realized event by time t is the event of interest - by jackknife pseudo-observation regression of the Aalen-Johansen estimator. We relate the three contrasts: rho equals omega+ exactly under proportional cause-specific hazards, and equals omega-tilde+ only in the small-time limit under proportional subdistribution hazards, drifting toward 1 thereafter. The orthogonality that makes omega+ efficient is lost on both cumulative-incidence scales - omega tilde+ through overlapping weighted risk sets and shared censoring weights, rho through the shared all-cause survivor - so each carries a covariance term that must be handled and that bounds efficiency relative to the hazard-scale test. We derive the corresponding variances, study operating characteristics by simulation, illustrate on hypothetical prostate and head-and-neck cohorts, and provide an implementation in the gcemod R package.

In biomedical research, multiple types of high-dimensional data, such as genomic, transcriptomic, proteomic, and metabolomic data, are increasingly collected from the same subjects. Integrating these multiple data views can improve prediction by exploiting shared or complementary information across the views. Cooperative learning provides an agreement-based framework for multiview supervised learning by encouraging predictions obtained from individual views to be similar. However, the original framework assumes independent observations and therefore does not account for clustered structures, such as repeated measurements obtained from the same subject. To address this limitation, we propose Cooperative Learning with a penalized Linear Mixed Model (CL-pLMM) for high-dimensional multiview data with a clustered structure. CL-pLMM replaces the ordinary prediction loss in cooperative learning with a covariance-weighted loss that accounts for within-cluster dependence, while retaining the agreement penalty between views and a Lasso penalty for variable selection. We further show that its objective function can be represented as a penalized linear mixed-effects model applied to augmented data, allowing existing estimation procedures to be used. The performance of CL-pLMM is evaluated through simulation studies under various signal and dependence settings and an application to longitudinal proteomic and metabolomic data for predicting the time to spontaneous labor.

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.

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.

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).

Constructing microbial association networks is a common strategy for exploring relationships among taxa in microbiome studies. Although marginal correlation methods are easy to implement and allow formal inference, they can produce spurious edges driven by indirect associations through other taxa. Conditional graphical-modeling methods aim to recover direct associations, but many rely on Gaussian or linear assumptions and often provide limited uncertainty quantification. We propose a conditional, nonparametric approach based on the scaled expected conditional covariance (SEcov). SEcov measures population-level conditional association by residualizing each taxon with respect to the remaining taxa and scaling the resulting expected conditional covariance. The resulting estimator can incorporate flexible machine-learning methods for conditional-mean estimation and admits asymptotic normal inference, enabling p-values and confidence intervals for taxon-pair associations. We demonstrate through simulation studies that our proposed approach improves network recovery relative to other methods, and we illustrate the new method via construction of a co-occurrence network for the vaginal microbiome during pregnancy. IMPORTANCEHigh-throughput sequencing has made it possible to characterize microbial communities at large scale, and network analysis is widely used to summarize relationships among taxa. However, networks based on marginal correlations may include indirect associations, whereas many conditional graphical models rely on assumptions that may be difficult to justify for sparse, zero-inflated, compositional microbiome data. SEcov offers a practical alternative by estimating conditional associations nonparametrically and attaching inferential uncertainty to individual edges. This allows investigators to construct microbiome networks using statistically interpretable evidence for taxon-pair associations, rather than relying solely on arbitrary correlation cutoffs or regularization tuning parameters.

Abstract: Compartmental epidemic models conventionally treat the probability of moving between dis ease states as fixed over time, an assumption that sits uneasily with the reality of a pandemic in which lockdowns, mask mandates, vaccination roll-out, and the arrival of new variants continually reshape transmission. This paper develops a time-inhomogeneous Markov chain framework for the Susceptible-Exposed-Infectious-Removed (SEIR) process, in which each transition probability pab(t) is allowed to vary with calendar time while respecting the struc tural zeros implied by the SEIR compartmental flow. We derive the constrained maximum likelihood estimator of pab(t) under these structural constraints, establish its finite sample efficiency, asymptotic normality, and Wilson score confidence intervals, and construct a like lihood ratio test of the null hypothesis that a compartments exit probability is constant over time. We further propose a stochastic machine learning hybrid extension in which the raw, kernel smoothed transition probabilities are regressed on policy and mobility covariates using both a logistic generalized linear model and a random forest, allowing the framework to attribute time-inhomogeneity to observable interventions. The methodology is applied to a compiled daily, district level COVID-19 surveillance panel for Sierra Leone spanning March 2020 to December 2023 (16 districts, 1,401 days). The likelihood ratio test rejects time-homogeneity of the exposed to infectious transition in 15 of 16 districts and of the infectious-to-removed transition in 8 of 16 districts ( = 0.05), and the covariate augmented logistic model achieves an out of sample Brier score roughly 76 times smaller than a time homogeneous pooled baseline, with healthcare capacity and the time trend emerging as the most influential predictors in the random-forest component. These results provide statisti cal evidence that time-inhomogeneous, covariate informed Markov models offer a materially better description of district-level COVID-19 transmission in Sierra Leone than classical time homogeneous compartmental models, with implications for sub-national outbreak monitoring in resource constrained settings

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.

Spatial transcriptomics technologies enable measurement of gene expression while preserving spatial tissue organization, but they remain highly sensitive to technical variability such as library size differences, slide-level effects, and spatial artifacts. Most existing normalization approaches treat normalization as a preprocessing step and perform downstream analyses on normalized values as fixed inputs, ignoring the uncertainty introduced during normalization. We introduce GPSNorm (Gaussian Process Spatial Normalization), a Bayesian spatial normalization framework that jointly models technical variation, spatial structure, and biological signal within a unified hierarchical model. GPSNorm represents gene expression counts using a negative binomial latent Gaussian model whose spatial component is a Gaussian Markov random field approximating a Gaussian process and performs efficient approximate Bayesian inference using the Integrated Nested Laplace Approximation (INLA), producing posterior estimates that propagate normalization uncertainty into downstream differential expression analysis. In simulations anchored to empirical spatial transcriptomics data, GPSNorm accurately recovers spatial technical structure and improves log-fold change estimation compared with existing normalization methods. Applications to three spatial transcriptomics datasets--including human dorsolateral prefrontal cortex Visium data, a GeoMx COVID-19 lung damage study, and the Spatial Organ Atlas kidney dataset--demonstrate improved preservation of biologically expected spatial patterns and marker gene contrasts. These results show that jointly modeling normalization and downstream inference can improve the robustness and interpretability of spatial transcriptomics analyses.An open-source R implementation of GPSNorm is available at https://github.com/Tiny-Quant/GPSNorm.

Accurate normalization is essential for differential expression analysis of RNA-sequencing data. Popular normalization methods such as the median-of-ratios and trimmed mean of M-values do not leverage information from the experimental design. This may be inefficient in experiments with large-scale systematic expression changes or complex designs. Here, we introduce design-informed size factor estimation (disize), a normalization method that uses information from the experimental design to improve accuracy. disize uses a modified generalized linear mixed model to robustly distinguish between biological signal and sample-specific size factors. We also propose a mechanistically justified data-generating process for RNA-sequencing counts that is derived from previous models of transcription and sequencing. Through simulations based on this data-generating process and validating on true RNA-seq data, we show that disize recovers size factors more accurately than existing methods, particularly in challenging scenarios with low gene expression and a high proportion of differentially expressed genes; this in turn improves downstream analysis. disize provides a robust and accurate approach to normalization, highlighting the significant benefits of integrating experimental design information directly into normalization for transcriptomic datasets. Author summaryIn transcriptomic analysis, normalization adjusts for technical biases arising from library preparation and sequencing. Methods implemented in widely used packages like DESeq2 and edgeR ignore information in the experimental design during normalization. Incorporating information from the experimental design into a normalization method has the potential to yield more accurate results. To do this, we developed a new method, design-informed size factor estimation (disize), that uses a statistical model to jointly account for the biological signal defined by the design and the sample-specific batch effect. By separating the biological variation into its components, disize can more robustly estimate the batch effect. To validate our approach, we constructed a flexible simulation framework relying on a mechanistically justified data-generating process for RNA-seq data. Our benchmarks on both simulated and true RNA-seq data show that disize recovers the true size factors more accurately than existing methods, particularly in challenging scenarios with low counts or a high proportion of differentially expressed genes. This improved normalization yields more reliable downstream results in differential expression analysis.

Estimating differential gene expression is a common task in RNA-seq. Current methods mostly rely on normalization to reduce variance and increase accuracy. These methods are widely used and provide invaluable information about transcriptomic changes between biological conditions. However, widely used normalization methods are known to distort estimates of transcript differential expression when a majority of genes are upregulated or downregulated, or when the total RNA content per cell changes. Despite the presence of global expression shifts in a number of contexts, few methods exist that can provide accurate normalization and estimate linear models under this context without spike-in controls. Here, we present \textbf{GLORB} (\textbf{G}\textbf{L}\textbf{O}bal-shift \textbf{R}obust \textbf{B}ayesian model), a method for estimating generalized linear models under global upregulation. We develop two models that are able to recover differentially expressed genes and linear model coefficients with lower distortion of results. We show that our model's method of accounting for library size variance is consistent with DESeq2's median of ratios, and edgeR's trimmed mean of M-values under conditions when a minority of genes are upregulated and outperforms them under circumstances when most genes are either increased or decreased between groups. Finally, because our method does not rely on calculating geometric means for each gene it is able to work in datasets with much higher sparsity.

Selecting stochastic gene-expression models from single-cell counts requires balancing goodness of fit against unnecessary mechanistic complexity. The probability-generating-function-based Bayesian information criterion (PGF-BIC) combines covariance-weighted fitting in generating-function space with a complexity penalty, allowing candidate models to be compared without reconstructing their full count distributions. However, its conventional zero-threshold rule does not account for sampling uncertainty in the fitted score difference and may therefore favor overly complex models in finite samples. To address this limitation, we develop an uncertainty-aware PGF-BIC rule that selects the more complex model only when its score advantage exceeds a data-driven threshold. We use Cantelli's one-sided inequality to motivate a selection margin expressed in terms of a standard deviation. To determine this scale, we use influence functions to quantify sensitivity to small perturbations in the data distribution and obtain a first-order description of sampling fluctuations. The resulting variance estimate accounts for variability in both the empirical probability generating function and the estimated covariance weights, yielding a data-driven threshold for assessing the complex model's score advantage. A Poisson versus Bursty benchmark shows that the calibrated rule reduces incorrect selection of the more complex model. The calibration requires neither resampling nor additional optimization, incorporating sampling uncertainty into model selection while retaining the computational efficiency of PGF-BIC.

Antimicrobial-resistant organisms (AMROs) impose a major burden on healthcare systems, yet routine surveillance cannot readily distinguish colonization imported at admission from transmission acquired within hospitals. This gap is especially consequential because both processes may vary substantially across wards, while asymptomatic carriage, incomplete testing, imperfect diagnostic sensitivity, and patient movement obscure the underlying transmission dynamics. To address this challenge, we developed a blockwise agent-based iterated filter (BAIF) for inference in a patient-level transmission model on a dynamic ward co-location network. The model tracks susceptible and colonized patients as they move across wards, represents unobserved colonization histories, and incorporates the recorded testing schedule and imperfect diagnostic sensitivity. BAIF uses blockwise likelihood evaluation and resampling to estimate ward-block-specific transmission rates and importation probabilities in this high-dimensional latent system. Synthetic experiments showed that BAIF recovered these parameters from partially observed outbreaks. We then applied the framework to hospitalization and microbiological surveillance data collected from 2012 to 2016 at an urban quaternary care hospital in New York City for four AMROs. Transmission and importation were highly heterogeneous across ward blocks. Elevated transmission was repeatedly concentrated in the same ward groups, whereas blocks with the highest importation varied by pathogen. By distinguishing importation-dominated from transmission-dominated ward blocks, the framework can inform more targeted surveillance and infection-control strategies. More broadly, BAIF provides an effective inference framework for high-dimensional, partially observed agent-based models on dynamic contact networks.

Background: Descriptive mapping of intensive care unit (ICU) and high care unit (HCU) capacity across Japan's secondary medical areas (SMAs) characterizes where beds exist, but medical planning also requires answers to prospective questions: how likely is a capacity shortfall under demand surge, which assumptions drive that risk, and how much protection do inter-zone transfer arrangements provide. No openly available tool addresses these questions at the SMA level, the geographic unit at which Japanese medical plans are written. Methods: We developed MeshScope-Scenario, a probabilistic capacity-demand framework operating on the MeshScope-Region platform. For a selected SMA, observed inputs (notified ICU/HCU beds from the Hospital Bed Function Reports; resident population) are combined with four explicitly flagged assumption parameters - effective staffed-bed rate, concurrent severe-care demand per 100,000 population, surge multiplier, and net cross-boundary inflow - each with a user-specified distribution. A seeded Monte Carlo engine (deterministic reproduction under a fixed seed) estimates the distribution of bed shortfall; interventions are compared under common random numbers. Parameter dependence is introduced by a Gaussian copula with automatic positive-semidefinite correction; global sensitivity is quantified by Sobol first-order and total-order indices (Saltelli sampling, Jansen estimators) alongside a deterministic one-at-a-time tornado analysis. A two-zone extension transfers unmet demand to the nearest ICU-holding SMA using road-network travel times measured in MeshScope-Region, with a transfer time limit and an acceptance cap; because both zones share the same systemic draws, correlated exhaustion of donor capacity under surge ("shared-fate" risk) is represented structurally. A seasonal layer applies twelve monthly surge multipliers and reports the distribution of annual maximum shortfall and month-specific shortfall probabilities. The engine is a dependency-free pure-function module verified by 40 statistical tests. Results: The framework reproduces identical output under identical seed and input; a flat seasonal profile reproduces the non-seasonal model exactly; copula factorization error is below 1e-9; and 10,000 iterations across three intervention variants complete in approximately 50 ms in a standard browser, permitting fully interactive use. Applied to three archetypal SMAs from the observed FY2024 supply map (seed 42, 20,000 iterations, demand prior 5 per 100,000), an ICU-zero zone with a transfer partner 48 minutes away has shortfall probability 66.0% (P50 4.8, P90 20.7 beds); a median metropolitan zone, 32.7% (P90 5.9), with Sobol indices ranking demand density and surge dominant; a high-supply zone shows zero shortfall up to approximately 2.5x surge. For the ICU-zero zone, independent-donor reasoning credits the transfer arrangement with a 2.27-bed reduction in expected shortfall, of which shared-fate correlation removes 93%; the probability of severe shortfall under the arrangement equals that with no arrangement at all, while a committed pool of five donor beds (6% of donor effective supply) restores a 13-point reduction and more than halves the arrangement's correlation exposure. Conclusions: MeshScope-Scenario extends SMA-level capacity mapping from description to prospective risk assessment. All demand-side inputs are declared assumptions with adjustable distributions rather than estimates presented as fact; the framework's value is to make the consequences of those assumptions, and their interaction with observed supply, explicit, reproducible, and inspectable for planning deliberation.

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.

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.

Causal mediation analysis is widely used to identify biological pathways linking exposures to outcomes, but most methods assume homogeneous mediation effects across individuals. In high-dimensional omics settings, this assumption can mask important heterogeneity driven by demographic, genetic, or environmental factors. We propose the M-high-learner, a flexible framework for detecting heterogeneous mediation effects with high-dimensional mediators. The method identifies mediators with subgroup-specific indirect effects while distinguishing them from null or homogeneous signals and controlling the type I error rate. It is computationally efficient, scalable, and yields interpretable sub-types. Simulation studies show that the proposed approach achieves high power while maintaining accurate error control. Applications to the Framingham Heart Study and the Multi-Ethnic Study of Atherosclerosis reveal that the mediation role of gene expression in sexs effect on high-density lipoprotein varies across subgroups defined by body mass index and age. Our framework provides a practical tool for uncovering heterogeneous biological mechanisms in high-dimensional genomic studies. Author SummaryBiological processes linking risk factors to disease often differ across individuals, but many existing methods assume these processes are the same for everyone. This can hide important differences between groups. We developed a powerful method to identify when these pathways vary across subgroups using large-scale molecular data. Our approach detects differences in how intermediate biological factors contribute to outcomes in populations defined by characteristics such as age and body mass index. Applying our method to population studies, we found that some biological pathways operate differently across groups, suggesting that key mechanisms may be missed when differences are ignored. Our work provides a tool to better understand how disease-related processes vary across individuals, which may support more targeted and personalized approaches to health research.

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.