Phylodynamics is a set of population genetics tools that aim at reconstructing demographic history of a population based on molecular sequences of individuals sampled from the population of interest. One important task in phylodynamics is to estimate changes in (effective) population size. When applied to infectious disease sequences, such estimation of population size trajectories can provide information about changes in the number of infections. To model changes in the number of infected individuals, current phylodynamic methods use nonparametric approaches (e.g., Bayesian curve-fitting based on change-point models or Gaussian process priors), parametric approaches (e.g., based on differential equations), and stochastic modeling in conjunction with likelihood-free Bayesian methods. The first class of methods yields results that are hard to interpret epidemiologically. The second class of methods provides estimates of important epidemiological parameters, such as infection and removal/recovery rates, but ignores variation in the dynamics of infectious disease spread. The third class of methods is the most advantageous statistically but relies on computationally intensive particle filtering techniques that limits its applications. We propose a Bayesian model that combines phylodynamic inference and stochastic epidemic models and achieves computational tractability by using a linear noise approximation (LNA)—a technique that allows us to approximate probability densities of stochastic epidemic model trajectories. LNA opens the door for using modern Markov chain Monte Carlo tools to approximate the joint posterior distribution of the disease transmission parameters and of high dimensional vectors describing unobserved changes in the stochastic epidemic model compartment sizes (e.g., numbers of infectious and susceptible individuals). In a simulation study we show that our method can successfully recover parameters of stochastic epidemic models. We apply our estimation technique to Ebola genealogies estimated using viral genetic data from the 2014 epidemic in Sierra Leone and Liberia.

Social distancing is widely acknowledged as an effective public health policy combating the novel coronavirus. But extreme forms of social distancing, like isolation and quarantine, have costs, and it is not clear how much social distancing is needed to achieve public health effects. In this article we develop a design-based framework to test the causal null hypothesis and make inference about the dose-response relationship between reduction in social mobility and COVID-19 related public health outcomes. We first discuss how to embed observational data with a time-independent, continuous treatment dose into an approximate randomized experiment and develop a randomization-based procedure that tests if a structured dose-response relationship fits the data. We then generalize the design and testing procedure to a longitudinal setting and apply them to investigate the effect of social distancing during the first phased reopening in the United States on public health outcomes using data compiled from $\mathsf{Unacast}{}^{\text{™}}$, the United States Census Bureau, and the County Health Rankings and Roadmaps Program. We rejected a primary analysis null hypothesis that stated the social distancing from April 27, 2020 to June 28, 2020, had no effect on the COVID-19-related death toll from June 29, 2020 to August 2, 2020 (p-value < 0.001), and found that it took more reduction in mobility to prevent exponential growth in case numbers for nonrural counties compared to rural counties.

Air pollution is a major global public health risk factor. Among all air pollutants, PM${}_{2.5}$ is especially harmful. It has been well demonstrated that chronic exposure to PM${}_{2.5}$ can cause many health problems, including asthma, lung cancer and cardiovascular diseases. To tackle problems caused by air pollution, governments have put a huge amount of resources to improve air quality and reduce the impact of air pollution on public health. In this effort it is extremely important to develop an air pollution surveillance system to constantly monitor the air quality over time and to give a signal promptly once the air quality is found to deteriorate so that a timely government intervention can be implemented. To monitor a sequential process, a major statistical tool is the statistical process control (SPC) chart. However, traditional SPC charts are based on the assumptions that process observations at different time points are independent and identically distributed. These assumptions are rarely valid in environmental data because seasonality and serial correlation are common in such data. To overcome this difficulty, we suggest a new control chart in this paper, which can properly accommodate dynamic temporal pattern and serial correlation in a sequential process. Thus, it can be used for effective air pollution surveillance. This method is demonstrated by an application to monitor the daily average PM${}_{2.5}$ levels in Beijing and shown to be effective and reliable in detecting the increase of PM${}_{2.5}$ levels.

Regression discontinuity designs are extensively used for causal inference in observational studies. However, they are usually confined to settings with simple treatment rules and determined by a single running variable with a single cutoff. Motivated by the problem of estimating the impact of grade retention on educational and juvenile crime outcomes in Chile, we propose a framework and methods for complex discontinuity designs that encompass multiple treatment rules. In this framework the observed covariates play a central role for identification, estimation, and generalization of causal effects. Identification is nonparametric and relies on a local strong ignorability assumption. Estimation proceeds, as in any observational study, under strong ignorability, yet in a neighborhood of the cutoffs of the running variables. We discuss estimation approaches based on matching and weighting, including complementary regression modeling adjustments. We present assumptions for generalization, that is, for identification and estimation of average treatment effects for target populations. We also describe two approaches to select the neighborhood for analysis. We find that grade retention in Chile has a negative impact on future grade retention but is not associated with dropping out of school or committing a juvenile crime.

In the National Longitudinal Study of Adolescent to Adult Health (Add Health) in the United States during the 1994–95 school year, sampled high school students were asked to rank 17 romantic relationship activities in ideal chronological order and to nominate up to five male and five female best friends within the school. It is natural to see that the students’ rank-order progression of social, romantic and sexual relationships may be influenced strongly by their peers or friends. So far, traditional ranking models do not account for such social network dependency. In this article we introduce a new class of models, called social order statistics (SOS) models, to learn ranking data in social networks. The new models combine the order statistics models and spatial autoregressive models to account for social dependencies among the individuals. A flexible formulation of weight matrices in the spatial model is adopted to provide diverse network effects among the individuals for different items. Efficient and scalable MCMC algorithms are developed to perform Bayesian inference in a parallel manner for large networks with even a few thousand nodes. Analysis of the Add Health dataset reveals that social network effects are different for students’ preferences toward different activities in a relationship. In particular, students’ preferences on romantic and sexual events generally have stronger peer effects than those on social events, and opinions of close friends of the same gender tend to have a larger impact on the preferences on sexual events.

Accurate HIV incidence estimation, based on individual recent infection status (recent vs long-term infection), is important for monitoring the epidemic, targeting interventions to those at greatest risk of new infection, and evaluating existing programs of prevention and treatment. Starting from 2015, the population-based HIV impact assessment (PHIA) individual-level surveys are implemented in the most-affected countries in sub-Saharan Africa. PHIA is a nationally-representative HIV-focused survey that combines household visits with key questions and cutting-edge technologies, such as biomarker tests for HIV antibody and HIV viral load which offer the unique opportunity of distinguishing between recent infection and long-term infection, and providing relevant HIV information by age, gender, and location. In this article we propose a semisupervised logistic regression model for estimating individual level HIV recency status. It incorporates information from multiple data sources—the PHIA survey, where the true HIV recency status is unknown, and the cohort studies provided in the literature where the relationship between HIV recency status and the covariates are presented in the form of a contingency table. It also utilizes the national level HIV incidence estimates from the epidemiology model. Applying the proposed model to Malawi PHIA data, we demonstrate that our approach is more accurate for the individual level estimation and more appropriate for estimating HIV recency rates at aggregated levels than the current practice—the binary classification tree (BCT).

Analyzing the escape direction of animals subject to covariates is a problem that requires statistical techniques beyond classical regression methods. Apart from the periodicity of the angle of direction, which demands the use of circular statistics, animal escape data usually call for the exploration of the preferred orientations rather than the expected orientation. In this paper we propose the use of a nonparametric method to estimate the conditional local modes of the escape directions of animals from a regression perspective. We present the estimation algorithms and study the asymptotic properties of the estimators as well as its finite sample performance through some simulation experiments. Our proposal is used to model the escape behavior of a group of larval zebrafish escaping from a robot predator. More broadly, the approach presented in this paper can be applied to many existing problems related to animal behavior or other fields.

Many factors must be taken into account when designing an observational study. Unlike controlled studies, observational studies cannot mitigate the effects of confounding through randomization, and such factors should be incorporated into both the study analysis and the study design. Unfortunately, there is often little data available on most of these factors at the design stage, rendering it infeasible to reliably postulate the impact of these factors on the treatment effect estimate and precision of such an effect. In this work we demonstrate how failure to account for adjustment covariates in the design stage of observational studies utilizing group sequential designs has deleterious effects on the study’s observed power. We do this by constructing a hypothetical study of Alzheimer’s disease biomarkers on cognition, using data from the Alzheimer’s Disease Neuroimaging Initiative. We then outline a constrained boundaries procedure that uses data collected at interim analyses to better estimate variance and correct stopping boundaries to maintain power and type I error, while allowing for early study termination.

Epilepsy is a chronic neurological disorder; it affects more than 50 million people globally. An epileptic seizure acts like a temporary shock to the neuronal system, disrupting normal electrical activity in the brain. Epilepsy is frequently diagnosed with electroencephalograms (EEGs). Current methods study only the time-varying spectra and coherence but do not directly model changes in extreme behavior, neglecting the fact that neuronal oscillations exhibit non-Gaussian heavy-tailed probability distributions. To overcome this limitation, we propose a new approach to characterize brain connectivity based on the joint tail (i.e., extreme) behavior of the EEGs. Our proposed method, the conditional extremal dependence for brain connectivity (Conex–Connect), is a pioneering approach that links the association between extreme values of higher oscillations at a reference channel with the other brain network channels. Using the Conex–Connect method, we discover changes in the extremal dependence driven by the activity at the foci of the epileptic seizure. Our model-based approach reveals that, preseizure, the dependence is notably stable for all channels when conditioning on extreme values of the focal seizure area. By contrast, the dependence between channels is weaker during the seizure, and dependence patterns are more “chaotic.” Using the Conex–Connect method, we identified the high-frequency oscillations as the most relevant features, explaining the conditional extremal dependence of brain connectivity.

While the study of a single network is well established, technological advances now allow for the collection of multiple networks with relative ease. Increasingly, anywhere from several to thousands of networks can be created from brain imaging, gene coexpression data, or microbiome measurements. And these networks, in turn, are being looked to as potentially powerful features to be used in modeling. However, with networks being non-Euclidean in nature, how best to incorporate them into standard modeling tasks is not obvious. In this paper we propose a Bayesian modeling framework that provides a unified approach to binary classification, anomaly detection, and survival analysis with network inputs. We encode the networks in the kernel of a Gaussian process prior via their pairwise differences, and we discuss several choices of provably positive definite kernel that can be plugged into our models. Although our methods are widely applicable, we are motivated here, in particular, by microbiome research (where network analysis is emerging as the standard approach for capturing the interconnectedness of microbial taxa across both time and space) and its potential for reducing preterm delivery and improving personalization of prenatal care.

Detection and attribution analyses play a central role in establishing the causal effect of human activities on global warming. The most commonly used method in such analyses, optimal fingerprinting, is a multiple regression where each covariate has a measurement error whose covariance matrix is the same as that of the regression error up to a known scale. Inferences about the regression coefficients are critical not only for making statements about detection and attribution but also for quantifying the uncertainty in important outcomes derived from detection and attribution analyses. When there are no errors-in-variables (EIV), the optimal weight matrix in estimating the regression coefficients is the precision matrix of the regression error. This matrix, however, is never known and has to be estimated from climate model simulations with regularization. The consequence is that the optimal fingerprinting method is not optimal, as believed in practice. We construct a weight matrix by inverting a nonlinear shrinkage estimate of the error covariance matrix that minimizes loss functions directly targeting the uncertainty of the resulting regression coefficient estimator. The resulting estimator of the regression coefficients is asymptotically optimal as the sample size of the climate model simulations and the matrix dimension go to infinity together with a limiting ratio. When EIVs are present, the estimator of the regression coefficients, based on the proposed weight matrix, is asymptotically more efficient than that based on the inverse of the existing linear shrinkage estimator of the error covariance matrix. The performance of the method is confirmed in finite sample simulation studies mimicking realistic situations in terms of the length of the confidence intervals and empirical coverage rates for the regression coefficients. In an application to detection and attribution analyses of the mean temperature at different spatial scales, the method yielded shorter confidence intervals which are important for such analyses in practice.

Data breaches in healthcare have become a substantial concern in recent years and cause millions of dollars in financial losses each year. It is fundamental for government regulators, insurance companies, and stakeholders to understand the breach frequency and the number of affected individuals in each state, as these are directly related to the federal Health Insurance Portability and Accountability Act (HIPAA) and state data breach laws. However, an obstacle to studying data breaches in healthcare is the lack of suitable statistical approaches. We develop a novel multivariate frequency-severity framework to analyze breach frequency and the number of affected individuals at the state level. A mixed effects model is developed to model the square root transformed frequency, and the log-gamma distribution is proposed to capture the skewness and heavy tail exhibited by the distribution of numbers of affected individuals. We further discover a positive nonlinear dependence between the transformed frequency and the log-transformed numbers of affected individuals (i.e., severity). In particular, we propose to use a D-vine copula to capture the multivariate dependence among conditional severities, given frequencies due to its inherent temporal structure and rich bivariate copula families. The rejection sampling technique is developed to simulate the predictive distributions. Both the in-sample and out-of-sample studies show that the proposed multivariate frequency-severity model that accommodates nonlinear dependence has satisfactory fitting and prediction performances.

Discovering association patterns of items from a collection of baskets composed of different items is an important problem in various fields. Assuming that each basket is composed of themes of items randomly sampled from a theme dictionary, the theme dictionary model provides a general framework to achieve efficient association pattern discovery with statistical inference. This paper extends the original theme dictionary model by allowing more than one category of items in a basket and only presence/absence of items is observed for each basket with all quantitative information missing. The extended models can solve a larger range of practical problems that cannot be handled by the original theme dictionary model. Both simulation studies and real data applications confirm the superiority of the proposed methods over the existing ones.

Noncompliance is a common intercurrent event in randomized clinical trials that raises important questions about analytical objectives and approaches. Motivated by the Multiple Risk Factor Intervention Trial (MRFIT), we consider how to estimate the average treatment effect (ATE) in randomized trials with all-or-none compliance. Confounding is a major challenge in estimating the ATE, and conventional methods for confounding adjustment typically require the assumption of no unmeasured confounders which may be difficult to justify. Using randomized treatment assignment as an instrumental variable, the ATE can be identified in the presence of unmeasured confounders under suitable assumptions, including an assumption that limits the effect-modifying activities of unmeasured confounders. We describe and compare several estimation methods based on different modeling assumptions. Some of these methods are able to incorporate information from auxiliary covariates for improved efficiency without introducing bias. The different methods are compared in a simulation study and applied to the MRFIT.

Hypertensive disorders of pregnancy occur in about 10% of pregnant women around the world. Though there is evidence that hypertension impacts maternal cardiac functions, the relation between hypertension and cardiac dysfunctions is only partially understood. The study of this relationship can be framed as a joint inferential problem on multiple populations, each corresponding to a different hypertensive disorder diagnosis, that combines multivariate information provided by a collection of cardiac function indexes. A Bayesian nonparametric approach seems particularly suited for this setup, and we demonstrate it on a dataset consisting of transthoracic echocardiography results of a cohort of Indian pregnant women. We are able to perform model selection, provide density estimates of cardiac function indexes and a latent clustering of patients: these readily interpretable inferential outputs allow to single out modified cardiac functions in hypertensive patients, compared to healthy subjects, and progressively increased alterations with the severity of the disorder. The analysis is based on a Bayesian nonparametric model that relies on a novel hierarchical structure, called symmetric hierarchical Dirichlet process. This is suitably designed so that the mean parameters are identified and used for model selection across populations, a penalization for multiplicity is enforced, and the presence of unobserved relevant factors is investigated through a latent clustering of subjects. Posterior inference relies on a suitable Markov chain Monte Carlo algorithm, and the model behaviour is also showcased on simulated data.

A major issue in the clinical management of epilepsy is the unpredictability of seizures. Yet, traditional approaches to seizure forecasting and risk assessment in epilepsy rely heavily on raw seizure frequencies which are a stochastic measurement of seizure risk. We consider a Bayesian nonhomogeneous hidden Markov model for unsupervised clustering of zero-inflated seizure count data. The proposed model allows for a probabilistic estimate of the sequence of seizure risk states at the individual level. It also offers significant improvement over prior approaches by incorporating a variable selection prior for the identification of clinical covariates that drive seizure risk changes and accommodating highly granular data. For inference, we implement an efficient sampler that employs stochastic search and data augmentation techniques. We evaluate model performance on simulated seizure count data. We then demonstrate the clinical utility of the proposed model by analyzing daily seizure count data from 133 patients with Dravet syndrome collected through the Seizure Tracker${}^{\mathit{TM}}$ system, a patient-reported electronic seizure diary. We report on the dynamics of seizure risk cycling, including validation of several known pharmacologic relationships. We also uncover novel findings characterizing the presence and volatility of risk states in Dravet syndrome which may directly inform counseling to reduce the unpredictability of seizures for patients with this devastating cause of epilepsy.

The ocean is filled with microscopic microalgae, called phytoplankton, which together are responsible for as much photosynthesis as all plants on land combined. Our ability to predict their response to the warming ocean relies on understanding how the dynamics of phytoplankton populations is influenced by changes in environmental conditions. One powerful technique to study the dynamics of phytoplankton is flow cytometry which measures the optical properties of thousands of individual cells per second. Today, oceanographers are able to collect flow cytometry data in real time onboard a moving ship, providing them with fine-scale resolution of the distribution of phytoplankton across thousands of kilometers. One of the current challenges is to understand how these small- and large-scale variations relate to environmental conditions, such as nutrient availability, temperature, light and ocean currents. In this paper we propose a novel sparse mixture of multivariate regressions model to estimate the time-varying phytoplankton subpopulations while simultaneously identifying the specific environmental covariates that are predictive of the observed changes to these subpopulations. We demonstrate the usefulness and interpretability of the approach using both synthetic data and real observations collected on an oceanographic cruise conducted in the northeast Pacific in the spring of 2017.

One formulation of forensic identification of source problems is to determine the source of trace evidence, for instance, glass fragments found on a suspect for a crime. The current state of the science is to compute a Bayes factor comparing the marginal distribution of measurements of trace evidence under two competing propositions for whether or not the unknown source evidence originated from a specific source. The obvious problem with such an approach is the ability to tailor the prior distributions (placed on the features/parameters of the statistical model for the measurements of trace evidence) in favor of the defense or prosecution which is further complicated by the fact that the typical number of measurements of trace evidence is typically sufficiently small that prior choice/specification has a strong influence on the value of the Bayes factor. To remedy this problem of prior specification and choice, we develop an alternative to the Bayes factor, within the framework of generalized fiducial inference, that we term a generalized fiducial factor. Furthermore, we demonstrate empirically, on synthetic and real Netherlands Forensic Institute casework data, deficiencies in Bayes factor and classical/frequentist likelihood ratio approaches.

This paper proposes a novel topological learning framework that integrates networks of different sizes and topology through persistent homology. Such challenging task is made possible through the introduction of a computationally efficient topological loss. The use of the proposed loss bypasses the intrinsic computational bottleneck associated with matching networks. We validate the method in extensive statistical simulations to assess its effectiveness when discriminating networks with different topology. The method is further demonstrated in a twin brain imaging study where we determine if brain networks are genetically heritable. The challenge here is due to the difficulty of overlaying the topologically different functional brain networks obtained from resting-state functional MRI onto the template structural brain network obtained through diffusion MRI.

Preoperative opioid use has been reported to be associated with higher preoperative opioid demand, worse postoperative outcomes, and increased postoperative healthcare utilization and expenditures. Understanding the risk of preoperative opioid use helps establish patient-centered pain management. In the field of machine learning, deep neural network (DNN) has emerged as a powerful means for risk assessment because of its superb prediction power; however, the blackbox algorithms may make the results less interpretable than statistical models. Bridging the gap between the statistical and machine learning fields, we propose a novel interpretable neural network regression (INNER) which combines the strengths of statistical and DNN models. We use the proposed INNER to conduct individualized risk assessment of preoperative opioid use. Intensive simulations and an analysis of 34,186 patients expecting surgery in the Analgesic Outcomes Study (AOS) show that the proposed INNER not only can accurately predict the preoperative opioid use using preoperative characteristics as DNN but also can estimate the patient-specific odds of opioid use without pain and the odds ratio of opioid use for a unit increase in the reported overall body pain, leading to more straightforward interpretations of the tendency to use opioids than DNN. Our results identify the patient characteristics that are strongly associated with opioid use and is largely consistent with the previous findings, providing evidence that INNER is a useful tool for individualized risk assessment of preoperative opioid use.

Subclassification and matching are often used in empirical studies to adjust for observed covariates; however, they are largely restricted to relatively simple study designs with a binary treatment and less developed for designs with a continuous exposure. Matching with exposure doses is particularly useful in instrumental variable designs and in understanding the dose-response relationships. In this article we propose two criteria for optimal subclassification based on subclass homogeneity, in the context of having a continuous exposure dose, and propose an efficient polynomial-time algorithm that is guaranteed to find an optimal subclassification with respect to one criterion and serves as a 2-approximation algorithm for the other criterion. We discuss how to incorporate dose and use appropriate penalties to control the number of subclasses in the design. Via extensive simulations, we systematically compare our proposed design to optimal nonbipartite pair matching and demonstrate that combining our proposed subclassification scheme with regression adjustment helps to reduce model dependence for parametric causal inference with a continuous dose. We apply the new design and associated randomization-based inferential procedure to study the effect of transesophageal echocardiography (TEE) monitoring during coronary artery bypass graft (CABG) surgery on patients’ postsurgery clinical outcomes, using Medicare and Medicaid claims data, and find evidence that TEE monitoring lowers patients’ all-cause 30-day mortality rate.

We consider a large-scale, cross-classified nested (CRON) joint model for modeling customer responses to opening, clicking, and purchasing from promotion emails. Our logistic regression-based joint model contains crossing of promotions and customer effects and allows estimation of the heterogeneous effects of different promotion emails, after adjusting for customer preferences, attributes, and historical behaviors. Using data from an email marketing campaign of an apparel retailer, we exhibit the varying effects of promotions not only based on the contents of the email but also across the three different stages, viz. open, click, and purchase of the conversion funnel. We conduct Bayesian estimation of the parameters in the joint model by using a block Metropolis–Hastings algorithm that not only incorporates nested subsampling to tackle the severe imbalance between conversions and no conversions but also uses additive transformation-based modifications of random walk Metropolis to scale estimation for large numbers of customers. We extend our approach to a segmented cross-classified nested (SCRON) joint model that encompasses the possibility of varying promotion effects across different customer segments. The resultant high-dimensional model is estimated using spike-and-slab priors on the promotion and customer segment interactions. Our nested joint model accounts for the correlations in customer preferences across the conversion funnel. Based on the promotion estimates from the model, we demonstrate how marketers can use different priced, nonpriced, and combination of price and nonprice promotions to increase brand awareness or increase purchases. Comparing estimates from CRON and SCRON models, we display the benefits of targeted marketing by using email promotion lists which are separately optimized for the different customer segments.

Extreme value applications commonly employ regression techniques to capture cross-sectional heterogeneity or time variation in the data. Estimation of the parameters of an extreme value regression model is notoriously challenging due to the small number of observations that are usually available in applications. When repeated extreme measurements are collected on the same individuals, that is, a panel of extremes is available, pooling the observations in groups can improve the statistical inference. We study three data sets related to risk assessment in finance, climate science, and hydrology. In all three cases the problem can be formulated as an extreme value panel regression model with a latent group structure and group-specific parameters. We propose a new algorithm that jointly assigns the individuals to the latent groups and estimates the parameters of the regression model inside each group. Our method efficiently recovers the underlying group structure without prior information, and for the three data sets it provides improved return level estimates and helps answer important domain-specific questions.

The brain is a high-dimensional directed network system, as it consists of many regions as network nodes that exert influence on each other. The directed influence exerted by one region on another is referred to as directed connectivity. We aim to reveal whole-brain directed networks based on resting-state functional magnetic resonance imaging (fMRI) data of many subjects. However, it is both statistically and computationally challenging to produce scientifically meaningful estimates of whole-brain directed networks. To address the statistical modeling challenge, we assume modular brain networks which reflect functional specialization and functional integration of the brain. We address the computational challenge by developing a variational Bayesian method to estimate the new model. We apply our method to resting-state fMRI data of many subjects and identify modules and directed connections in whole-brain directed networks. The identified modules are accordant with functional brain systems specialized for different functions. We also detect directed connections between functionally specialized modules, which is not attainable by existing network methods, based on functional connectivity. In summary, this paper presents a new computationally efficient and flexible method for directed network studies of the brain as well as new scientific findings regarding the functional organization of the human brain.

Many environments within the human body host a collection of micro-organisms called microbiota. Recent findings have linked the composition of the microbiota to the development of different human diseases, including cancer. Motivated by a recent colorectal cancer (crc) study, we investigate the effect of clinical factors and diet-related covariates on the microbiota compositions; for the patients enrolled in this study, microbiota abundance counts are collected from three different districts, namely, tumor, fecal and salivary samples. Building upon the Dirichlet-multinomial regression framework, we develop a high-dimensional Bayesian hierarchical model that exploits subject-specific regression coefficients to simultaneously borrow strength across districts and include complex interactions between diet and clinical factors if supported by the data. The proposed method identifies relevant associations through model selection priors and thresholding mechanisms. Posterior inference is performed via a Markov chain Monte Carlo algorithm. We use simulation studies to assess the performance of our method, and found our approach to outperform competing methods that do not account for complex interactions. Finally, a thorough analysis of the crc data illustrates the benefits of the proposed approach.

Accurate spatiotemporal modeling of conditions leading to moderate and large wildfires provides better understanding of mechanisms driving fire-prone ecosystems and improves risk management. Here, we develop a joint model for the occurrence intensity and the wildfire size distribution, by combining extreme-value theory and point processes within a novel Bayesian hierarchical model, and use it to study daily summer wildfire data for the French Mediterranean basin during 1995–2018. The occurrence component models wildfire ignitions as a spatiotemporal log-Gaussian Cox process. Burnt areas are numerical marks attached to points and are considered as extreme if they exceed a high threshold. The size component is a two-component mixture varying in space and time that jointly models moderate and extreme fires. We capture nonlinear influence of covariates (Fire Weather Index, forest cover) through component-specific smooth functions which may vary with season. We propose estimating shared random effects between model components to reveal and interpret common drivers of different aspects of wildfire activity. This increases parsimony and reduces estimation uncertainty, giving better predictions. Specific stratified subsampling of zero counts is implemented to cope with large observation vectors. We compare and validate models through predictive scores and visual diagnostics. Our methodology provides a holistic approach to explaining and predicting the drivers of wildfire activity and associated uncertainties.

The coronavirus (COVID-19) global pandemic has made a significant impact on people’s social activities. Cell phone mobility data provide unique and rich information on studying this impact. The motivating dataset of this study is the daily leaving-home index data at Harris County in Texas provided by SafeGraph. To study changes in daily leaving-home index and how they relate to public policy and sociodemographic variables, we propose a new Bayesian wavelet model for modeling and clustering spatial functional data, where domain partitioning is achieved by operating on the spanning trees. The resulting clusters can have arbitrary shapes and are spatially contiguous in the input domain. An efficient tailored reversible jump Markov chain Monte Carlo algorithm is proposed to implement the model. The method is applied to the spatial functional data of the daily percentages of people who left home. We focus on the time period covering both lockdown and phased reopening in Texas during the COVID-19 pandemic and study the changing behaviors of those functional curves. By linking the clustering results with the sociodemographic information, we identify several covariates of census blocks that have a noticeable impact on the clustering patterns of people’s mobility behaviors.

Jam formation is a problem that may occur when granular material is discharged by gravity from a silo. The estimation of the minimum outlet size, which guarantees that the time to the next jamming event is long enough, can be crucial in the industry. The time is modeled by an exponential distribution with two unknown parameters, and this goal translates to precise estimation of a nonlinear transformation of the parameters. We obtain c-optimum experimental designs with that purpose, applying the graphic Elfving method. Because the optimal experimental designs depend on the nominal values of the parameters, we conduct a sensitivity analysis on our dataset. Finally, a simulation study checks the performance of the approximations, first with the Fisher Information matrix, then with the linearization of the function to be estimated. The results are useful for experimenting in a laboratory and then translating the results to a real scenario. From the application we develop a general methodology for estimating a one-dimensional transformation of the parameters of a nonlinear model.

In immunology studies, flow cytometry is a commonly used multivariate single-cell assay. One key goal in flow cytometry analysis is to detect the immune cells responsive to certain stimuli. Statistically, this problem can be translated into comparing two protein expression probability density functions (pdfs) before and after the stimulus; the goal is to pinpoint the regions where these two pdfs differ. Further screening of these differential regions can be performed to identify enriched sets of responsive cells. In this paper we model identifying differential density regions as a multiple testing problem. First, we partition the sample space into small bins. In each bin we form a hypothesis to test the existence of differential pdfs. Second, we develop a novel multiple testing method, called TEAM (testing on the aggregation tree method), to identify those bins that harbor differential pdfs while controlling the false discovery rate (FDR) under the desired level. TEAM embeds the testing procedure into an aggregation tree to test from fine- to coarse-resolution. The procedure achieves the statistical goal of pinpointing density differences to the smallest possible regions. TEAM is computationally efficient, capable of analyzing large flow cytometry data sets in much shorter time compared with competing methods. We applied TEAM and competing methods on a flow cytometry data set to identify T cells responsive to the cytomegalovirus (CMV)-pp65 antigen stimulation. With additional downstream screening, TEAM successfully identified enriched sets containing monofunctional, bifunctional, and polyfunctional T cells. Competing methods either did not finish in a reasonable time frame or provided less interpretable results. Numerical simulations and theoretical justifications demonstrate that TEAM has asymptotically valid, powerful, and robust performance. Overall, TEAM is a computationally efficient and statistically powerful algorithm that can yield meaningful biological insights in flow cytometry studies.

A risk-limiting election audit (RLA) offers a statistical guarantee: if the reported electoral outcome is incorrect, the audit has a known maximum chance (the risk limit) of not correcting it before it becomes final. BRAVO (Lindeman, Stark and Yates (In Proceedings of the 2011 Electronic Voting Technology Workshop/Workshop on Trustworthy Elections (EVT/WOTE’11) (2012) USENIX)), based on Wald’s sequential probability ratio test for the Bernoulli parameter, is the simplest and most widely tried method for RLAs, but it has limitations. It cannot accommodate sampling without replacement or stratified sampling which can improve efficiency and are sometimes required by law. It applies only to ballot-polling audits which are less efficient than comparison audits. It applies to plurality, majority, supermajority, proportional representation, and instant-runoff voting (IRV, using RAIRE (Blom, Stuckey and Teague (In Electronic Voting (2018) 17–34 Springer))) but not to other social choice functions for which there are RLA methods. And while BRAVO has the smallest expected sample size among sequentially valid ballot-polling-with-replacement methods when the reported vote shares are exactly correct, it can require arbitrarily large samples when the reported reported winner(s) really won but the reported vote shares are incorrect. ALPHA is a simple generalization of BRAVO that: (i) works for sampling with and without replacement, with and without weights, with and without stratification, and for Bernoulli sampling; (ii) works not only for ballot polling but also for ballot-level comparison, batch polling, and batch-level comparison audits; (iii) works for all social choice functions covered by SHANGRLA (Stark (In Financial Cryptography and Data Security (2020) Springer)), including approval voting, STAR-Voting, proportional representation schemes, such as D’Hondt and Hamilton, IRV, Borda count, and all scoring rules, and (iv) in situations where both ALPHA and BRAVO apply, requires smaller samples than BRAVO when the reported vote shares are wrong but the outcome is correct—five orders of magnitude in some examples. ALPHA includes the family of betting martingale tests in RiLACS (Waudby-Smith, Stark and Ramdas (In Electronic Voting. E-Vote-ID 2021 (2021) Springer)) with a different betting strategy parametrized as an estimator of the population mean and explicit flexibility to accommodate sampling weights and population bounds that change with each draw. A Python implementation is provided.

Statistical model development is a central feature of many scientific investigations with a vast methodological landscape. However, uncertainty in the model development process has received less attention and is frequently resolved nonrigorously through beliefs about generalizability, practical usefulness, and computational ease. This is particularly problematic in settings of abundant data, such as clickstream data, as model selection routinely admits multiple models and imposes a source of uncertainty, unacknowledged and unknown by many, on all postselection conclusions. Regression models, based on the beta distribution, are a class of nonlinear models, attractive because of their great flexibility and potential explanatory power, but have not been investigated from the standpoint of multimodel uncertainty and model averaging. For this reason a formalized tool that can combine model selection uncertainty and beta regression modeling is presented in this work. The tool combines bootstrap model averaging, model selection, and asymptotic theory to yield a procedure that can perform joint modeling of the mean and precision parameters, capture sources of variability in the data, and achieve more accurate claims of estimate precision, variable importance, generalization performance, and model stability. Practical utility of the tool is demonstrated through a study of model selection consistency and variable importance in average exit and bounce rate statistical models. This work emphasizes the necessity of a departure from the all-too-common practice of ignoring model selection uncertainty and introduces an accessible technique to handle frequently neglected aspects of the modeling pipeline.

It is widely speculated that auditors’ public forecasts of bankruptcy are, at least in part, self-fulfilling prophecies in the sense that they actually cause bankruptcies that would not have otherwise occurred. This conjecture is hard to prove, however, because the strong association between bankruptcies and bankruptcy forecasts could simply indicate that auditors are skillful forecasters with unique access to highly predictive covariates. In this paper we investigate the causal effect of bankruptcy forecasts on bankruptcy using nonparametric sensitivity analysis. We contrast our analysis with two alternative approaches: a linear bivariate probit model with an endogenous regressor and a recently developed bound on risk ratios called E-values. Additionally, our machine learning approach incorporates a monotonicity constraint corresponding to the assumption that bankruptcy forecasts do not make bankruptcies less likely. Finally, a tree-based posterior summary of the treatment effect estimates allows us to explore which observable firm characteristics moderate the inducement effect.

Continuous time discrete state models are a valuable tool for explaining animal movement. However, data collection to fit such models over a specified window of time can be misaligned with the actual realization of the movement process. This necessitates approximate model fitting, at present, through approximate imputation distributions (AIDs). Here, we propose a direct time-discretization approximation to the likelihood. The approach employs familiar ideas from hidden Markov modeling. Computation is implemented through the induced infinitesimal generator matrix. Linearization of this matrix expedites computation time. Through simulation and a real data application involving whale movement, we demonstrate that this model fitting strategy can outperform AID approaches.

Risk assessments to help inform criminal justice decisions have been used in the United States since the 1920s. Over the past several years, statistical learning risk algorithms have been introduced amid much controversy about fairness, transparency and accuracy. In this paper we focus on accuracy for a large department of probation and parole that is considering a major revision of its current, statistical learning risk methods. Because the content of each offender’s supervision is substantially shaped by a forecast of subsequent conduct, forecasts have real consequences. Here, we consider the probability that risk forecasts are correct. We augment standard statistical learning estimates of forecasting uncertainty (i.e., confusion tables) with uncertainty estimates from nested conformal prediction sets. In a demonstration of concept using data from the department of probation and parole, we show that the standard uncertainty measures and uncertainty measures from nested conformal prediction sets can differ dramatically in formulation and output. We also provide a modification of nested conformal, called the localized conformal method, to match confusion tables more closely. A strong case can be made favoring the nested conformal approach. As best, we can tell our formulation of such comparisons and consequent recommendations is novel.

A common difficulty in data analysis is how to handle categorical predictors with a large number of levels or categories. Few proposals have been developed to tackle this important and frequent problem. We introduce a generative model that simultaneously carries out the model fitting and the aggregation of the categorical levels into larger groups. We represent the categorical predictor by a graph where the nodes are the categories and establish a probability distribution over meaningful partitions of this graph. Conditionally on the observed data, we obtain a posterior distribution for the levels aggregation, allowing the inference about the most probable clustering for the categories. Simultaneously, we extract inference about all the other regression model parameters. We compare our and state-of-art methods showing that it has equally good predictive performance and more interpretable results. Our approach balances out accuracy vs. interpretability, a current important concern in statistics and machine learning.

The growing number of available treatment options has led to urgent needs for reliable answers when choosing the best course of treatment for a patient. As it is often infeasible to compare a large number of treatments in a single randomized controlled trial, multivariate network meta-analyses (NMAs) are used to synthesize evidence from trials of a subset of the treatments, where both efficacy and safety related outcomes are considered simultaneously. However, these large-scale multiple-outcome NMAs have created challenges to existing methods due to the increasing complexity of the unknown correlations between outcomes and treatment comparisons. In this paper, we proposed a new framework for PAtient-centered treatment ranking via Large-scale Multivariate network meta-analysis, termed as PALM, which includes a parsimonious modeling approach, a fast algorithm for parameter estimation and inference, a novel visualization tool for presenting multivariate outcomes, termed as the origami plot, as well as personalized treatment ranking procedures taking into account the individual’s considerations on multiple outcomes. In application to an NMA that compares 14 treatment options for labor induction, we provided a comprehensive illustration of the proposed framework and demonstrated its computational efficiency and practicality, and we obtained new insights and evidence to support patient-centered clinical decision making.

In social and economic surveys, it can be difficult to directly reach units of the target population, and indirect sampling is often advocated to solve this issue. In indirect sampling, the sample is drawn from a frame population that is linked to the target population, and estimation of target population parameters is typically achieved through the generalized weight share method (GWSM). This method provides a weight, for every unit of the target population, that depends on the one hand, on the sampling weights in the frame population and, on the other hand, on the link weights between the frame population and the target population. In the present study, we focus on the situation in which the units from the frame population are linked to one and only one unit from the target population (Many-to-One case). This situation is encountered at the French postal service where addresses are sampled instead of postman rounds. We aim at understanding of the impact of the link weights on the efficiency of the GWSM estimators. We derive variance expressions and optimality results for a large class of sampling designs. Moreover, we note that the Many-to-One case can lead to too many links to observe. We alleviate the problem by introducing an intermediate population and double indirect sampling. The question of the loss of precision in this situation is discussed in detail through theoretical results and simulations. These findings help to explain the loss of precision of double GWSM estimators observed recently at the French postal service.

Measurements of many biological processes are characterized by an initial trend period followed by an equilibrium period. Scientists may wish to quantify features of the two periods as well as the timing of the change point. Specifically, we are motivated by problems in the study of electrical cell-substrate impedance sensing (ECIS) data. ECIS is a popular new technology which measures cell behavior noninvasively. Previous studies using ECIS data have found that different cell types can be classified by their equilibrium behavior. However, it can be challenging to identify when equilibrium has been reached and to quantify the relevant features of cells’ equilibrium behavior. In this paper we assume that measurements during the trend period are independent deviations from a smooth nonlinear function of time, and that measurements during the equilibrium period are characterized by a simple long memory model. We propose a method to simultaneously estimate the parameters of the trend and equilibrium processes and locate the change point between the two. We find that this method performs well in simulations and in practice. When applied to ECIS data, it produces estimates of change points and measures of cell equilibrium behavior which offer improved classification of infected and uninfected cells.

Insurance claims’ data is an increasingly important health policy research resource, given its longitudinal assessment of cancer care clinical outcomes. Population-level information on medical cost trajectory from disease diagnosis to terminal events, such as death, specifically interests policy makers. Estimating the mean cost trajectory has statistical challenges. The shape of the trajectory is usually highly nonlinear with varying durations, depending on the diagnosis-to-death population time distribution. The terminal event may be right censored, resulting in missing subsequent costs. Medical costs often have skewed distributions with zero inflation and heteroscedasticity which may not fit well with the commonly used parametric family of distributions. In this paper we propose a flexible semiparametric model to address challenges without imposing a cost data distributional assumption. The estimation procedure is based on generalized estimating equations with censored covariates. The proposed model adopts a bivariate surface that quantifies the interrelationship between longitudinal medical costs and survival, and results in the nonlinear population mean cost trajectories given survival time. We develop a novel generalized estimating equations algorithm to accommodate covariates subject to right censoring without fully specifying the joint distribution of the cost and survival data. We provide theoretical and simulation-based justification for the proposed approach and apply the methods to estimate prostate cancer patient cost trajectories from the Surveillance, Epidemiology, and End Results (SEER)-Medicare linked database.