The three dimensional (3D) spatial organization of the genome is closely linked to biological functions and can be captured by Hi-C assays through interrogating genome-wide chromatin interactions. Methodologies for inferring 3D structures from Hi-C data summarized as a two-dimensional (2D) contact matrix can be broadly placed within the paradigms of optimization-based and sampling-based. Many optimization-based methods are capable of constructing whole genome 3D structures but do not account for spatial dependency in the 2D data matrix nor cell heterogeneity in bulk Hi-C data, which provide an average over millions of cells. Sampling-based methods, on the other hand, are probabilistic model-based and can account for not only dependency, heterogeneity, but also other features inherent in Hi-C data, such as over-dispersion and sparsity. However, whole-genome 3D structure recapitulation is too computationally expensive for sampling-based methods, while chromosome-by-chromosome strategies for sampling-based methods ignore important information on inter-chromosomal contacts. To address these issues, we propose the truncated Random effect EXpression-cut and paste (tREX-cap) method, which applies the tREX model within a divide and conquer strategy. The resulting method inherits the good data-feature-cognizant properties of tREX and, in the meantime, can efficiently infer the whole genome 3D structure. We demonstrate the performance of tREX-cap through an extensive simulation study and analyses of a Hi-C lymphoblastoid dataset and a Hi-C IMR90 dataset.
Accelerometer studies typically include repeated 24-h observations over several days and more recently have collected information on participants for weeks, months, or years. Meanwhile, there is growing evidence to suggest that components of physical activity beyond average or expected values are important for quantifying and understanding activity. Motivated by these emerging aspects, we introduce Multilevel Functional Quantile Principal Component Analysis (MFQPCA), a novel dimension-reduction technique that extends Functional Quantile Principal Component Analysis to hierarchical functional data. We examine accelerometry data collected for 4 to 9 months from congestive heart failure patients, where daily activity profiles are nested within individuals and exhibit substantial variability: participants often differ in the timing and intensity of physical activity behaviors across days and from each other, and capturing information beyond expected-value curves is necessary for a robust quantification of diurnal patterns. To this end, MFQPCA decomposes quantile-specific variability into between-participant and within-participant components, estimating shared patterns and corresponding scores at each level and quantile of interest. Our method robustly captures complex distributional features in hierarchical data, disentangles between- and within-participant sources of variability at different quantile levels, and facilitates the study of longitudinal changes in these distributions. We produce participant- and participant-day-level 10%, 30%, 50%, 70% and 90% quantile curves over 24 h, revealing that day-to-day variability dominates sedentary periods, whereas between-participant differences grow in vigorous activity. An efficient alternating minimization algorithm is implemented in the open-source R package FunQ, making MFQPCA readily applicable to a wide range of hierarchical functional data settings.
Advances in data acquisition technologies and computational resources have significantly improved the analysis of large datasets across various domains. These datasets often feature a high number of variables but a limited number of observations, commonly referred to as "high-dimensional data." Analyzing such data requires innovative statistical methods to further scientific insights. Mean vector testing plays a vital role in many scientific investigations. Although considerable research has been devoted to developing efficient methods for two-sided mean vector tests in high-dimensional data-such as those used to identify differentially expressed genes-one-sided high-dimensional mean vector tests have not been as thoroughly explored. These tests are particularly valuable for detecting significantly up-regulated or down-regulated gene sets within predefined groups. Addressing this need, we introduce a new method: the Sum Max-Component (SMC) test. We have explored the asymptotic behavior of the SMC test statistic as both the sample size and the data dimensions grow indefinitely. SMC test method has undergone extensive validation in finite sample scenarios, where it has proven to be effective, showing competitive rates for both type I error and power. To illustrate its practical applications, we have employed the SMC test in a gene set enrichment analysis on the proteomic data of ovarian cancer patients from the National Cancer Institute's CPTAC study, underscoring its potential value in ovarian cancer biology. Multivariate one-sided test, High-dimensional data, Gaussian-type tails, Exponential-type tails, enrichment analysis, high-grade serous ovarian cancer.
MicroRNAs play a central role in the regulation of gene expression and the modulation of diseases. Despite their importance, statistical methods for analyzing microRNAs have received less attention compared to messenger RNAs. Critically, messenger RNA sequencing methods are often applied to analyze microRNA sequencing data without considering the unique characteristics of microRNAs. This study examines the assumptions of messenger RNA-based methods and shows that they may incur high false discovery rates. We propose a negative binomial softmax regression (NBSR) model for microRNA sequencing data. Our approach has several advantages over existing methods. First, differential expression across experimental conditions is interpreted using the log relative abundance ratio (log-RAR). Second, it achieves greater statistical power and narrower confidence intervals by modeling the relationship between the biological coefficient of variation and relative abundance. Third, NBSR effectively handles highly variable and sparsely expressed microRNAs, resulting in improved sensitivity for detecting differential expression. Additionally, we show that debiasing the log-RAR enables accurate inference of fold changes in absolute abundance, particularly when only a small subset of microRNAs differ in expression between conditions. We demonstrate the efficacy of our approach using both real and simulated data.
A live clinical question is: which patients benefit from intensive care unit (ICU) transfer? Motivated by this question we address the problem of estimating conditional average treatment effects (CATE) in the presence of unmeasured confounding, by leveraging instrumental variables (IVs). CATE-learners (eg R-learner) developed for settings without unmeasured confounding can readily incorporate IVs by substituting the observed exposure with first-stage predictions from a regression of treatment on IVs and covariates. Such predictions may be obtained via flexible data-adaptive methods (eg statistical or machine learning procedures) to alleviate concerns about model misspecification. However, the large regularization bias typical of data-adaptive predictions may propagate into the CATE estimates, resulting in poor accuracy. Neyman-orthogonal learners have therefore been developed, which prevent this by "insulating" the resulting CATE estimates against bias and estimation errors in these predictions. However, synthetic data simulations reveal that previously proposed Neyman-orthogonal learners for IV regression perform poorly. We remedy this problem using infinite-dimensional targeted learning, which strategically tailors first-stage predictions to perform well in their ultimate task: delivering accurate, precise CATE estimates. The resulting targeted Neyman-orthogonal learner is easy to construct based on arbitrary, off-the-shelf learners. It can handle continuous or discrete exposures, and arbitrary types and numbers of IVs and covariates. Simulation studies and a re-analysis of the benefits of ICU transfer show substantial enhancements in performance, underscoring the importance of the proposed IV-learner.
Regression analysis of correlated data, where multiple correlated responses are recorded on the same unit, is ubiquitous in many scientific areas. With the advent of new technologies, in particular high-throughput omics profiling assays, such correlated data increasingly consist of a large number of variables compared with the available sample size. Motivated by recent longitudinal proteomics studies of COVID-19, we propose a novel inference procedure for linear functionals of high-dimensional regression coefficients in generalized estimating equations, which are widely used to analyze correlated data. Our estimator for this more general inferential target, obtained via constructing projected estimating equations, is shown to be asymptotically normally distributed under mild regularity conditions. We also introduce a data-driven cross-validation procedure to select the tuning parameter for estimating the projection direction, which is not addressed in the existing procedures. We illustrate the utility of the proposed procedure in providing confidence intervals for associations of individual proteins and severe COVID risk scores obtained based on high-dimensional proteomics data, and demonstrate its robust finite-sample performance, especially in estimation bias and confidence interval coverage, via extensive simulations.
In situations where high-quality external data are available or when it is challenging to recruit participants to the control arm of a randomized and controlled clinical trial (eg rare or pediatric diseases), it is desirable to borrow information from external data to augment the control arm. However, a main challenge in borrowing information from external data is to accommodate potential heterogeneous subpopulations across the external and trial data. We apply a Bayesian nonparametric model called the Shared Atoms Model (SAM) to identify overlapping and unique subpopulations across datasets, with which we restrict the information borrowing to the common subpopulations. This forms a hybrid control (HC) that leads to more precise estimation of treatment effects. The degree of information borrowing is confined by the sample size and degree of similarity in outcomes. Simulation studies demonstrate the robustness of the new method, and an application to an Atopic Dermatitis dataset shows improved treatment effect estimation.
New SARS-CoV-2 variants arise frequently with different viral properties that can impact the effectiveness of the vaccines. Updating estimates of vaccine effectiveness (VE) in public health surveillance can be limited by the necessity of conducting a distinct study that entails analysis of prospective cohort data or using a test-negative design. We introduce a method for dynamically updating estimates of VE using data that accumulate in real time. Our method uses dynamic case-control sampling to estimate VE against a newly emerging variant relative to a previous variant. Dynamic case-control sampling is a technique that continuously updates VE estimates by comparing individuals infected with a newly emerging variant (defined as "cases") to those infected with a previously circulating variant (defined as "controls"). We use this estimate in combination with information about VE from the previous variant (these estimates are typically available from larger, traditional studies) to infer VE against the emerging variant. We demonstrate the utility of this method on the BA.1 and BA.2 sub-lineages of the Omicron variant. The method produces estimates of VE comparable to those produced using traditional methods, although with increased SE. The increase in error, however, is reasonable given a much smaller sample size than other studies, and error ranges of the estimates could be significantly improved by sequencing a larger proportion of identified cases. Our method, which assumes only a fraction of the new cases are being sequenced, can be applied by health departments using routinely collected data to produce timely, rigorous VE estimates to rapidly identify potential changes in VE.
High-dimensional longitudinal data are increasingly available in biomedical research, especially from omics platforms, but pose substantial challenges for joint modeling with survival outcomes. These challenges include modeling complex temporal dynamics, accommodating cross-feature dependencies, and maintaining computational feasibility. We propose a novel joint modeling framework that addresses these issues using supervised low-rank functional tensor decomposition to capture latent structure in multivariate longitudinal data and proportional hazards modeling for time-to-event outcomes. The longitudinal process is represented as a multivariate functional tensor, with a low-rank approximation that incorporates supervision from baseline covariates. Estimation is performed using a likelihood-based Monte Carlo Expectation-Maximization algorithm, enabling coherent inference and individualized prediction. Our method produces dynamic predictions of both longitudinal feature trajectories and survival probabilities. Simulation studies demonstrate substantial improvements in estimation accuracy and predictive performance over a standard two-stage approach, particularly under high censoring and limited sample sizes. In application to the Alzheimer's Disease Neuroimaging Initiative lipidomics data, the proposed model explains over 99% of variation with four components, and identifies significant subject-level latent predictors of dementia onset. This framework provides a scalable and interpretable strategy for integrating high-dimensional longitudinal biomarkers into joint models for disease progression and risk stratification.
The infant microbiome undergoes rapid changes in composition over time and is associated with long-term risks of conditions such as immune strength, allergy, asthma, and other health outcomes. Modeling the associations between exposures or treatments and microbial composition over time is essential for understanding the factors that drive these changes. Estimating these temporal dynamics has several challenges including repeated measures, overdispersion, compositionality, high-dimensional parameter spaces, and zero-inflation. Many longitudinal regression models used in human microbiome research assume constant effects over time that cannot capture time-varying or functional effects of exposures, ignore the compositional structure of the data by modeling each taxon separately, and are not equipped to handle potential zero-inflation. Dirichlet-multinomial (DM) regression models inherently accommodate overdispersion and the compositional structure of the data and have been extended to account for excess zeros. However, existing DM-based regression models are unable to additionally handle repeated measures designs. To fill this gap, we propose a functional concurrent zero-inflated Dirichlet-multinomial regression model which is designed to model time-varying relations between observed covariates and microbial taxa while accounting for zero-inflation, compositionality, and repeated measures. Through simulation, we demonstrate that the model can accurately estimate the underlying functional relations and scale to large compositional spaces. We apply our model to investigate time-varying associations between infant microbiome composition and observed covariates during the 11-wk postnatal period. We found that $ \boldsymbol{\alpha} $-diversity (ie the diversity of the microbiome within an individual) is positively associated with a higher gestational age and percentage of breast milk in the diet. We provide an accompanying R package and shiny app to implement the method and generate plots.
Egocentric-Network Randomized Trials (ENRTs) are increasingly used to estimate causal effects under interference when measuring complete sociocentric network data is infeasible. ENRTs rely on egocentric network sampling, where a set of egos is first sampled, and each ego recruits a subset of its neighbors as alters. Treatments are then randomized across egos. While the observed ego-networks are disjoint by design, the underlying population network may contain edges connecting them, leading to contamination. Under a design-based framework, we show that the Horvitz-Thompson estimators of direct and indirect effects are biased whenever contamination is present. To address this, we derive bias-corrected estimators and propose a novel sensitivity analysis framework based on sensitivity parameters representing the probability or expected number of missing edges. This framework is implemented via both grid sensitivity analysis and probabilistic bias analysis, providing researchers with a flexible tool to assess the robustness of the causal estimators to contamination. We apply our methodology to the HIV Prevention Trials Network 037 study, finding that ignoring contamination may lead to underestimation of indirect effects and overestimation of direct effects.
Mortality risk estimated from studies that ascertain date of death through linkage to vital statistics registries may be subject to outcome measurement error. As a result, some deaths among study participants may not be captured, some study participants who are alive may be falsely categorized as deceased, and some deaths may be recorded at incorrect times, leading to bias in estimates of mortality risk and survival. Here, we illustrate an extension of the Rogan-Gladen estimator to account for outcome measurement error in risk and survival functions in settings with right censoring. As a motivating application, we consider and account for outcome measurement error that could be induced by incomplete and/or incorrect linkage to death registries when estimating mortality risk among people entering care for HIV in the University of North Carolina Center for AIDS Research HIV Clinical Cohort between 2001 and 2022. A series of simulation studies demonstrates that the approach performed well even when participants selected into the validation study were at higher mortality risk than the main study. The proposed approach may be parameterized using internal or external validation data or used as a form of quantitative bias analysis.
A number of domains in biomedical research use data with a large number of predictors all representing the same type of measurement. Often, an important summary is the within-person distribution of these predictors. Here we focus on settings where the mean relationship between outcome and predictors is fully captured by this distribution and, more generally, on problems where the goal is to learn a mapping that is invariant under permutations of the input vector. We compare unstructured neural networks, which do not explicitly incorporate the permutation invariance property, versus networks that we call ordered predictors neural networks. We show in simulations that the unstructured deep learning approach can yield higher prediction errors, compared to the approach that explicitly leverages the invariance to simplify the learning task. Additionally, in the context of neural Bayes estimation, in which neural networks are used to construct point estimators, we show that ordered predictors neural networks can yield substantially more precise estimators. We therefore recommend that, when permutation invariance is known or suspected to hold, investigators use a learning or statistical modeling approach that can leverage the invariance, rather than an unstructured deep learning approach.
Large-scale neuroimaging studies often collect data from multiple scanners across different sites, where variations in scanners, scanning procedures, and other conditions across sites can introduce artificial site effects. These effects may bias brain connectivity measures, such as functional connectivity, which quantify functional network organization derived from functional magnetic resonance imaging. How to leverage high-dimensional network structures to effectively mitigate site effects has yet to be addressed. In this paper, we propose Sparse LAtent Covariate-driven Connectome (SLACC) factorization, a multivariate method that explicitly parameterizes covariate effects in latent subject scores corresponding to sparse rank-1 latent patterns derived from brain connectivity. The proposed method identifies localized site-driven variability within and across brain networks, enabling targeted correction. We develop a penalized Expectation-Maximization algorithm for parameter estimation, incorporating the Bayesian Information Criterion to guide optimization. Extensive simulations validate SLACC's robustness in recovering the true parameters and underlying connectivity patterns. Applied to the Autism Brain Imaging Data Exchange dataset, SLACC demonstrates its ability to reduce site effects.
Existing approaches to modelling antibody concentration data are mostly based on finite mixture models that rely on the assumption that individuals can be divided into 2 distinct groups: seronegative and seropositive. Here, we challenge this dichotomous modelling assumption and propose a latent variable modelling framework in which the immune status of each individual is represented along a continuum of latent seroreactivity, ranging from minimal to strong immune activation. This formulation provides greater flexibility in capturing age-related changes in antibody distributions while preserving the full information content of quantitative measurements. We show that the proposed class of models can accommodate a large variety of model formulations, both mechanistic and regression-based, and also includes finite mixture models as a special case. We also propose a computationally efficient $ L_{2} $-based estimator as an alternative to maximum likelihood estimation, which substantially reduces computational cost, and we establish its consistency. Through a case study on malaria serology, we demonstrate how the flexibility of the novel framework enables joint analyses across all ages while accounting for changes in transmission patterns. We conclude by outlining extensions of the proposed modelling framework and its relevance to other omics applications.
Understanding the pathways through which diet affects human metabolism is a central task in nutritional epidemiology. This article proposes novel methodology to identify food items associated with blood metabolites in 2 cohorts of healthcare professionals. We analyze 244 metabolites characterized by statistical complexities that include skewness, left-censoring, and structural missingness. Though existing methods can address such factors in low-dimensional settings, they cannot exploit the nutritional or statistical relationships among the 30 considered food intake variables, and they are unsuitable for performing high-dimensional inference. To address these challenges, we develop a novel Bayesian variable selection framework for metabolite response variables based on a skew-normal censored mixture model, while exploiting substantive information on the considered food items via a Markov random field prior. Applying this methodology to the cohort data identifies multiple metabolite-diet associations that are consistent with previous research as well as several potentially novel associations that were not detected using standard methods. The proposed approach is implemented in the R package multimetab, facilitating its use in high-dimensional metabolomic analyses.
Serum prostate-specific antigen (PSA) is widely used for prostate cancer screening. While the genetics of PSA levels has been studied to enhance screening accuracy, the genetic basis of PSA velocity, the rate of PSA change over time, remains unclear. The Prostate, Lung, Colorectal, and Ovarian (PLCO) Cancer Screening Trial, a large, randomized study with longitudinal PSA data (15,260 cancer-free males, averaging 5.34 samples per subject) and genome-wide genotype data, provides a unique opportunity to estimate PSA velocity heritability. We developed a mixed model to jointly estimate heritability of PSA levels at age 54 and PSA velocity. To accommodate the large dataset, we implemented two efficient computational approaches: a partitioning and meta-analysis strategy using average information restricted maximum likelihood (AI-REML), and a fast restricted Haseman-Elston (REHE) regression method. Simulations showed that both methods yield unbiased estimates of both heritability metrics, with AI-REML providing smaller variability in the estimation of velocity heritability than REHE. Applying AI-REML to PLCO data, we estimated heritability at 0.32 (s.e. = 0.07) for baseline PSA and 0.45 (s.e. = 0.18) for PSA velocity. These findings reveal a substantial genetic contribution to PSA velocity, supporting future genome-wide studies to identify variants affecting PSA dynamics and improve PSA-based screening.
Advances in cellular imaging technologies, especially those based on fluorescence in situ hybridization (FISH) now allow detailed visualization of the spatial organization of human or bacterial cells. Quantifying this spatial organization is crucial for understanding the function of multicellular tissues or biofilms, with implications for human health and disease. To address the need for better methods to achieve such quantification, we propose a flexible multivariate point process model that characterizes and estimates complex spatial interactions among multiple cell types. The proposed Bayesian framework is appealing due to its unified estimation process and the ability to directly quantify uncertainty in key estimates of interest, such as those of inter-type correlation and the proportion of variance due to inter-type relationships. To ensure stable and interpretable estimation, we consider shrinkage priors for coefficients associated with latent processes. Model selection and comparison are conducted by using a deviance information criterion designed for models with latent variables, effectively balancing the risk of overfitting with that of oversimplifying key quantities. Furthermore, we develop a hierarchical modeling approach to integrate multiple image-specific estimates from a given subject, allowing inference at both the global and subject-specific levels. We apply the proposed method to microbial biofilm image data from the human tongue dorsum and find that specific taxon pairs, such as Streptococcus mitis-Streptococcus salivarius and Streptococcus mitis-Veillonella, exhibit strong positive spatial correlations, while others, such as Actinomyces-Rothia, show slight negative correlations. For most of the taxa, a substantial portion of spatial variance can be attributed to inter-taxon relationships.
To address the challenges in modeling time-to-event outcomes in small-sample settings, we propose a novel transfer learning approach, termed CoxTL, based on the widely used Cox proportional hazards model, accounting for potential covariate and concept shifts between source and target datasets. CoxTL utilizes a combination of density ratio weighting and importance weighting techniques to address multi-level data heterogeneity, including covariate and coefficient shifts between source and target datasets. Additionally, it accounts for potential model misspecification, ensuring robustness across a wide range of settings. We assess the performance of CoxTL through extensive simulation studies, considering data under various types of distributional shifts. In addition, we apply CoxTL to predict End-Stage Renal Disease (ESRD) in the Hispanic population using electronic health record-derived features from the All of Us Research Program. Data from non-Hispanic White and non-Hispanic Black populations are leveraged as source cohorts. Model performance is evaluated using the C-index and Integrated Brier Score (IBS). In simulation studies, CoxTL demonstrates higher predictive accuracy, particularly in scenarios involving multi-level heterogeneity between target and source datasets. In other scenarios, CoxTL performs comparably to alternative methods specifically designed to address only a single type of distributional shift. For predicting the 2-year risk of ESRD in the Hispanic population, CoxTL achieves an increase in the C-index of up to 6.76% compared to the model trained exclusively on target data. Furthermore, it demonstrates up to a 17.94% increase in the C-index compared to the state-of-the-art transfer learning method based on Cox model. The proposed method effectively utilizes source data to enhance time-to-event predictions in target populations with limited samples. Its ability to handle various sources and levels of data heterogeneity ensures robustness, making it particularly well-suited for real-world applications involving target populations with small sample sizes, where traditional Cox models often struggle.
Single-cell RNA sequencing allows the quantification of gene expression at the individual cell level, enabling the study of cellular heterogeneity and gene expression dynamics. Dimensionality reduction is a common preprocessing step critical for the visualization, clustering, and phenotypic characterization of samples. This step, often performed using principal component analysis or closely related methods, is challenging because of the size and complexity of the data. In this work, we present a generalized matrix factorization model assuming a general exponential dispersion family distribution and we show that many of the proposed approaches in the single-cell dimensionality reduction literature can be seen as special cases of this model. Furthermore, we propose a scalable adaptive stochastic gradient descent algorithm that allows us to estimate the model efficiently, enabling the analysis of millions of cells. We benchmark the proposed algorithm through extensive numerical experiments against state-of-the-art methods and showcase its use in real-world biological applications. The proposed method systematically outperforms existing methods of both generalized and non-negative matrix factorization, demonstrating faster execution times and parsimonious memory usage, while maintaining, or even enhancing, matrix reconstruction fidelity and accuracy in biological signal extraction. On real data, we show that our method scales seamlessly to millions of cells, enabling dimensionality reduction in large single-cell datasets. Finally, all the methods discussed here are implemented in an efficient open-source R package, sgdGMF, available on CRAN.
