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.

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

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.

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.

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.

Importance. Observational studies dominate evidence in cardiovascular outcomes research, yet no widely adopted reporting standard requires authors to disclose how much of an observed association is attributable to systematic error. The negative-control calibration framework of the Observational Health Data Sciences and Informatics consortium returns a binary answer: a calibrated P-value either crosses 0.05 or it does not. That answer conflates "is there a signal" with "is the signal credible." Objective. To define the Bias Attribution Fraction (BAF), a continuous estimator in [0, 1] that quantifies the fraction of an observed log relative risk that is attributable to systematic bias, together with its 95% credible interval, and to evaluate whether reporting calibrated P-value plus BAF plus interval width produces more defensible decisions than reporting the calibrated P-value alone. Design. Three-part study: (i) metric definition and interval construction; (ii) Monte Carlo simulation across 1920 conditions by 1000 repetitions = 1,920,000 simulated target estimates (a primary grid plus a sign-flipped companion grid of identical size, pooled for all analyses), with truth known by construction; (iii) two-arm beta-blocker target trial emulation in MIMIC IV, adhering to the TARGET guideline, with targeted maximum likelihood estimation (TMLE) and 12 negative-control outcomes. Setting. MIMIC IV (critical care) and MIMIC IV ECG, 2008 to 2019. Participants. Adults with a first intensive care unit stay and a recorded 12-lead ECG. Exposures. Two clinical questions are evaluated: (1) reaching at least 50% of guideline-recommended target dose of an index beta-blocker within 7 days of admission, sustained to 12 months, versus lower or no exposure; (2) receiving guideline-directed medical therapy (GDMT) with at least two of three foundational drug classes at adequate dose for 12 months versus fewer. Main Outcomes and Measures. One-year all-cause mortality. Primary analytic outputs: (i) the calibrated P-value from a negative-control empirical null distribution with 12 outcomes; (ii) BAF-hat with a Markov chain Monte Carlo 95% credible interval; (iii) a verdict drawn from a truth-calibrated lookup table mapping the BAF-hat by interval-width cell to a probability of being bias-dominated. Results. Across 1,920,000 simulated repetitions, the calibrated P-value alone yielded a misuse rate (declared usable while truth was bias-dominated) of 36.02%. Adding the BAF point estimate reduced misuse to 16.48%, and requiring the whole credible interval of BAF to lie below 0.5, clear of the bias-dominated zone, reduced misuse to 10.16%. Requiring in addition a calibrated probability of bias dominance below 15% reduced misuse to 3.02% at the cost of declaring 22.1% of estimates usable (versus 61.8% for the calibrated P-value alone). Read from the bias-dominated view, meaning that withholding a repetition is a positive call and being truly bias-dominated is a positive truth, the BAF estimator achieved an area under the receiver operating characteristic curve of 0.922, compared with 0.787 for the calibrated P-value. In MIMIC IV, question 1 produced a calibrated P-value of 0.894 and BAF-hat = 0.875 (95% credible interval 0.385 to 0.992); the apparent one-year mortality reduction of 16.9% was entirely attributable to systematic bias and the association was not usable as effect evidence. Question 2 also failed layer 1, with a calibrated P-value of 0.1026, and BAF-hat = 0.460 placed it in the "competitive, no verdict" zone. Conclusions and Relevance. Pairing the calibrated P-value with BAF and its 95% credible interval turns a binary credibility judgement into a three-channel report (signal existence, effect size, and bias fraction). Despite TARGET-adherent design, doubly robust TMLE, and 12 negative controls, both clinical questions in this study still failed to qualify as effect evidence, with question 1 demonstrating complete bias domination and BAF quantifying the systematic error that remained. The reporting standard proposed here can be applied to any observational study that includes negative controls, without requiring additional data.

Doubly robust survival estimators require nuisance models that the available information can support. We examine design and estimator dependence in TCGA breast cancer data. Treatment indicators disagree for 337 of 1076 patients; 67.4% have unobserved five-year binary outcomes; and receptor imputation expands the triple-negative subgroup from 156 directly measured patients to 226. We compare chemotherapy initiation by day 90 with non-initiation, using adjusted restricted mean survival time (RMST) over the next 1095 days. Weighting-only estimates are 55 days (95% bootstrap interval: -24, 133) in triple- negative disease (16 deaths) and 38 days (12, 72) in the whole cohort (72 deaths). A Cox outcome-model benchmark yields 108 and 35 days, respectively, demonstrating material model dependence in the smaller population. A parameter-to-information ratio summarizes the burden on complete-outcome regressions within cross-fitting. Its simulation-calibrated warning level is 0.31 (0.28 - 0.35) for the original regime; additional simulations examine covariate-dependent censoring, delayed benefit and nonlinear misspecification. The ratio does not certify estimator validity or replace weight diagnostics. The contribution is a reproducible workflow that separates design, nuisance-model support and estimator stability before interpreting observational survival contrasts.

Accurate prediction of individual medical outcomes is essential for optimizing treatment allocation amid rising costs, coverage denials, and limited clinical resources. Traditional predictive models, including linear regression and neural networks, rely on average effects and cannot tailor predictions to the specific circumstances of individual cases. We present relevance-based prediction (RBP), a model-free and task-specific method that predicts outcomes as weighted averages of observed cases, with weights determined by a rigorously defined and theoretically justified measure of relevance. Unlike model-based methods that rely on fixed calibrated parameters, relevance-based prediction revisits the original data for each prediction and customizes both the cases and variables used. We illustrate RBP by applying it to predict opioid treatment outcomes and demonstrate that it provides case-specific insights unavailable from conventional models, including how each prior case informs a prediction, how each variable affects each predictions reliability and value, and how reliable each prediction is before it is made. These individualized insights may prevent misleading average-based decisions and reduce harmful or suboptimal treatment.

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.

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.

Background: Excess mortality estimates depend on an expected baseline, an unobservable quantity typically derived from models with prescribed functional forms and short reference periods that exclude exceptional events. Singular Spectrum Analysis (SSA) offers a nonparametric alternative but requires subjective component selection and is sensitive to outliers. We combined SSA with automated component selection and outlier interpolation to estimate long-term baselines and compile a catalog of excess mortality events (EMEs). Methods: We analyzed weekly all-cause mortality (2015-2024) for six metropolitan areas with contrasting climatic, demographic, and epidemiological contexts: Brussels, Antwerp, and Liege (Belgium), and Athens, Thessaloniki, and Larisa (Greece). SSA was applied to standardized log-mortality rates; components were selected using Monte Carlo SSA, and COVID-19 periods were interpolated iteratively. EMEs were defined as weeks exceeding a quasi-Poisson two-z-score threshold, then screened for coincidence with heat waves, cold waves, COVID-19 periods, and influenza epidemics. Results: Baselines were consistent within countries: annual cycles in Belgium, with an additional significant semiannual cycle in Greece reflecting summer mortality. We cataloged 170 EMEs spanning 456 weeks and 14,210 excess deaths. Pandemic weeks accounted for 77% of excess mortality, and 96 weeks coincided with temperature extremes; 56 EMEs, representing 3.7% of excess mortality, remained unmatched. Conclusions: SSA yields reproducible, low-assumption baselines that capture country-specific seasonality. Because flexible baselines can absorb recurrent seasonal mortality, component selection is a substantive choice that should be reported explicitly, together with its sensitivity.

Amongst measures of linkage disequilibrium, the correlation measure r2, due originally to Sewall Wright, has several advantages. Although basically a haploid statistic, it can readily be estimated using covariance and correlation from diploid data. Its population expectation can also be found from probability methods. We show that the probability of linked identity by descent, LIBD, directly estimates the expectation for r2 for Wright-Fisher models. The simplest such model for two loci is the haploid model, where each haplotype is chosen independently from haplotypes of the previous generation. This gives an equilibrium between recombination (recombination frequency c) and drift (population size N) where the expected value of r2 is approximately 1/(1+4Nc). For diploid models, the expectation is similar for low values of c, but contains a term in c2 and rises to double the haploid value for the limiting value of c = 0.5. An increase in r2 is also expected for the model of monogamous mating, as opposed to random mating for each offspring. The expectation for r2 is unaffected by the use of diploid rather than haploid data, an exception being the case of monogamous mating. Computer simulation shows close agreement of r2 with probability expectations, with minor discrepancies observed for unlinked loci. Using the probabilistic definition, r2 can be intuitively understood as a parameter shaping the effective sample size in genome-wide association studies.

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.

Microbial sequencing data are subject to sample-specific and species-specific distortions that complicate inference about true microbial abundances. Recent work established assumptions under which relative fold-differences in microbial abundances can be estimated from sequencing data, accompanied by an estimation and testing method called radEmu. While radEmu has excellent error rates, its computational cost can be prohibitive for large-scale analyses. Here we introduce fastEmu, a scalable approximation to radEmu's testing methodology. fastEmu uses a small, data-driven subset of taxa to approximate the computationally expensive test statistic. In a metagenomic analysis, we show a reduction in per-taxon runtime from 4-8 hours to 0.5-45 seconds, and simulations illustrate no adverse impact on estimation error or Type I error rate control. fastEmu enables testing differential abundance for thousands of microbial taxa while retaining clear parameter interpretation and error rate control.

Spatial transcriptomics enables high-throughput gene expression profiling while preserving spatial information, offering valuable insights into tissue architecture and cellular organization. Methods that achieve single-cell or subcellular spatial resolution typically rely on predefined gene panels, which limits genome-wide discovery. In contrast, sequencing-based spatial transcriptomics platforms provide broad transcriptome coverage but measure expression at lower spatial resolution, capturing mixtures of multiple cell types within each spatial location and thereby requiring accurate cell-type deconvolution. Many existing deconvolution methods do not jointly model molecular expression similarity and spatial context, fail to align predicted cell-type compositions with biological priors, and may produce spatially inconsistent results that do not reflect the smooth and structured organization of real tissues. We present ST-DISTAL, a dual-branch graph convolutional framework for cell-type deconvolution in spatial transcriptomics data. ST-DISTAL integrates conventional spatial graph convolution with spectral Chebyshev filtering through an adaptive attention mechanism, enabling the model to capture both local spatial dependencies and global structural patterns. The framework is trained using a composite objective that jointly enforces accurate cell-type proportion estimation, global distributional alignment with reference profiles, and spatial smoothness. Evaluations on three simulated benchmark datasets demonstrate that ST-DISTAL consistently outperforms state-of-the-art methods in both prediction accuracy and distributional alignment. Further validation on a real human embryonic heart dataset and a colorectal cancer dataset shows biologically plausible and spatially coherent cell-type organization.

High-dimensional mediators are routinely collapsed into low-dimensional scores before mediation analysis, but aggregation can change the mediated estimand rather than merely reduce noise. We characterize when a linear score preserves the full mediated effect. Exact preservation for every outcome model occurs precisely when the score contains the covariance-adjusted exposure-predictive mediator direction, the classical sufficient reduction for exposure under a Gaussian inverse-regression model. We introduce a distribution-free minimax adequacy loss equal to the fraction of exposure-relevant mediator information discarded by the score and exactly linked to overlap recovered by aggregation. A split-sample inverse-Wishart correction yields unbiased estimators of retained and full exposure information in moderate dimensions. Under spherically symmetric outcome-effect benchmarks, relative aggregation error has a Cauchy location-scale profile, and paired influence-function inference targets the realized gap for an observed outcome. In simulations, paired interval coverage ranged from 93.3% to 97.5%. In an outcome-free TCGA lung adenocarcinoma audit of 1,314 smoking-responsive genes, the pooled gene-mean loss was 0.811 (descriptive 95% gene-bootstrap interval 0.797 to 0.823), versus a dimension-only random-subspace benchmark of 0.913; retaining promoter and non-promoter means separately reduced pooled loss to 0.625. Independent half-sample screening was unstable, but top-100 and top-250 cross-screen pooled gene-mean losses ranged from 0.827 to 0.932. The framework separates structural adequacy, benchmark fragility, and outcome-specific discrepancy.