Applications of single-cell RNA sequencing in various biomedical research areas have been blooming. This new technology provides unprecedented opportunities to study disease heterogeneity at the cellular level. However, unique characteristics of scRNA-seq data, including large dimensionality, high dropout rates, and possibly batch effects, bring great difficulty into the analysis of such data. Not appropriately addressing these issues obstructs true scientific discovery. Herein we propose a unified Regularized Zero-inflated Mixture Model framework, designed for scRNA-seq data (RZiMM-scRNA), to simultaneously detect cell subgroups and identify gene differential expression based on a developed importance score, accounting for both dropouts and batch effects. We conduct extensive simulation studies in which we evaluate the performance of RZiMM-scRNA and compare it with several popular methods, including Seurat, SC3, K-means, and hierarchical clustering. Simulation results show that RZiMM-scRNA demonstrates superior clustering performance and enhanced biomarker detection accuracy, compared to alternative methods, especially when cell subgroups are less distinct, verifying the robustness of our method.

Our empirical investigations focus on two brain tumor studies dealing with astrocytoma of various grades, including the most malignant of all brain tumors, glioblastoma multiforme (GBM). Our goal is to delineate cell heterogeneity and identify driving biomarkers associated with these tumors. Notably, RZiMM-scNRA successfully identifies a small group of oligodendrocyte cells, which has drawn much attention in biomedical literature on brain cancers. In addition, our method discovers several new biomarkers which are not discussed in the original studies, including PLP1, BCAN, and PTPRZ1—all associated with the development and malignant growth of glioma—as well as CAMK2B, which is downregulated in glioma and GBM and implicated in neurodevelopment, brain function, learning and memory processes.

The transition from nonrenewable to renewable energies represents a global societal challenge, and developing a sustainable energy portfolio is an especially daunting task for developing countries where little to no information is available regarding the abundance of renewable resources such as wind. Weather model simulations are key to obtain such information when observational data are scarce and sparse over a country as large and geographically diverse as Saudi Arabia. However, output from such models is uncertain, as it depends on inputs such as the parametrization of the physical processes and the spatial resolution of the simulated domain. In such situations a sensitivity analysis must be performed, and the input may have a spatially heterogeneous influence of wind. In this work we propose a latent Gaussian functional analysis of variance (ANOVA) model that relies on a nonstationary Gaussian Markov random field approximation of a continuous latent process. The proposed approach is able to capture the local sensitivity of Gaussian and non-Gaussian wind characteristics such as speed and threshold exceedances over a large simulation domain, and a continuous underlying process also allows us to assess the effect of different spatial resolutions. Our results indicate that: (1) the nonlocal planetary boundary layer scheme and high spatial resolution are both instrumental in capturing wind speed and energy (especially over complex mountainous terrain), and (2) the impact of planetary boundary layer scheme and resolution on Saudi Arabia’s planned wind farms is small (at most 1.4%). Thus, our results lend support for the construction of these wind farms in the next decade.

An experimental comparison of two or more optimization algorithms requires the same computational resources to be assigned to each algorithm. When a maximum runtime is set as the stopping criterion, all algorithms need to be executed in the same machine if they are to use the same resources. Unfortunately, the implementation code of the algorithms is not always available, which means that running the algorithms to be compared in the same machine is not always possible. And even if they are available, some optimization algorithms might be costly to run, such as training large neural-networks in the cloud.

In this paper we consider the following problem: how do we compare the performance of a new optimization algorithm B with a known algorithm A in the literature if we only have the results (the objective values) and the runtime in each instance of algorithm A? Particularly, we present a methodology that enables a statistical analysis of the performance of algorithms executed in different machines. The proposed methodology has two parts. First, we propose a model that, given the runtime of an algorithm in a machine, estimates the runtime of the same algorithm in another machine. This model can be adjusted so that the probability of estimating a runtime longer than what it should be is arbitrarily low. Second, we introduce an adaptation of the one-sided sign test that uses a modified p-value and takes into account that probability. Such adaptation avoids increasing the probability of type I error associated with executing algorithms A and B in different machines.

Real-world evidence synthesis through integration of data from distributed research networks has gained increasing attention in recent years. Due to privacy concerns and restrictions of sharing patient-level data, distributed algorithms that do not require sharing patient level information are in great need for facilitating multisite collaborations. On the other hand, data collected at multiple sites often come from diverse populations, and there exists a substantial amount of heterogeneity across sites in patient characteristics. Most of the existing distributed algorithms have ignored such between-site heterogeneity. In this paper we aim to fill this methodological gap by proposing a general distributed algorithm. We develop our distributed algorithm based on a general semiparametric model, namely, the proportional likelihood ratio model (Biometrika 99 (2012) 211–222), which is a semiparametric extension of generalized linear model. We devise the proportional likelihood ratio model with site-specific baseline function, to account for between-site heterogeneity, and shared regression parameters to borrow information across sites. Under this flexible formulation, our distributed algorithm is designed to be privacy-preserving and communication-efficient (i.e., only one round of communication across sites is needed). We validate our method via simulation studies and demonstrate the utility of our method via a multisite study of pediatric avoidable hospitalization based on electronic health record data from a total of 354,672 patients across 26 different clinical sites within the Children’s Hospital of Philadelphia health system.

In multiple instance learning (MIL), the response variable is predicted by features (or covariates) of one or more instances, which are collectively denoted as a bag. Learning the relationship between bags and instances is challenging because of the unknown and possibly complicated data generating mechanism regarding how instances contribute to the bag label. MIL has been applied to solve a variety of real-world problems, which have been mostly focused on supervised tasks, such as molecule activity prediction, protein binding affinities prediction, object detection, and computer-aided diagnosis. However, to date, the majority of the off-the-shelf MIL methods are developed in the computer science domain, and they focus on improving the prediction performance while spending little effort on explainability of the algorithm. In this article a Bayesian multiple instance learning model, based on probit regression (MICProB), is proposed, which contributes a significant portion to the suite of statistical methodologies for MIL. MICProB is composed of two nested probit regression models, where the inner model is estimated for predicting primary instances, which are considered as the “important” ones that determine the bag label, and the outer model is for predicting bag-level responses based on the primary instances estimated by the inner model. The posterior distribution of MICProB can be conveniently approximated using a Gibbs sampler, and the prediction for new bags can be performed in a fully integrated Bayesian way. We evaluate the performance of MICProB against 15 benchmark methods and demonstrate its competitiveness in simulation and real-data examples. In addition to its capability of identifying primary instances, as compared to existing optimization-based approaches, MICProB also enjoys great advantages in providing a transparent model structure, straightforward statistical inference of quantities related to model parameters, and favorable interpretability of covariate effects on the bag-level response.

Diffusion-weighted magnetic resonance imaging (D-MRI) is an in vivo and noninvasive imaging technology for characterizing tissue microstructure in biological samples. A major application of D-MRI is for white matter fiber tract reconstruction in brains. It begins by estimating the water molecule movements (serving as proxies for fiber directions) in the brain voxels and then combines the results to form fiber tracts. The voxel-level fiber direction information can be modeled by a fiber orientation distribution (FOD) function, and in this paper, we propose a computationally scalable FOD estimator, the blockwise James–Stein (BJS) estimator. We then apply BJS to the D-MRI data from the Human Connectome Project (HCP) to study brain lateralization, an important topic in neuroscience. Specifically, we focus on the association between lateralization of the superior longitudinal fasciculus (SLF)—a major association tract and handedness. For each subject from the HCP data, we extract voxel-level directional information by BJS and then reconstruct the SLF in each brain hemisphere through a tractography algorithm. Finally, we derive a lateralization score that quantifies hemispheric asymmetry of the reconstructed SLF. We then relate this lateralization score to gender and handedness through an ANOVA model, where significant handedness effects are found. The results indicate that the SLF lateralization is likely to be different in right-handed and left-handed individuals. Codes and example scripts for both synthetic experiments and HCP data application can be found at https://github.com/vic-dragon/BJS.

Systematic measurement error in self-reported data creates important challenges in association studies between dietary intakes and chronic disease risks, especially when multiple dietary components are studied jointly. The joint regression calibration method has been developed for measurement error correction when objectively measured biomarkers are available for all dietary components of interest. Unfortunately, objectively measured biomarkers are only available for very few dietary components, which limits the application of the joint regression calibration method. Recently, for single dietary components, controlled feeding studies have been performed to develop new biomarkers for many more dietary components. However, it is unclear whether the biomarkers separately developed for single dietary components are valid for joint calibration. In this paper we show that biomarkers developed for single dietary components cannot be used for joint regression calibration. We propose new methods to utilize controlled feeding studies to develop valid biomarkers for joint regression calibration to estimate the association between multiple dietary components simultaneously with the disease of interest. Asymptotic distribution theory for the proposed estimators is derived. Extensive simulations are performed to study the finite sample performance of the proposed estimators. We apply our methods to examine the joint effects of sodium and potassium intakes on cardiovascular disease incidence using the Women’s Health Initiative cohort data. We identify positive associations between sodium intake and cardiovascular diseases as well as negative associations between potassium intake and cardiovascular disease.

Online health communities (OHCs) provide free, open, and well-resourced platforms for patients, family members, and others to discuss illnesses, express feelings, and connect with others. Linguistic analysis of OHC posts can assist in better understanding disease conditions as well as monitoring the emotional and mental status of patients and those who are closely related. Many existing OHC linguistic analyses are limited by focusing on individual words. There are a handful of cooccurrence network analyses, which have multiple methodological limitations. In this article we analyze posts that are publicly available at the LUNGevity Foundation’s Lung Cancer Support Community (LCSC). The analyzed data contains 21,028 posts published between April 2018 and February 2022. For word cooccurrence network analysis, we develop a two-part latent space model, which advances from the existing ones by accommodating network weights. Further, we consider the scenario where there are change points in time, networks remain the same between two change points but differ on the two sides of a change point, and the number and locations of change points are unknown. A penalized fusion approach is developed to data-dependently determine change points and estimate networks. In data analysis multiple change points are identified, which reflect significant changes in lung cancer patients’ and their close affiliates’ emotional/mental status and mostly align with the changes in COVID-19. The obtained network structures and other findings are also sensible.

Ensemble weather forecasts based on multiple runs of numerical weather prediction models typically show systematic errors and require postprocessing to obtain reliable forecasts. Accurately modeling multivariate dependencies is crucial in many practical applications, and various approaches to multivariate postprocessing have been proposed where ensemble predictions are first postprocessed separately in each margin and multivariate dependencies are then restored via copulas. These two-step methods share common key limitations, in particular, the difficulty to include additional predictors in modeling the dependencies. We propose a novel multivariate postprocessing method based on generative machine learning to address these challenges. In this new class of nonparametric data-driven distributional regression models, samples from the multivariate forecast distribution are directly obtained as output of a generative neural network. The generative model is trained by optimizing a proper scoring rule, which measures the discrepancy between the generated and observed data, conditional on exogenous input variables. Our method does not require parametric assumptions on univariate distributions or multivariate dependencies and allows for incorporating arbitrary predictors. In two case studies on multivariate temperature and wind speed forecasting at weather stations over Germany, our generative model shows significant improvements over state-of-the-art methods and particularly improves the representation of spatial dependencies.

Correlated data are ubiquitous in today’s data-driven society. While regression models for analyzing means and variances of responses of interest are relatively well developed, the development of these models for analyzing the correlations is largely confined to longitudinal data, a special form of sequentially correlated data. This paper proposes a new method for the analysis of correlations to fully exploit the use of covariates for general correlated data. In a renewed analysis of the classroom data, a highly unbalanced multilevel clustered data with within-class and within-school correlations, our method reveals informative insights on these structures not previously known. In another analysis of the malaria immune response data in Benin, a longitudinal study with time-dependent covariates where the exact times of the observations are not available, our approach again provides promising new results. At the heart of our approach is a new generalized z-transformation that converts correlation matrices, constrained to be positive definite, to vectors with unrestricted support and is order-invariant. These two properties enable us to develop regression analysis incorporating covariates for the modelling of correlations via the use of maximum likelihood.

The prevalence of security threats to organizational defense demands models that support real-world policymaking. Security games are a potent tool in this regard; however, although canonical models effectively allocate limited resources, they generally do not consider adaptive, boundedly rational adversaries. Empirical findings suggest this characterization describes real-world human behavior, so the development of decision-support frameworks against such adversaries is a critical need. We examine a family of policies applicable to repeated games in which a boundedly rational adversary is modeled using a behavioral-economic theory of learning, that is, experience-weighted attraction learning. These policies take into account realistic uncertainty about the competition by adopting the perspective of adversarial risk analysis. Using Bayesian reasoning, these repeated games are decomposed into multiarm bandit problems. A collection of cost-function approximation policies are given to solve these problems. The efficacy of our approach is shown via extensive computational testing on a defense-related case study.

The COVID-19 pandemic has been a worldwide health crisis for the past three years, casting unprecedented challenges for policymakers in different countries and regions. While one country or region can only implement one social mobility restriction policy at a given time, it is of great interest for policy makers to decide whether to elevate or deelevate the restriction policy from time to time. This article proposes a novel nonnegative tensor completion method to predict the potential counterfactual outcomes of multifaceted social mobility restriction policies over time. The proposed method builds upon a low-rank tensor decomposition of the pandemic data, which also explicitly characterizes the ordinal nature of the mobility restriction strength and the smooth trend of the pandemic evolution over time. Its application to the COVID-19 pandemic data reveals some interesting facts regarding the impact of social mobility restriction policy on the spread of the virus. The effectiveness of the proposed method is also supported by its asymptotic estimation consistency and extensive numerical experiments on the synthetic datasets.

Dynamic risk prediction that incorporates longitudinal measurements of biomarkers is useful in identifying high-risk patients for better clinical management. Our work is motivated by the prediction of cervical precancers. Currently, Pap cytology is used to identify HPV-positive (HPV+) women at high-risk of cervical precancer, but cytology lacks accuracy and reproducibility. Molecular markers, like HPV DNA methylation, that are closely linked to the carcinogenic process show promise of improved risk stratification. We are interested in developing a dynamic risk model that uses all longitudinal biomarker information to improve precancer risk estimation. We propose a joint model to link both the continuous methylation biomarker and a binary cytology biomarker to the time to precancer outcome using shared random effects. The model uses a discretization of the time scale to allow for closed-form likelihood expressions, thereby avoiding potential high dimensional integration of the random effects. The method handles an interval-censored time-to-event outcome, due to intermittent clinical visits, incorporates sampling weights to deal with stratified sampling data and can provide immediate and five-year risk estimates that may inform clinical decision-making. Applying the method to longitudinally measured HPV methylation data improves risk stratification for triage of HPV+ women.

There is increasing appetite for analysing populations of network data due to the fast-growing body of applications demanding such methods. While methods exist to provide readily interpretable summaries of heterogeneous network populations, these are often descriptive or ad hoc, lacking any formal justification. In contrast, principled analysis methods often provide results difficult to relate back to the applied problem of interest. Motivated by two complementary applied examples, we develop a Bayesian framework to appropriately model complex heterogeneous network populations, while also allowing analysts to gain insights from the data and make inferences most relevant to their needs. The first application involves a study in computer science measuring human movements across a university. The second analyses data from neuroscience investigating relationships between different regions of the brain. While both applications entail analysis of a heterogeneous population of networks, network sizes vary considerably. We focus on the problem of clustering the elements of a network population, where each cluster is characterised by a network representative. We take advantage of the Bayesian machinery to simultaneously infer the cluster membership, the representatives, and the community structure of the representatives, thus allowing intuitive inferences to be made. The implementation of our method on the human movement study reveals interesting movement patterns of individuals in clusters, readily characterised by their network representative. For the brain networks application, our model reveals a cluster of individuals with different network properties of particular interest in neuroscience. The performance of our method is additionally validated in extensive simulation studies.

The natural cycles of the surface-to-atmosphere fluxes of carbon dioxide (CO₂) and other important greenhouse gases are changing in response to human influences. These changes need to be quantified to understand climate change and its impacts, but this is difficult to do because natural fluxes occur over large spatial and temporal scales and cannot be directly observed. Flux inversion is a technique that estimates the spatiotemporal distribution of a gas’ fluxes using observations of the gas’ mole fraction and a chemical transport model. To infer trends in fluxes and identify phase shifts and amplitude changes in flux seasonal cycles, we construct a flux-inversion system that uses a novel spatially-varying time-series decomposition of the fluxes. We incorporate this decomposition into the Wollongong Methodology for Bayesian Assimilation of Trace-gases (WOMBAT, Zammit-Mangion et al., Geosci. Model Dev., 15, 2022), a Bayesian hierarchical flux-inversion framework that yields posterior distributions for all unknowns in the underlying model. We also extend WOMBAT to accommodate physical constraints on the fluxes and to take direct in situ and flask measurements of trace-gas mole fractions as observations. We apply the new method, which we call WOMBAT v2.0, to a mix of satellite observations of CO₂ mole fraction from the Orbiting Carbon Observatory-2 (OCO-2) satellite and direct measurements of CO₂ mole fraction from a variety of sources. We estimate the changes in the natural cycles of CO₂ fluxes that occurred from January 2015 to December 2020, and compare our posterior estimates to those from an alternative method based on a bottom-up understanding of the physical processes involved. We find substantial trends in the fluxes, including that tropical ecosystems trended from being a net source to a net sink of CO₂ over the study period. We also find that the amplitude of the global seasonal cycle of ecosystem CO₂ fluxes increased over the study period by 0.11 PgC/month (an increase of 8%) and that the seasonal cycle of ecosystem CO₂ fluxes in the northern temperate and northern boreal regions shifted earlier in the year by 0.4–0.7 and 0.4–0.9 days, respectively (2.5th to 97.5th posterior percentiles), consistent with expectations for the carbon cycle under a warming climate.

We propose a novel analysis of power (ANOPOW) model for analyzing replicated nonstationary time series commonly encountered in experimental studies. Based on a locally stationary ANOPOW Cramér spectral representation, the proposed model can be used to compare the second-order time-varying frequency patterns among different groups of time series and to estimate group effects as functions of both time and frequency. Formulated in a Bayesian framework, independent two-dimensional second-order random walk (RW2D) priors are assumed on each of the time-varying functional effects for flexible and adaptive smoothing. A piecewise stationary approximation of the nonstationary time series is used to obtain localized estimates of time-varying spectra. Posterior distributions of the time-varying functional group effects are then obtained via integrated nested Laplace approximations (INLA) at a low computational cost. The large-sample distribution of local periodograms can be appropriately utilized to improve estimation accuracy since INLA allows modeling of data with various types of distributions. The usefulness of the proposed model is illustrated through two real-data applications: analyses of seismic signals and pupil diameter time series in children with attention deficit hyperactivity disorder. Simulation studies, Supplementary Material (Li, Yue and Bruce (2024a)), and $\mathtt{R}$ code (Li, Yue and Bruce (2024b)) for this article are also available.

Assessing heterogeneity in the effects of treatments has become increasingly popular in the field of causal inference and carries important implications for clinical decision-making. While extensive literature exists for studying treatment effect heterogeneity when outcomes are fully observed, there has been limited development in tools for estimating heterogeneous causal effects when patient-centered outcomes are truncated by a terminal event, such as death. Due to mortality occurring during study follow-up, the outcomes of interest are unobservable, undefined, or not fully observed for many participants in which case principal stratification is an appealing framework to draw valid causal conclusions. Motivated by the Acute Respiratory Distress Syndrome Network (ARDSNetwork) ARDS respiratory management (ARMA) trial, we developed a flexible Bayesian machine learning approach to estimate the average causal effect and heterogeneous causal effects among the always-survivors stratum when clinical outcomes are subject to truncation. We adopted Bayesian additive regression trees (BART) to flexibly specify separate mean models for the potential outcomes and latent stratum membership. In the analysis of the ARMA trial, we found that the low tidal volume treatment had an overall benefit for participants sustaining acute lung injuries on the outcome of time to returning home but substantial heterogeneity in treatment effects among the always-survivors, driven most strongly by biologic sex and the alveolar-arterial oxygen gradient at baseline (a physiologic measure of lung function and degree of hypoxemia). These findings illustrate how the proposed methodology could guide the prognostic enrichment of future trials in the field.

Women’s educational attainment and contraceptive prevalence are two mechanisms identified as having an accelerating effect on fertility decline and that can be directly impacted by policy. Quantifying the potential accelerating effect of education and family planning policies on fertility decline in a probabilistic way is of interest to policymakers, particularly in high-fertility countries. We propose a conditional Bayesian hierarchical model for projecting fertility, given education and family planning policy interventions. To illustrate the effect policy changes could have on future fertility, we create probabilistic projections of fertility that condition on scenarios such as achieving the sustainable development goals (SDGs) for universal secondary education and universal access to family planning by 2030.

Intense wildfires impact nature, humans, and society, causing catastrophic damage to property and the ecosystem as well as the loss of life. Forecasting wildfire front propagation and understanding the behavior of wildfire propagation within a formal uncertainty quantification framework are essential in order to support fire fighting efforts and plan evacuations. The level set method has been widely used to analyze the change in surfaces, shapes, and boundaries. In particular, a signed distance function used in level set methods can readily be interpreted to represent complicated boundaries and their changes in time. While there is substantial literature on the level set method in wildfire applications, these implementations have relied on a heavily-parameterized formula for the rate of spread. These implementations have not typically considered uncertainty quantification, incorporated data-driven learning, nor summarized the effect of the environmental covariates. Here we present a Bayesian spatio-temporal dynamic model, based on level sets, which can be utilized for inference and forecasting the boundary of interest in the presence of uncertain data and lack of knowledge about the boundary velocity. The methodology relies on both a mechanistically-motivated dynamic model for level sets and a stochastic spatio-temporal dynamic model for the front velocity. We show the effectiveness of our method via simulation and with forecasting the fire front boundary evolution of two classic California megafires—the 2017–2018 Thomas fire and the 2017 Haypress fire.

Traffic state prediction is an essential component and an underlying backbone of intelligent transportation systems, especially in the context of smart city framework. Its significance is mainly twofold in modern transportation systems: supporting advanced traffic operations and management for highways and urban road networks to mitigate traffic congestion and enabling individual drivers with connected vehicles in the traffic system to dynamically optimize their routes to improve travel time. Traffic state prediction with interval-based pointwise methods at 15-minute or hourly intervals is common in traffic literature. However, because traffic dynamics are a continuous process over time, the discrete-time pointwise methods for traffic prediction at a fixed time interval hardly meet the advanced demands of continuous prediction in modern transportation systems. To close the gap, we propose functional approaches to intraday and day-by-day continuous-time prediction for traffic volume. This research focuses on network-level traffic flow predictions concurrently for all locations of interest. Two functional approaches are introduced, namely, the network-integrated functional time-series model and the functional neural network model. With functional approaches a 24-hour intraday traffic profile is modeled as a functional curve over time, and sequences of historical traffic curves are used to predict traffic curves for near future days in a row and multiple locations of interest. We also include the functional varying coefficient model, Sparse VAR and traditional AR models in the comparative study; empirical results show that the network-integrated functional time-series model outperforms other approaches in terms of the accuracy of predictions at network-scale.

Multitask learning is frequently used to model a set of related response variables from the same set of features, improving predictive performance and modeling accuracy relative to methods that handle each response variable separately. Despite the potential of multitask learning to yield more powerful inference than single-task alternatives, prior work in this area has largely omitted uncertainty quantification. Our focus in this paper is a common multitask problem in neuroimaging, where the goal is to understand the relationship between multiple cognitive task scores (or other subject-level assessments) and brain connectome data collected from imaging. We propose a framework for selective inference to address this problem, with the flexibility to: (i) jointly identify the relevant predictors for each task through a sparsity-inducing penalty and (ii) conduct valid inference in a model based on the estimated sparsity structure. Our framework offers a new conditional procedure for inference, based on a refinement of the selection event that yields a tractable selection-adjusted likelihood. This gives an approximate system of estimating equations for maximum likelihood inference, solvable via a single convex optimization problem, and enables us to efficiently form confidence intervals with approximately the correct coverage. Applied to both simulated data and data from the Adolescent Brain Cognitive Development (ABCD) study, our selective inference methods yield tighter confidence intervals than commonly used alternatives, such as data splitting. We also demonstrate through simulations that multitask learning with selective inference can more accurately recover true signals than single-task methods.

Image-on-scalar regression has been a popular approach to modeling the association between brain activities and scalar characteristics in neuroimaging research. The associations could be heterogeneous across individuals in the population, as indicated by recent large-scale neuroimaging studies, for example, the Adolescent Brain Cognitive Development (ABCD) Study. The ABCD data can inform our understanding of heterogeneous associations and how to leverage the heterogeneity and tailor interventions to increase the number of youths who benefit. It is of great interest to identify subgroups of individuals from the population such that: (1) within each subgroup the brain activities have homogeneous associations with the clinical measures; (2) across subgroups the associations are heterogeneous, and (3) the group allocation depends on individual characteristics. Existing image-on-scalar regression methods and clustering methods cannot directly achieve this goal. We propose a latent subgroup image-on-scalar regression model (LASIR) to analyze large-scale, multisite neuroimaging data with diverse sociodemographics. LASIR introduces the latent subgroup for each individual and group-specific, spatially varying effects, with an efficient stochastic expectation maximization algorithm for inferences. We demonstrate that LASIR outperforms existing alternatives for subgroup identification of brain activation patterns with functional magnetic resonance imaging data via comprehensive simulations and applications to the ABCD study. We have released our reproducible codes for public use with the software package available on Github.

Many genetic studies contain rich information on longitudinal phenotypes that require powerful analytical tools for optimal analysis. Genetic analysis of longitudinal data that incorporates temporal variation is important for understanding the genetic architecture and biological variation of complex diseases. Most of the existing methods assume that the contribution of genetic variants is constant over time and fail to capture the dynamic pattern of disease progression. However, the relative influence of genetic variants on complex traits fluctuates over time. In this study, we propose a retrospective varying coefficient mixed model association test, RVMMAT, to detect time-varying genetic effect on longitudinal binary traits. We model dynamic genetic effect using smoothing splines, estimate model parameters by maximizing a double penalized quasi-likelihood function, design a joint test using a Cauchy combination method, and evaluate statistical significance via a retrospective approach to achieve robustness to model misspecification. Through simulations we illustrated that the retrospective varying-coefficient test was robust to model misspecification under different ascertainment schemes and gained power over the association methods assuming constant genetic effect. We applied RVMMAT to a genome-wide association analysis of longitudinal measure of hypertension in the Multi-Ethnic Study of Atherosclerosis. Pathway analysis identified two important pathways related to G-protein signaling and DNA damage. Our results demonstrated that RVMMAT could detect biologically relevant loci and pathways in a genome scan and provided insight into the genetic architecture of hypertension.

Estimation of heterogeneous treatment effects (HTE) is of prime importance in many disciplines, from personalized medicine to economics among many others. Random forests have been shown to be a flexible and powerful approach to HTE estimation in both randomized trials and observational studies. In particular “causal forests” introduced by Athey, Tibshirani and Wager (Ann. Statist. 47 (2019) 1148–1178), along with the $\mathsf{R}$ implementation in package grf were rapidly adopted. A related approach, called “model-based forests” that is geared toward randomized trials and simultaneously captures effects of both prognostic and predictive variables, was introduced by Seibold, Zeileis and Hothorn (Stat. Methods Med. Res. 27 (2018) 3104–3125) along with a modular implementation in the $\mathsf{R}$ package model4you.

Neither procedure is directly applicable to the estimation of individualized predictions of excess postpartum blood loss caused by a cesarean section in comparison to vaginal delivery. Clearly, randomization is hardly possible in this setup, and thus model-based forests lack clinical trial data to address this question. On the other hand, the skewed and interval-censored postpartum blood loss observations violate assumptions made by causal forests. Here we present a tailored model-based forest for skewed and interval-censored data to infer possible predictive prepartum characteristics and their impact on excess postpartum blood loss caused by a cesarean section.

As a methodological basis, we propose a unifying view on causal and model-based forests that goes beyond the theoretical motivations and investigates which computational elements make causal forests so successful and how these can be blended with the strengths of model-based forests. To do so, we show that both methods can be understood in terms of the same parameters and model assumptions for an additive model under ${\mathit{L}}_2$ loss. This theoretical insight allows us to implement several flavors of “model-based causal forests” and dissect their different elements in silico.

The original causal forests and model-based forests are compared with the new blended versions in a benchmark study exploring both randomized trials and observational settings. In the randomized setting, both approaches performed akin. If confounding was present in the data-generating process, we found local centering of the treatment indicator with the corresponding propensities to be the main driver for good performance. Local centering of the outcome was less important and might be replaced or enhanced by simultaneous split selection with respect to both prognostic and predictive effects. This lays the foundation for future research combining random forests for HTE estimation with other types of models.

The algorithms used for the optimal management of an ambulance fleet require an accurate description of the spatiotemporal evolution of the emergency events. In the last years, several authors have proposed sophisticated statistical approaches to forecast ambulance dispatches, typically modelling the data as a point pattern occurring on a planar region. Nevertheless, ambulance interventions can be more appropriately modelled as a realisation of a point process occurring on a linear network. The constrained spatial domain raises specific challenges and unique methodological problems that cannot be ignored when developing a proper statistical approach. Hence, this paper proposes a spatiotemporal model to analyse ambulance dispatches focusing on the interventions that occurred in the road network of Milan (Italy) from 2015 to 2017. We adopt a nonseparable first-order intensity function with spatial and temporal terms. The temporal dimension is estimated semiparametrically using a Poisson regression model, while the spatial dimension is estimated nonparametrically using a network kernel function. A set of weights is included in the spatial term to capture space-time interactions, inducing nonseparability in the intensity function. A series of tests show that our approach successfully models the ambulance interventions and captures the space-time patterns more accurately than planar or separable point process models.

Reliable estimates of volatility and correlation are fundamental in economics and finance for understanding the impact of macroeconomics events on the market and guiding future investments and policies. Dependence across financial returns is likely to be subject to sudden structural changes, especially in correspondence with major global events, such as the COVID-19 pandemic. In this work we are interested in capturing abrupt changes over time in the conditional dependence across U.S. industry stock portfolios, over a time horizon that covers the COVID-19 pandemic. The selected stocks give a comprehensive picture of the U.S. stock market. To this end, we develop a Bayesian multivariate stochastic volatility model based on a time-varying sequence of graphs capturing the evolution of the dependence structure. The model builds on the Gaussian graphical models and the random change points literature. In particular, we treat the number, the position of change points, and the graphs as object of posterior inference, allowing for sparsity in graph recovery and change point detection. The high dimension of the parameter space poses complex computational challenges. However, the model admits a hidden Markov model formulation. This leads to the development of an efficient computational strategy, based on a combination of sequential Monte-Carlo and Markov chain Monte-Carlo techniques. Model and computational development are widely applicable, beyond the scope of the application of interest in this work.

This paper estimates the curvature of the Earth, defined as one over its radius, without relying on physical measurements. The orthodox model states that the Earth is (nearly) spherical with a curvature of $\mathit{\pi }/20,000\text{km}$. By contrast, the heterodox flat-Earth model stipulates a curvature of zero. Abstracting from the well-worn arguments for and against both models, rebuttals and counter-rebuttals ad infinitum, we propose a novel statistical methodology based on verifiable flight times along regularly scheduled commercial airline routes; this methodology allows for both estimating and making inference for Earth’s curvature. In particular, a formal hypothesis test resolutely rejects the flat-Earth model, whereas it does not reject the orthodox spherical-Earth model.

Air pollution surveillance is critically important for public health. One air pollutant, ozone, is extremely challenging to analyze properly, as it is a secondary pollutant caused by complex chemical reactions in the air and does not emit directly into the atmosphere. Numerous environmental studies confirm that ozone concentration levels are associated with meteorological conditions, and long-term exposure to high ozone concentration levels is associated with the incidence of many diseases, including asthma, respiratory, and cardiovascular diseases. Thus, it is important to develop an air pollution surveillance system to collect both air pollution and meteorological data and monitor the data continuously over time. To this end, statistical process control (SPC) charts provide a major statistical tool. But most existing SPC charts are designed for cases when the in-control (IC) process observations at different times are assumed to be independent and identically distributed. The air pollution and meteorological data would not satisfy these conditions due to serial data correlation, high dimensionality, seasonality, and other complex data structure. Motivated by an application to monitor the ground ozone concentration levels in the Houston–Galveston–Brazoria (HGB) area, we developed a new process monitoring method using principal component analysis and sequential learning. The new method can accommodate high dimensionality, time-varying IC process distribution, serial data correlation, and nonparametric data distribution. It is shown to be a reliable analytic tool for online monitoring of air quality.

Sensor arrays are often used to identify chemicals by measuring properly chosen chemical interactions. Machine learning techniques are of vital importance to accurately recognize a chemical based on the sensor array measurements. However, sensor array data often take the form of matrices (i.e, two-way tensors), and the concentration levels may have a complex impact on the measurements. Hence, existing linear and/or vector classification methods may be inadequate for sensor array data. In this article we propose a novel tensor mixture discriminant analysis (TMDA) model carefully tailored for the classification of sensor array data. We model the distribution of each chemical by a mixture of tensor normal distributions. TMDA leverages the tensor structure for better estimation and prediction, while the mixed tensor normal component accounts for the possibly varying concentration levels. The TMDA model can also be viewed as an approximation of the potentially nonnormal measurements. An efficient expectation-maximization algorithm is developed to fit the TMDA model. The application of TMDA on two sensor array datasets demonstrates its superior performance to many popular competitors.

Pollutant exposure during gestation is a known and adverse factor for birth and health outcomes. However, the links between prenatal air pollution exposures and educational outcomes are less clear, in particular, the critical windows of susceptibility during pregnancy. Using a large cohort of students in North Carolina, we study the link between prenatal daily ${\text{PM}}_{2.5}$ exposure and fourth end-of-grade reading scores. We develop and apply a locally adaptive and highly scalable Bayesian regression model for scalar responses with functional and scalar predictors. The proposed model pairs a B-spline basis expansion with dynamic shrinkage priors to capture both smooth and rapidly-changing features in the regression surface. The model is accompanied by a new decision analysis approach for functional regression that extracts the critical windows of susceptibility and guides the model interpretations. These tools help to identify and address broad limitations with the interpretability of functional regression models. Simulation studies demonstrate more accurate point estimation, more precise uncertainty quantification, and far superior window selection than existing approaches. Leveraging the proposed modeling, computational, and decision analysis framework, we conclude that prenatal ${\text{PM}}_{2.5}$ exposure during early and late pregnancy is most adverse for fourth end-of-grade reading scores.

For valid scientific discoveries, it is fundamental to evaluate whether research findings are replicable across different settings. While large-scale replication projects across broad research topics are not feasible, systematic reviews and meta-analyses (SRMAs) offer viable alternatives to assess replicability. Due to subjective inclusion and exclusion of studies, SRMAs may contain nonreplicable study findings. However, there is no consensus on rigorous methods to assess the replicability of SRMAs or to explore sources of nonreplicability. Nonreplicability is often misconceived as high heterogeneity. This article introduces a new measure, the externally standardized residuals from a leave-m-studies-out procedure, to quantify replicability. It not only measures the impact of nonreplicability from unknown sources on the conclusion of an SRMA but also differentiates nonreplicability from heterogeneity. A new test statistic for replicability is derived. We explore its asymptotic properties and use extensive simulations and real data to illustrate this measure’s performance. We conclude that replicability should be routinely assessed for all SRMAs and recommend sensitivity analyses, once nonreplicable study results are identified in an SRMA.

Respondent-driven sampling (RDS) is used throughout the world to estimate prevalence and population size for hidden populations. Although RDS is an effective method for enrolling people from key populations in studies, it relies on a partially unknown sampling mechanism, and thus each individual’s inclusion probability is unknown. Current estimators for population prevalence, population size, and other outcomes rely on a participant’s network size (degree) to approximate their inclusion probability in the sample from the networked population. However, in most RDS studies, a participant’s network size is attained via a self-report and is subject to many types of misreporting and bias. Because design-based inclusion probabilities cannot be exactly computed, we instead use the term visibility to describe how likely a person is to be selected to participate in the study. The commonly used successive sampling population size estimation (SS-PSE) framework to estimate population sizes from RDS data relies on self-reported network sizes in the model for the sampling mechanism. We propose an enhancement of the SS-PSE framework that adds a measurement error model for visibility used in place of the self-reported network size and a model for the number of recruits an individual can enroll. Inferred visibilities are a way to smooth the degree distribution and bring in outliers as well as a mechanism to deal with missing and invalid network sizes. We demonstrate the performance of visibility SS-PSE on three populations from Kosovo sampled in 2014 using RDS. We also discuss how the visibility modeling framework could be extended to prevalence estimation.

This paper is motivated by the joint analysis of genetic, imaging, and clinical (GIC) data collected in the Alzheimer’s Disease Neuroimaging Initiative (ADNI) study. We propose a partially functional linear regression (PFLR) framework to map high-dimensional GIC-related pathways for Alzheimer’s disease (AD). We develop a joint model selection and estimation procedure by embedding imaging data in the reproducing kernel Hilbert space and imposing the ${\ell }_0$ penalty for the coefficients of genetic variables. We apply the proposed method to the ADNI dataset to identify important features from tens of thousands of genetic polymorphisms (reduced from millions using a preprocessing step) and study the effects of a certain set of informative genetic variants and the baseline hippocampus surface on 13 future cognitive scores. We also explore the shared and distinct heritability patterns of these cognitive scores. Analysis results suggest that both the hippocampal and genetic data have heterogeneous effects on different scores, with the trend that the value of both hippocampi are negatively associated with the severity of cognition deficits. Polygenic effects are observed for all the thirteen cognitive scores. The well-known APOE4 genotype only explains a small part of the cognitive function. Shared genetic etiology exists; however, greater genetic heterogeneity exists within disease classifications after accounting for the baseline diagnosis status. These analyses are useful in further investigation of functional mechanisms for AD progression.

Risk-adjusted quality measures are used to evaluate healthcare providers with respect to national norms while controlling for factors beyond their control. Existing healthcare provider profiling approaches typically assume that the between-provider variation in these measures is entirely due to meaningful differences in quality of care. However, in practice, much of the between-provider variation will be due to trivial fluctuations in healthcare quality or unobservable confounding risk factors. If these additional sources of variation are not accounted for, conventional methods will disproportionately identify larger providers as outliers, even though their departures from the national norms may not be “extreme” or clinically meaningful. Motivated by efforts to evaluate the quality of care provided by transplant centers, we develop a composite evaluation score based on a novel individualized empirical null method, which robustly accounts for overdispersion due to unobserved risk factors, models the marginal variance of standardized scores as a function of the effective sample size, and only requires the use of publicly-available center-level statistics. The evaluations of United States kidney transplant centers based on the proposed composite score are substantially different from those based on conventional methods. Simulations show that the proposed empirical null approach more accurately classifies centers in terms of quality of care, compared to existing methods.

When searching for exoplanets, one wants to count how many planets orbit a given star, and to determine what their characteristics are. If the estimated planet characteristics are too far from those of a planet truly present, this should be considered as a false detection. This setting is a particular instance of a general one: aiming to retrieve parametric components in a dataset corrupted by nuisance signals, with a certain accuracy on their parameters. We exhibit a detection criterion minimizing false and missed detections, either as a function of their relative cost or when the expected number of false detections is bounded. If the components can be separated in a technical sense discussed in detail, the optimal detection criterion is a posterior probability obtained as a by-product of Bayesian evidence calculations. Optimality is guaranteed within a model, and we introduce model criticism methods to ensure that the criterion is robust to model errors. We show on two simulations emulating exoplanet searches that the optimal criterion can significantly outperform other criteria. Finally, we show that our framework offers solutions for the identification of components of mixture models and Bayesian false discovery rate control when hypotheses are not discrete.

We consider matching problems where the goal is to determine whether two observations randomly drawn from a population with multiple (sub)groups are from the same (sub)group. This is a key question in forensic science, where items with unidentified origins from suspects and crime scenes are compared to objects from a known set of sources to see if they originated from the same source. We derive the optimal matching rule under known density functions of data that minimizes the decision error probabilities. Empirically, the proposed matching rule is computed by plugging parametrically estimated density functions using training data into the formula of the optimal matching rule. The connections between the optimal matching rule and existing methods in forensic science are explained. In particular, we contrast the optimal matching rule to classification and also compare it to a score-based approach that relies on similarity features extracted from paired items. Numerical simulations are conducted to evaluate the proposed method and show that it outperforms the existing methods in terms of a higher ROC curve and higher power to identify matched pairs of items. We also demonstrate the utility of the proposed method by applying it to a real forensic data analysis of glass fragments.

Collective efficacy—the capacity of communities to exert social control toward the realization of their shared goals—is a foundational concept in the urban sociology and neighborhood effects literature. Traditionally, empirical studies of collective efficacy use large sample surveys to estimate collective efficacy of different neighborhoods within an urban setting. Such studies have demonstrated an association between collective efficacy and local variation in community violence, educational achievement, and health. Unlike traditional collective efficacy measurement strategies, the Adolescent Health and Development in Context (AHDC) Study implemented a new approach, obtaining spatially-referenced, place-based ratings of collective efficacy from a representative sample of individuals residing in Columbus, OH. In this paper we introduce a novel nonstationary spatial model for interpolation of the AHDC collective efficacy ratings across the study area, which leverages administrative data on land use. Our constructive model specification strategy involves dimension expansion of a latent spatial process and the use of a filter defined by the land-use partition of the study region to connect the latent multivariate spatial process to the observed ordinal ratings of collective efficacy. Careful consideration is given to the issues of parameter identifiability, computational efficiency of an MCMC algorithm for model fitting, and fine-scale spatial prediction of collective efficacy.

The rapid growth of high-throughput genomic and epigenomic data enables exploration of biological mechanisms underlying diseases causing processes beyond traditional association studies. Using the causal mediation analysis framework, we develop the kernel machine difference (KMD) method, which provides a testing procedure for detecting the mediation effects of a set of mediators, for example, the DNA methylation probes within a region or a gene. Our method extends the difference method in single mediator analysis to jointly model the mediatory role of the methylation of multiple neighboring probes, as they often work together in a collaborative fashion. Kernel machine regression is employed to accommodate flexible parametric and nonparametric effects of multiple mediators on the outcome and to allow for robust testing for the joint natural indirect effect (NIE). The proposed testing procedure does not require explicit modeling of the dependence of multiple mediators on exposure and confounders and the correlation among multiple mediators. It hence provides a robust and computationally efficient tool, especially for genomic regions with moderate to high-dimensional probes. We evaluate the performance of the proposed test using extensive simulations and demonstrate its gain in robustness and power when the effects are nonlinear. We apply the proposed test to the analysis of the epigenome-wide Normative Aging Study (NAS) to investigate the mediatory role of DNA methylation in the causal pathway between smoking behavior and lung function.

In this study we present a data-driven method called false negative control (FNC) screening to address the challenge of detecting weak signals in underpowered genome-wide association studies (GWASs), where true signals are often obscured by a large amount of noise. Our approach focuses on controlling false negatives and efficiently regulates the proportion of false negatives at a user-specified level in realistic settings with arbitrary covariance dependence between variables. We calibrate overall dependence using a parameter that aligns with the existing phase diagram in high-dimensional sparse inference, allowing us to asymptotically explicate the joint effect of covariance dependence, signal sparsity, and signal intensity on the proposed method. Our new phase diagram shows that FNC screening can efficiently select a set of candidate variables to retain a high proportion of signals, even when the signals are not individually separable from noise. We compare the performance of FNC screening to several existing methods in simulation studies, and the proposed method outperforms the others in adapting to a user-specified false negative control level. Moreover, we apply FNC screening to 145 GWAS datasets, obtained from the UK Biobank, and demonstrate a substantial increase in power to retain true signals for downstream analyses.

In scientific studies involving analyses of multivariate data, basic but important questions often arise for the researcher: Is the sample exchangeable, meaning that the joint distribution of the sample is invariant to the ordering of the units? Are the features independent of one another, or perhaps the features can be grouped so that the groups are mutually independent? In statistical genomics these considerations are fundamental to downstream tasks such as demographic inference and the construction of polygenic risk scores. We propose a nonparametric approach, which we call the V test, to address these two questions, namely, a test of sample exchangeability given dependency structure of features and a test of feature independence given sample exchangeability. Our test is conceptually simple, yet fast and flexible. It controls the Type I error across realistic scenarios and handles data of arbitrary dimensions by leveraging large-sample asymptotics. Through extensive simulations and a comparison against unsupervised tests of stratification based on random matrix theory, we find that our test compares favorably in various scenarios of interest. We apply the test to data from the 1000 Genomes Project, demonstrating how it can be employed to assess exchangeability of the genetic sample or find optimal linkage disequilibrium (LD) splits for downstream analysis. For exchangeability assessment we find that removing rare variants can substantially increase the p-value of the test statistic. For optimal LD splitting, the V test reports different optimal splits than previous approaches not relying on hypothesis testing. Software for our methods is available in R (CRAN: $\mathsf{flintyR}$) and Python (PyPI: $\mathsf{flintyPy}$).

Knowledge networks, such as the healthcare delivery network (HDN), describing relationships among different medical encounters, are useful summaries of state-of-art medical knowledge. The increasing availability of longitudinal electronic health records (EHR) data promises a rich data source for learning HDN. Most existing methods for inferring knowledge networks are based on cooccurrence patterns that do not account for temporal effects or patient-level heterogeneity. In this article, building upon the multivariate Hawkes process (mvHP), we propose a flexible covariate-adjusted random effects (CARE) mvHP modeling strategy for HDN construction. Our model allows for patient-specific time-varying background intensity functions via random effects, which can also adjust for effects of important covariates. We adopt a penalized approach to select fixed effects, yielding a sparse network structure, and to remove unnecessary random effects from the model. Through extensive simulation studies, we show that our proposed method performs well in recovering the network structure and that it is essential to account for patient heterogeneities. We further illustrate our CARE mvHP method in an EHR study of type 2 diabetes patients to learn an HDN for these patients and demonstrate that our results are consistent with current clinical practice in healthcare systems.