The cyclical and heterogeneous nature of many substance use disorders highlights the need to adapt the type and/or the dose of treatment to accommodate the specific and changing needs of individuals. The Adaptive Treatment for Alcohol and Cocaine Dependence study (ENGAGE) is a sequential multiple assignment randomized trial (SMART) that provided longitudinal data for constructing dynamic treatment regimes (DTRs) to improve patients’ engagement in therapy. However, the high rate of noncompliance and lack of analytic tools to account for noncompliance has impeded researchers from using the data to achieve the main goal of the trial; namely, construction of individually tailored DTRs. We address this by defining our target parameter as the mean outcome under different DTRs for potential compliance strata and propose a marginal structural model with principal stratification to estimate this quantity. We model the principal strata using a Bayesian semiparametric approach. An important feature of our work is that we consider partial rather than binary compliance strata, which is more relevant in longitudinal studies. We assess the performance of our method through simulation. We illustrate its application on ENGAGE and demonstrate the optimal DTRs depend on compliance strata compared with ignoring compliance information as in intention-to-treat analyses.

Data-driven network models to measure systemic risk in the financial sector and identify “too-connected-to-fail” institutions are becoming increasingly common in financial applications. Existing statistical methods for building such networks either take a pairwise approach of fitting many bivariate models or a system-wide approach of fitting penalized regression models. The former strategy is prone to large false positive selection, while the latter suffers from shrinkage bias and lack of formal inference machinery. These issues are accentuated in small sample, low signal-to-noise settings common in financial data. Building up on recent advances in high-dimensional inference, we propose debiased lasso Penalized Vector Autoregression (DLVAR), a method for building financial networks that addresses these limitations. Our empirical analysis highlights the importance of debiasing in a way that increases power of the algorithm in finite samples. We also provide formal inference guarantees of Granger causality tests in high-dimension to justify our method. We apply DLVAR to the stock returns of U.S. large financial institutions covering the period 1990–2021 and illustrate its usefulness in detecting systemically risky periods and institutions, especially during the Great Financial Crisis of 2008–2009 and the most recent Covid-19 related market shock.

The field of retail analytics has been transformed by the availability of rich data, which can be used to perform tasks such as demand forecasting and inventory management. However, one task which has proved more challenging is the forecasting of demand for products which exhibit very few sales. The sparsity of the resulting data limits the degree to which traditional analytics can be deployed. To combat this, we represent sales data as a structured sparse multivariate point process, which allows for features such as autocorrelation, cross-correlation, and temporal clustering, known to be present in sparse sales data. We introduce a Bayesian point process model to capture these phenomena, which includes a hurdle component to cope with sparsity and an exciting component to cope with temporal clustering within and across products. We then cast this model within a Bayesian hierarchical framework, to allow the borrowing of information across different products, which is key in addressing the data sparsity per product. We conduct a detailed analysis, using real sales data, to show that this model outperforms existing methods in terms of predictive power, and we discuss the interpretation of the inference.

Mendelian randomization (MR) is a popular method in epidemiology and genetics that uses genetic variation as instrumental variables for causal inference. Existing MR methods usually assume most genetic variants are valid instrumental variables that identify a common causal effect. There is a general lack of awareness that this effect homogeneity assumption can be violated when there are multiple causal pathways involved, even if all the instrumental variables are valid. In this article we introduce a latent mixture model MR-Path that groups instruments that yield similar causal effect estimates together. We develop a Monte Carlo EM algorithm to fit this mixture model, derive approximate confidence intervals for uncertainty quantification, and adopt a modified Bayesian Information Criterion (BIC) for model selection. We verify the efficacy of the Monte Carlo EM algorithm, confidence intervals, and model selection criterion using numerical simulations. We identify potential mechanistic heterogeneity when applying our method to estimate the effect of high-density lipoprotein cholesterol on coronary heart disease and the effect of adiposity on type II diabetes.

The online chauffeured service demand (OCSD) research is an exploratory market study of designated driver services in China. Researchers are interested in the influencing factors of chauffeured service adoption and usage and have collected relevant data using a self-reported questionnaire. As self-reported count measure data is typically inflated, there exist challenges to its validity, which may bias estimation and increase error in empirical research. Motivated by the analysis of self-reported data with multiple inflated values, we propose a novel approach to simultaneously achieve data-driven inflated value selection and identification of important influencing factors. In particular, the regularization technique is applied to the mixing proportions of inflated values and the regression parameters to obtain shrinkage estimates. We analyze the OCSD data with the proposed approach, deriving insights into the determinants impacting service demand. The proper interpretations and implications contribute to service promotion and related policy optimization. Extensive simulation studies and consistent asymptotic properties further establish the effectiveness of the proposed approach.

Traditionally, readability prediction has relied on readability formulas, which are based on shallow text characteristics such as average word and sentence length. With recent advances in text mining and natural language processing, more complex text properties can be incorporated into readability prediction models, with papers in the literature suggesting to use up to 200 features for predicting text readability. However, many of the features generated using natural language processing tools are highly correlated and can be thought to measure similar latent text properties. When dealing with a high-dimensional space of correlated features, removing the redundant variables has two advantages: (1) improving interpretability and (2) increasing the predictive power of the model. In this paper we propose an ordinal version of the averaged lasso, which combines hierarchical clustering with the lasso, in order to identify relevant features for readability prediction. We illustrate the approach on two corpora and show improved prediction accuracy when benchmarking against a set of competing models. The annotated corpora as well as the steps necessary for feature creation are freely available as $\mathsf{R}$ packages, thus allowing the obtained results to be directly incorporated into a readability estimation pipeline.

In biomedical research the analysis of disease prevalence is of critical importance. While most of the existing prevalence studies focus on individual diseases, there has been increasing effort that jointly examines the prevalence values and their trends of multiple diseases. Such joint analysis can provide valuable insights not shared by individual-disease analysis. A critical limitation of the existing analysis is that there is a lack of attention to existing information, which has been accumulated through a large number of studies and can be valuable especially when there are a large number of diseases but the number of prevalence values for a specific disease is limited. In this study we conduct the functional clustering analysis of prevalence trends for a large number of diseases. A novel approach based on the penalized fusion technique is developed to incorporate information mined from published articles. It is innovatively designed to take into account that such information may not be fully relevant or correct. Another significant development is that statistical properties are rigorously established. Simulation is conducted and demonstrates its competitive performance. In the analysis of data from Taiwan NHIRD (National Health Insurance Research Database), new and interesting findings that differ from the existing ones are made.

We consider the problem of predicting breast cancer risk using mammogram imaging data where the dimension of pixels greatly exceed the number of individuals in the cohort. The functional partial least squares (FPLS) is a popular dimensional reduction method in constructing latent explanatory components using linear combinations of the original predictor variables. While FPLS with scalar responses has been studied in the literature, the presence of right censoring under the survival framework poses challenges in modeling and estimation. Given several different representations for PLS with Cox regression in the literature, we unify and extend three formulations to deal with right censoring, that is, reweighing, mean imputation, and deviance residuals to the functional setting in this paper. We empirically investigate and compare the performance of the three proposed FPLS frameworks in the context of imaging predictor via intensive simulation studies. The proposed methods are applied to the Joanne Knight Breast Health Cohort where we show increased model discriminatory performance under the FPLS framework compared to competing models.

Numerical air quality models, such as the Community Multiscale Air Quality (CMAQ) system, play a critical role in characterizing pollution levels at fine spatial and temporal scales. The model outputs, however, tend to systematically over- or underestimate the real pollutant concentrations. In this study we propose a Bayesian hierarchical dynamic model to calibrate large-scale grid-level CMAQ model outputs using data from other sources, especially point-level observations from sparsely located monitoring stations. In our model a stochastic integro-differential equation (IDE) is implemented to account for space-time interactions of air pollutants. To better approximate the spatial pattern of pollutants, we employ nonregular meshes to discretize IDEs. A spatial partitioning procedure is embedded to improve the scalability of the approach for very large meshes. An algorithm based on variational Bayes and ensemble Kalman smoother is developed to accelerate the parameter estimation and calibration procedure. We apply the proposed approach to calibrate CMAQ outputs for China’s Beijing–Tianjin–Hebei region. In contrast to existing methods, the proposed approach captures space-time interactions, produces more accurate calibration results, and operates at a higher computational efficiency. A reanalysis dataset is also adopted to demonstrate the effectiveness and efficiency of our approach to large spatial data.

Brain connectome analysis commonly compresses high-resolution brain scans (typically composed of millions of voxels) down to only hundreds of regions of interest (ROIs) by averaging within-ROI signals. This significant dimension reduction improves computational speed and the morphological properties of anatomical structures; however, it comes at the cost of substantial losses in spatial specificity and sensitivity, especially when the signals exhibit high within-ROI heterogeneity. Oftentimes, abnormally expressed functional connectivity (FC) between a pair of ROIs, caused by a brain disease, is primarily driven by only small subsets of voxel pairs within the ROI pair. This article proposes a new network method for the detection of voxel-pair-level neural dysconnectivity with spatial constraints. Specifically, focusing on an ROI pair, our model aims to extract dense subareas that contain aberrant voxel-pair connections while ensuring that the involved voxels are spatially contiguous. In addition, we develop subcommunity-detection algorithms to realize the model, and we justify the consistency of these algorithms. Comprehensive simulation studies demonstrate our method’s effectiveness in reducing the false-positive rate while increasing statistical power, detection replicability, and spatial specificity. We apply our approach to reveal: (i) disrupted voxelwise FC patterns related to nicotine addiction between the basal ganglia, hippocampus, and insular gyrus in 3269 participants using UK Biobank data; (ii) voxelwise schizophrenia-altered FC patterns within the salience and temporal-thalamic network in 330 participants in a schizophrenia study. The detected results align with previous medical findings but include improved localized information.

Sex difference in allele frequency is an emerging topic that is crucial to our understanding of data quality and features, particularly when it comes to the largely overlooked X chromosome. To detect sex differences in allele frequency for both X chromosomal and autosomal variants, the existing method is conservative when applied to samples from multiple ancestral populations. Additionally, it remains unexplored whether the sex difference in allele frequency varies between populations, which is important for transancestral genetic studies. To answer these questions, we thus developed a novel, retrospective regression-based testing framework that led to interpretable and easy-to-implement solutions. We then applied the proposed methods to the high-coverage whole genome sequence data of the 1000 Genomes Project, robustly analyzing all samples available from the five super-populations. We had 97 novel findings by recognizing and modelling ancestral differences. Finally, we replicated the specific findings and overall conclusion using the gnomAD v3.1.2 data.

Understanding cause-specific mortality rates is crucial for monitoring population health and designing public health interventions. Worldwide, two-thirds of deaths do not have a cause assigned. Verbal autopsy (VA) is a well-established tool to collect information describing deaths outside of hospitals by conducting surveys to caregivers of a deceased person. It is routinely implemented in many low- and middle-income countries. Statistical algorithms to assign cause of death using VAs are typically vulnerable to the distribution shift between the data used to train the model and the target population. This presents a major challenge for analyzing VAs, as labeled data are usually unavailable in the target population. This article proposes a latent class model framework for VA data (LCVA) that jointly models VAs collected over multiple heterogeneous domains, assigns causes of death for out-of-domain observations and estimates cause-specific mortality fractions for a new domain. We introduce a parsimonious representation of the joint distribution of the collected symptoms using nested latent class models and develop a computationally efficient algorithm for posterior inference. We demonstrate that LCVA outperforms existing methods in predictive performance and scalability. Supplementary Material and reproducible analysis codes are available online. The R package LCVA implementing the method is available on GitHub ( https://github.com/richardli/LCVA).

Local field potentials (LFPs) are signals that measure electrical activities in localized cortical regions and are collected from multiple tetrodes implanted across a patch on the surface of cortex. Hence, they can be treated as multigroup functional data, where the trajectories collected across temporal epochs from one tetrode are viewed as a group of functions. In many cases multitetrode LFP trajectories contain both global variation patterns (which are shared by most groups, due to signal synchrony) and idiosyncratic variation patterns (common only to a small subset of groups), and such structure is very informative to the data mechanism. Therefore, one goal in this paper is to develop an efficient algorithm that is able to capture and quantify both global and idiosyncratic features. We develop the novel filtrated common functional principal components (filt-fPCA) method, which is a novel forest-structured fPCA for multigroup functional data. A major advantage of the proposed filt-fPCA method is its ability to extract the common components in a flexible “multiresolution” manner. The proposed approach is highly data-driven, and no prior knowledge of “ground-truth” data structure is needed, making it suitable for analyzing complex multigroup functional data. In addition, the filt-fPCA method is able to produce parsimonious, interpretable, and efficient functional reconstruction (low reconstruction error) for multigroup functional data with orthonormal basis functions. Here the proposed filt-fPCA method is employed to study the impact of a shock (induced stroke) on the synchrony structure of rat brain. The proposed filt-fPCA is general and inclusive that can be readily applied to analyze any multigroup functional data, such as multivariate functional data, spatial-temporal data, and longitudinal functional data.

Bulk transcriptomics in tissue samples reflects the average expression levels across different cell types and is highly influenced by cellular fractions. As such, it is critical to estimate cellular fractions to both deconfound differential expression analyses and infer cell type-specific differential expression. Since experimentally counting cells is infeasible in most tissues and studies, in silico cellular deconvolution methods have been developed as an alternative. However, existing methods are designed for tissues consisting of clearly distinguishable cell types and have difficulties estimating highly correlated or rare cell types. To address this challenge, we propose hierarchical deconvolution (HiDecon) that uses single-cell RNA sequencing references and a hierarchical cell-type tree, which models the similarities among cell types and cell differentiation relationships, to estimate cellular fractions in bulk data. By coordinating cell fractions across layers of the hierarchical tree, cellular fraction information is passed up and down the tree, which helps correct estimation biases by pooling information across related cell types. The flexible hierarchical tree structure also enables estimating rare cell fractions by splitting the tree to higher resolutions. Through simulations and real data applications with the ground truth of measured cellular fractions, we demonstrate that HiDecon outperforms existing methods and accurately estimates cellular fractions. Finally, we show the utility of HiDecon estimates in identifying the associations between cellular fractions and Alzheimer’s disease.

Multivariate longitudinal data are frequently encountered in practice such as in our motivating longitudinal microbiome study. It is of general interest to associate such high-dimensional, longitudinal measures with some univariate continuous outcome. However, incomplete observations are common in a regular study design, as not all samples are measured at every time point, giving rise to the so-called blockwise missing values. Such missing structure imposes significant challenges for association analysis and defies many existing methods that require complete samples. In this paper we propose to represent multivariate longitudinal data as a three-way tensor array (i.e., sample-by-feature-by-time) and exploit a parsimonious scalar-on-tensor regression model for association analysis. We develop a regularized covariance-based estimation procedure that effectively leverages all available observations without imputation. The method achieves variable selection and smooth estimation of time-varying effects. The application to the motivating microbiome study reveals interesting links between the preterm infant’s gut microbiome dynamics and their neurodevelopment. Additional numerical studies on synthetic data and a longitudinal aging study further demonstrate the efficacy of the proposed method.

Let a collection of networks represent interactions within several (social or ecological) systems. We pursue two objectives: identifying similarities in the topological structures that are held in common between the networks and clustering the collection into subcollections of structurally homogeneous networks. We tackle these two questions with a probabilistic model-based approach. We propose an extension of the stochastic block model (SBM) adapted to the joint modeling of a collection of networks. The networks in the collection are assumed to be independent realizations of SBMs. The common connectivity structure is imposed through the equality of some parameters.

The model parameters are estimated with a variational expectation-maximization (EM) algorithm. We derive an ad hoc penalized likelihood criterion to select the number of blocks and to assess the adequacy of the consensus found between the structures of the different networks. This same criterion can also be used to cluster networks on the basis of their connectivity structure. It thus provides a partition of the collection into subsets of structurally homogeneous networks.

The relevance of our proposition is assessed on two collections of ecological networks. First, an application to three stream food webs reveals the homogeneity of their structures and the correspondence between groups of species in different ecosystems playing equivalent ecological roles. Moreover, the joint analysis allows a finer analysis of the structure of smaller networks. Second, we cluster 67 food webs according to their connectivity structures and demonstrate that five mesoscale structures are sufficient to describe this collection.

Sports organizations often want to estimate athlete strengths. For games with scored outcomes, a common approach is to assume observed game scores follow a normal distribution conditional on athletes’ latent abilities, which may change over time. In many games, however, this assumption of conditional normality does not hold. To estimate athletes’ time-varying latent abilities using nonnormal game score data, we propose a Bayesian dynamic linear model with flexible monotone response transformations. Our model learns nonlinear monotone transformations to address nonnormality in athlete scores and can be easily fit using standard regression and optimization routines, which we implement in the $\mathtt{dlmt}$ package in $\mathtt{R}$. We demonstrate our method on data from several Olympic sports, including biathlon, diving, rugby, and fencing.

Many important policy decisions concerning policing hinge on our understanding of how likely various criminal offenses are to result in arrests. Since many crimes are never reported to law enforcement, estimates based on police records alone must be adjusted to account for the likelihood that each crime would have been reported to the police. In this paper we present a methodological framework for estimating the likelihood of arrest from police data that incorporates estimates of crime reporting rates computed from a victimization survey. We propose a parametric regression-based two-step estimator that: (i) estimates the likelihood of crime reporting using logistic regression with survey weights and then (ii) applies a second regression step to model the likelihood of arrest. Our empirical analysis focuses on racial disparities in arrests for violent crimes (sex offenses, robbery, aggravated and simple assaults) from 2006–2015 police records from the National Incident Based Reporting System (NIBRS), with estimates of crime reporting obtained using 2003–2020 data from the National Crime Victimization Survey (NCVS). We find that, after adjusting for unreported crimes, the likelihood of arrest computed from police records decreases significantly. We also find that, while incidents with white offenders, on average, result in arrests more often than those with black offenders, the disparities tend to be small after accounting for crime characteristics and unreported crimes.

The adrenocorticotropic hormone and cortisol play critical roles in stress regulation and the sleep-wake cycle. Most research has been focused on how the two hormones regulate each other in terms of short-term pulses. Few studies have been conducted on the circadian relationship between the two hormones and how it differs between normal and abnormal groups. The circadian patterns are difficult to model as parametric functions. Directly extending univariate functional mixed effects models would result in a large dimensional problem and a challenging nonparametric inference. In this article we propose a semiparametric bivariate hierarchical state space model in which each hormone profile is modeled by a hierarchical state space model with nonparametric population-average and subject-specific components. The bivariate relationship is constructed by concatenating two latent independent subject-specific random functions specified by a design matrix, leading to a parametric inference on the correlation. We propose a computationally efficient state-space EM algorithm for estimation and inference. We apply the proposed method to a study of chronic fatigue syndrome and fibromyalgia and discover an erratic regulation pattern in the patient group in contrast to a circadian regulation pattern conforming to the day–night cycle in the control group.

Neuroimaging studies often involve predicting a scalar outcome from an array of images collectively called tensor. The use of magnetic resonance imaging (MRI) provides a unique opportunity to investigate the structures of the brain. To learn the association between MRI images and human intelligence, we formulate a scalar-on-image quantile regression framework. However, the high dimensionality of the tensor makes estimating the coefficients for all elements computationally challenging. To address this, we propose a low-rank coefficient array estimation algorithm, based on tensor train (TT) decomposition, which we demonstrate can effectively reduce the dimensionality of the coefficient tensor to a feasible level while ensuring adequacy to the data. Our method is more stable and efficient compared to the commonly used canonic polyadic rank approximation-based method. We also propose a generalized lasso penalty on the coefficient tensor to take advantage of the spatial structure of the tensor, further reduce the dimensionality of the coefficient tensor, and improve the interpretability of the model. The consistency and asymptotic normality of the TT estimator are established under some mild conditions on the covariates and random errors in quantile regression models. The rate of convergence is obtained with regularization under the total variation penalty. Extensive numerical studies, including both synthetic and real MRI imaging data, are conducted to examine the empirical performance of the proposed method and its competitors.

Schizophrenia is a severe mental disorder that distorts patients’ perception of reality, and its treatment with antipsychotics can lead to significant side effects. Despite the heterogeneity in patient responses to treatments, most existing studies on individualized treatment regimes only focus on optimizing treatment efficacy, disregarding potential negative effects. To fill this gap, we propose a restricted outcome weighted learning method that optimizes efficacy outcomes while adhering to individual-level negative effect constraints. Our method is developed for multistage treatment decision problems that include single-stage decision as a special case. We propose an efficient learning algorithm that utilizes the difference-of-convex algorithm and the Lagrange multiplier to solve nonconvex optimization with nonconvex risk constraints. We also establish theoretical properties, including Fisher consistency and strong duality results, for the proposed method. We apply our method to a clinical study to design effective schizophrenia treatment [Stroup et al. (Schizophr. Bull. 29 (2003) 15–31)] and find that our approach reduces side-effect risk by at least 22.5% and improves efficacy by at least 26.3% compared to competing methods. In addition, we discover that certain covariates, such as the PANSS score, clinician global impressions severity score, and BMI, have a significant impact on controlling side effects and determining optimal treatment recommendations. These results are valuable in identifying subgroups of patients who need special attention when prescribing more aggressive treatment plans.

Accurate hospital performance measurement is important to both patients and providers but is challenging due to case-mix heterogeneity, differences in treatment guidelines, and data privacy regulations that preclude the sharing of individual patient data. Motivated to overcome these issues in the setting of hospital quality measurement, we develop a federated causal inference framework. We devise a doubly robust estimator of the mean potential outcome in a target population and show that it is consistent even when some models are misspecified. To enable real-world use, our proposed algorithm is privacy-preserving (requiring only summary statistics to be shared between hospitals) and communication-efficient (requiring only one round of communication between hospitals). We show that our estimator has good finite sample properties in simulation studies. We investigate the quality of hospital care provided by a diverse set of 51 candidate Cardiac Centers of Excellence, as measured by 30-day mortality and length of stay for acute myocardial infarction (AMI) patients. We find that our proposed federated global estimator improves the precision of treatment effect estimates by 34% to 86%, compared to using data from the target hospital alone. This precision gain results in qualitatively different conclusions about the estimated effect of percutaneous coronary intervention (PCI), compared to medical management (MM) in 43% (22 of 51) of hospitals.

Environmental exposures, such as cigarette smoking, influence health outcomes through intermediate molecular phenotypes, such as the methylome, transcriptome, and metabolome. Mediation analysis is a useful tool for investigating the role of potentially high-dimensional intermediate phenotypes in the relationship between environmental exposures and health outcomes. However, little work has been done on mediation analysis when the mediators are high-dimensional and the outcome is a survival endpoint, and none of it has provided a robust measure of total mediation effect. To this end, we propose an estimation procedure for Mediation Analysis of Survival outcome and High-dimensional omics mediators (MASH), based on a second-moment-based measure of total mediation effect for survival data analogous to the ${\mathit{R}}^2$ measure in a linear model. In addition, we propose a three-step mediator selection procedure to mitigate potential bias induced by nonmediators. Extensive simulations showed good performance of MASH in estimating the total mediation effect and identifying true mediators. By applying MASH to the metabolomics data of 1919 subjects in the Framingham Heart Study, we identified five metabolites as mediators of the effect of cigarette smoking on coronary heart disease risk (total mediation effect, 51.1%) and two metabolites as mediators between smoking and risk of cancer (total mediation effect, 50.7%). Application of MASH to a diffuse large B-cell lymphoma genomics data set identified copy-number variations for eight genes as mediators between the baseline International Prognostic Index score and overall survival.

We develop flexible multivariate spatiotemporal Hawkes process models to analyze patterns of terrorism. Previous applications of point process methods to political violence data mainly utilize temporal Hawkes process models, neglecting spatial variation in these attack patterns. This limits what can be learned from these models, as any effective counter-terrorism strategy requires knowledge on both when and where attacks are likely to occur. Even the existing work on spatiotemporal Hawkes processes imposes restrictions on the triggering function that are not well-suited for terrorism data. Therefore, we generalize the structure of the spatiotemporal triggering function considerably, allowing for nonseparability, nonstationarity, and cross-triggering (across multiple terror groups). To demonstrate the utility of our models, we analyze two samples of real-world terrorism data: Afghanistan (2002–2013) as a univariate analysis and Nigeria (2009–2017) as a bivariate analysis. Jointly, these two studies demonstrate that our generalized models outperform standard Hawkes process models, besting widely-used alternatives in overall model fit and revealing spatiotemporal patterns that are, by construction, masked in these models (e.g., increasing dispersion in cross-triggering over time).

Extreme weather events associated with climate change have caused significant damages. In particular, hail storms damage millions of properties in the U.S. and result in billion-dollar insured losses each year in the recent decade. To facilitate the insurance claims management operations in insurance companies, we construct a hierarchical dependence model, which accommodates the complex dependence within and between the outcomes of interests including the propensity of filing a claim, time to report a claim, and the claim amount. The storm-specific and property-specific characteristics are incorporated through marginal models, such as generalized linear models and survival analysis models. The dependence within the hail event is captured by spatial factor copula, while the dependence between different outcomes is captured by bivariate copula. For parameter estimation we develop a two-step procedure that first maximizes the marginal likelihood function and then maximizes the pairwise likelihood, which ensures computational feasibility for big data. We apply this modeling framework to analyze a large dataset involving hail storms in Colorado from 2011 to 2015 impacting hundreds of thousands of insured properties and demonstrate that the predictive performance can be improved by our proposed methodology.

The relationship between inflation and predictors, such as unemployment, is potentially nonlinear with a strength that varies over time, and prediction errors may be subject to large, asymmetric shocks. Inspired by these concerns, we develop a model for inflation forecasting that is nonparametric both in the conditional mean and in the error using Gaussian and Dirichlet processes, respectively. We discuss how both these features may be important in producing accurate forecasts of inflation. In a forecasting exercise involving CPI inflation, we find that our approach has substantial benefits, both overall and in the left tail, with nonparametric modeling of the conditional mean being of particular importance.

We propose and investigate a hidden Markov model (HMM) for the analysis of dependent, aggregated, superimposed two-state signal recordings. A major motivation for this work is that often these signals cannot be observed individually but only their superposition. Among others, such models are in high demand for the understanding of cross-talk between ion channels, where each single channel cannot be measured separately. As an essential building block, we introduce a parameterized vector norm dependent Markov chain model and characterize it in terms of permutation invariance as well as conditional independence. This building block leads to a hidden Markov chain sum process which can be used for analyzing the dependence structure of superimposed two-state signal observations within an HMM. Notably, the model parameters of the vector norm dependent Markov chain are uniquely determined by the parameters of the sum process and are, therefore, identifiable. We provide algorithms to estimate the parameters, discuss model selection and apply our methodology to real-world ion channel data from the heart muscle, where we show competitive gating.

Methods utilizing instrumental variables have been a fundamental statistical approach to causal estimation in the presence of unmeasured confounding, usually occurring in nonrandomized observational data common to fields such as economics and public health. However, such methods traditionally make constricting linearity and additivity assumptions that are inapplicable to the complex modeling challenges of today. The growing body of observational data being collected may benefit from flexible regression modeling while also retaining the ability to control for confounding using instrumental variables. Therefore, this article presents a flexible instrumental variable regression model based on Bayesian regression tree ensembles to estimate the causal exposure-outcome relationship, including interactions with covariates, in the presence of confounding. One exciting application of this method is to use genetic variants as instruments, known as Mendelian randomization. We present our flexible Bayesian instrumental variable regression tree method with an example from the UK Biobank where body mass index is related to blood pressure using genetic variants as the instruments. Body mass index is one factor that is hypothesized to have a nonlinear relationship with cardiovascular risk factors, such as blood pressure, while interacting with age. Heterogeneity in patient characteristics, such as age, could be clinically interesting from a precision medicine perspective where individualized treatment is emphasized.

High-dimensional biomarkers, such as gene expression profiles, are often collected longitudinally to monitor disease progression in clinical studies, where the primary endpoint of interest is often a survival outcome. It is of great interest to study the associations between high-dimensional longitudinal biomarkers and the survival outcome as well as to identify biomarkers related to the survival outcome. Joint models, which have been extensively studied in the past decades, are commonly used to study the associations between longitudinal biomarkers and the survival outcome. However, existing joint models only consider one or a few longitudinal biomarkers and cannot deal with high-dimensional longitudinal biomarkers. In this paper we propose a novel penalized joint model that can handle high-dimensional longitudinal biomarkers. Specifically, we impose an adaptive lasso penalty on the parameters for the effects of the longitudinal biomarkers on the survival outcome, which allows for variable selection. We also develop a computationally efficient algorithm for model estimation based on the Gaussian variational approximation method, which can be implemented using the HDJM package in R. Furthermore, based on the penalized joint model, we propose a two-stage selection procedure that can reduce the estimation bias, due to the penalization, and allows for inference. We conduct extensive simulation studies to evaluate the performance of our proposed method. The performance of our proposed method is further demonstrated on a longitudinal gene expression dataset of patients with idiopathic pulmonary fibrosis.

Test-negative designs have rapidly become an appealing approach to assess disease interventions when randomization is not feasible and specifically used to measure the effectiveness of vaccines in the field (Vaccine 31 (2013) 2165–2168). An innovative extension of the test-negative design was recently used to assess the impact of a mosquito intervention where the intervention was applied at a cluster level with cluster assignment chosen at random, the AWED (applying Wolbachia to eliminate dengue) trial. The primary analysis reported was intention-to-treat (ITT) (Trials 19 (2018) 302; N. Engl. J. Med. 384 (2021) 2177–2186). However, the level of uptake of the intervention on mosquitoes was routinely captured in all clusters over time, and, furthermore, participants’ mobility across clusters was measured in the time immediately preceding the onset of symptoms (whether test-positive or test-negative). Combinations of these measurements provide proxies for the true exposure to the intervention, thereby permitting an “as treated” assessment. We consider the use of marginal generalized estimating equations (GEE) and conditional generalized inear mixed models (GLMM) to estimate as treated efficacy, contrasting both with the ITT. We illustrate the strengths and challenges of these methods in the context of the AWED trial, highlighting several ways that common approaches to analysis of clustered data can yield incorrect results that can in turn be obscured and compounded by limitations in routine software. In addition, we estimate a greater level of intervention efficacy than shown in the ITT analysis.

Quantifying changes in the probability and magnitude of extreme flooding events is key to mitigating their impacts. While hydrodynamic data are inherently spatially dependent, traditional spatial models, such as Gaussian processes, are poorly suited for modeling extreme events. Spatial extreme value models with more realistic tail dependence characteristics are under active development. They are theoretically justified but give intractable likelihoods, making computation challenging for small datasets and prohibitive for continental-scale studies. We propose a process mixture model (PMM) which specifies spatial dependence in extreme values as a convex combination of a Gaussian process and a max-stable process, yielding desirable tail dependence properties but intractable likelihoods. To address this, we employ a unique computational strategy where a feed-forward neural network is embedded in a density regression model to approximate the conditional distribution at one spatial location, given a set of neighbors. We then use this univariate density function to approximate the joint likelihood for all locations by way of a Vecchia approximation. The PMM is used to analyze changes in annual maximum streamflow within the U.S. over the last 50 years and is able to detect areas which show increases in extreme streamflow over time.

Cancer screening facilitates the early detection of cancer at a stage when treatment is often most effective. However, it also brings the risk of overdiagnosis, where a diagnosis made through screening would not have led to symptoms or death during the patient’s lifetime. In this paper we tackle a significant unresolved issue in the evaluation of screening efficacy: selecting primary endpoints and inferential procedures that efficiently consider potential overdiagnosis in screening trials. This is motivated by the necessity to design and analyze a phase IV Early Detection Initiative (EDI) trial for evaluating a pancreatic cancer screening strategy. We introduce two novel approaches for assessing screening efficacy, grounded on cancer stage shift. These methods address potential overdiagnosis by: (i) borrowing information about clinical diagnosis from the control arm that hasn’t undergone screening (the BR approach) and (ii) performing sensitivity analysis, contingent upon a conservative bound of the overdiagnosis magnitude (the SEN-T approach). Analytical methods and extensive simulation studies underscore the superiority of our proposed methods, demonstrating enhanced efficiency in estimating and testing screening efficacy compared to existing methods. The latter either overlook overdiagnosis or adhere to a valid, yet conservative, cumulative incidence endpoint. We illustrate the practical application of these approaches using ovarian cancer data from the Prostate, Lung, Colorectal, and Ovarian Cancer Screening Trial. The results affirm that our methods bolster an efficient and robust study design for cancer screening trials.

Small area population counts are necessary for many epidemiological studies, yet their quality and accuracy are often not assessed. In the United States, small area population counts are published by the United States Census Bureau (USCB) in the form of the decennial census counts, intercensal population projections (PEP), and American Community Survey (ACS) estimates. Although there are significant relationships between these three data sources, there are important contrasts in data collection, data availability, and processing methodologies such that each set of reported population counts may be subject to different sources and magnitudes of error. Additionally, these data sources do not report identical small area population counts due to post-survey adjustments specific to each data source. Consequently, in public health studies, small area disease/mortality rates may differ depending on which data source is used for denominator data. To accurately estimate annual small area population counts and their associated uncertainties, we present a Bayesian population (BPop) model, which fuses information from all three USCB sources, accounting for data source specific methodologies and associated errors. We produce comprehensive small area race-stratified estimates of the true population, and associated uncertainties, given the observed trends in all three USCB population estimates. The main features of our framework are: (1) a single model integrating multiple data sources, (2) accounting for data source specific data generating mechanisms and specifically accounting for data source specific errors, and (3) prediction of population counts for years without USCB reported data. We focus our study on the Black and White only populations for 159 counties of Georgia and produce estimates for years 2006–2023. We compare BPop population estimates to decennial census counts, PEP annual counts, and ACS multi-year estimates. Additionally, we illustrate and explain the different types of data source specific errors. Lastly, we compare model performance using simulations and validation exercises. Our Bayesian population model can be extended to other applications at smaller spatial granularity and for demographic subpopulations defined further by race, age, and sex, and/or for other geographical regions.

In this paper we predict sea surface salinity (SSS) in the Arctic Ocean based on satellite measurements. SSS is a crucial indicator for ongoing changes in the Arctic Ocean and can offer important insights about climate change. We particularly focus on areas of water mistakenly flagged as ice by satellite algorithms. To remove bias in the retrieval of salinity near sea ice, the algorithms use conservative ice masks, which result in considerable loss of data. We aim to produce realistic SSS values for such regions to obtain more complete understanding about the SSS surface over the Arctic Ocean and benefit future applications that may require SSS measurements near edges of sea ice or coasts. We propose a class of scalable nonstationary processes that can handle large data from satellite products and complex geometries of the Arctic Ocean. Barrier overlap-removal acyclic directed graph GP (BORA-GP) constructs sparse directed acyclic graphs (DAGs) with neighbors conforming to barriers and boundaries, enabling characterization of dependence in constrained domains. The BORA-GP models produce more sensible SSS values in regions without satellite measurements and show improved performance in various constrained domains in simulation studies compared to state-of-the-art alternatives. An $\mathtt{R}$ package is available at https://github.com/jinbora0720/boraGP.

Standard models for longitudinal data ignore the stochastic nature of time-varying covariates and their stochastic evolution over time by treating them as fixed variables. There have been recent methods for modelling time-varying covariates; however, those methods cannot be applied to analyse longitudinal data when the longitudinal response and the time-varying covariates for each subject are measured at different time points. Moreover, it is difficult to study the temporal effects of a time-varying covariate on the longitudinal response and the temporal correlation between them. Motivated by data from an AIDS cohort study conducted over 26 years at the University Hospitals Leuven in which the measurements on the CD4 cell count and viral load for patients are not taken at the same time point, we present a framework to address those challenges by using joint multivariate mixed models to jointly model time-varying covariates and a longitudinal response, instead of including time-varying covariates in the response model. This approach also has the advantage that one can study the association between the covariate at any time point and the response at any other time point without having to explicitly model the conditional distribution of the response given the covariate. We use penalised spline functions of time to capture the evolutions of both the response and time-varying covariates over time.

Longitudinal biomarker data and cross-sectional outcomes are routinely collected in modern epidemiology studies, often with the goal of informing tailored early intervention decisions. For example, hormones, such as estradiol (E2) and follicle-stimulating hormone (FSH), may predict changes in womens’ health during the midlife. Most existing methods focus on constructing predictors from mean marker trajectories. However, subject-level biomarker variability may also provide critical information about disease risks and health outcomes. Current literature does not provide statistical models to investigate such relationships with valid uncertainty quantification. In this paper we develop a fully Bayesian joint model that estimates subject-level means, variances, and covariances of multiple longitudinal biomarkers and uses these as predictors to evaluate their respective associations with a cross-sectional health outcome. Simulations demonstrate excellent recovery of true model parameters. The proposed method provides less biased and more efficient estimates, relative to alternative approaches that either ignore subject-level differences in variances or perform two-stage estimation where estimated marker variances are treated as observed. Empowered by the model, analyses of women’s health data reveal, for the first time, that larger variability of E2 was associated with slower increases in waist circumference across the menopausal transition.

Multivariate functional data that are cross-sectionally compositional data are attracting increasing interest in the statistical modeling literature, a major example being trajectories over time of compositions derived from cause-specific mortality rates. In this work we develop a novel functional concurrent regression model in which independent variables are functional compositions. This allows us to investigate the relationship over time between life expectancy at birth and compositions derived from cause-specific mortality rates of four distinct age classes, namely, zero to four, five to 39, 40–64 and 65+ in 25 countries. A penalized approach is developed to estimate the regression coefficients and select the relevant variables. Then an efficient computational strategy, based on an augmented Lagrangian algorithm, is derived to solve the resulting optimization problem. The good performances of the model in predicting the response function and estimating the unknown functional coefficients are shown in a simulation study. The results on real data confirm the important role of neoplasms and cardiovascular diseases in determining life expectancy emerged in other studies and reveal several other contributions not yet observed.

This paper studies the dynamic patterns of the prelaunch online movie reviews, or movie electronic word-of-mouth (eWOM), over time and investigates their relations to the subsequent box office revenues. The volume and valence of prelaunch eWOM have been shown to be early indicators of strong or weak box office. The time patterns of prelaunch eWOM evolution, which are essentially functional data, on the other hand, tend to be overlooked. We apply the functional principal component analysis, a dimension reduction technique in functional data analysis, to analyze the dynamic patterns of various quantile trajectories of the movie eWOM, instead of directly studying the whole eWOM functional data. The functional principal component (FPC) scores of quantile trajectories at various quantile levels are used to predict the box office revenues. We use the sparse group lasso method to select the quantile levels and individual FPC scores that make significant contributions to the prediction of box office revenues. The results show that compared with other measures, such as valence and variance, the top-end quantiles would be a better measure in capturing the relations between the prelaunch product ratings time pattern and launch sales.

Weather observations are important for a wide range of applications although they do pose statistical challenges, such as missing values, errors, flawed outliers and poor spatial and temporal coverage to name a few. A Bayesian hierarchical spline framework is presented here to deal with such challenges in temperature time series. Motivated by a real-life problem, the approach uses penalised splines, constructed hierarchically, to pool the data, along with a discrete mixture distribution to deal with outliers and publicly available global reanalysis data sets (climate model data) to integrate physically constrained information. Efficient Bayesian implementation is achieved using conditional conjugacy, which allows thorough model checking and uncertainty quantification. Fitting the model to daily maximum temperature illustrates its flexibility in capturing temporal structures, in pooling of the information and in outlier detection. The model is used to hindcast the time series 50 years into the past while maintaining uncertainty at reasonable levels.

Mendelian randomization (MR) is a widely-used method to estimate the causal relationship between a risk factor and disease. A fundamental part of any MR analysis is to choose appropriate genetic variants as instrumental variables. Genome-wide association studies often reveal that hundreds of genetic variants may be robustly associated with a risk factor, but in some situations investigators may have greater confidence in the instrument validity of only a smaller subset of variants. Nevertheless, the use of additional instruments may be optimal from the perspective of mean squared error, even if they are slightly invalid; a small bias in estimation may be a price worth paying for a larger reduction in variance. For this purpose we consider a method for “focused” instrument selection whereby genetic variants are selected to minimise the estimated asymptotic mean squared error of causal effect estimates. In a setting of many weak and locally invalid instruments, we propose a novel strategy to construct confidence intervals for postselection focused estimators that guards against the worst case loss in asymptotic coverage. In empirical applications to: (i) validate lipid drug targets and (ii) investigate vitamin D effects on a wide range of outcomes, our findings suggest that the optimal selection of instruments does not involve only a small number of biologically-justified instruments but also many potentially invalid instruments.

Investigating technical skills of swimmers is a challenge for performance improvement that can be achieved by analyzing multivariate functional data recorded by inertial measurement units (IMU). To investigate technical levels of front-crawl swimmers, a new model-based approach is introduced to obtain two complementary partitions reflecting, for each swimmer, its swimming pattern and its ability to reproduce it. Contrary to the usual approaches for functional data clustering, the proposed approach also considers the information of the error terms resulting from the functional basis decomposition. Indeed, after decomposing into functional basis with finite number of elements both the original signal (measuring the swimming pattern) and the signal of squared error terms (measuring the ability to reproduce the swimming pattern), the method fits the joint distribution of the coefficients related to both decompositions by considering dependency between both partitions. Modeling this dependency is mandatory since the difficulty of reproducing a swimming pattern depends on its shape. Moreover, a sparse decomposition of the distribution within components that permits a selection of the relevant dimensions during clustering is proposed. The partitions obtained on the IMU data aggregate the kinematical stroke variability linked to swimming technical skills and allow relevant biomechanical strategy for front-crawl sprint performance to be identified.