Search

Comparative biologists are often interested in inferring covariation between multiple biological traits sampled across numerous related taxa. To properly study these relationships, we must control for the shared evolutionary history of the taxa to avoid spurious inference. An additional challenge arises as obtaining a full suite of measurements becomes increasingly difficult with increasing taxa. This generally necessitates data imputation or integration, and existing control techniques typically scale poorly as the number of taxa increases. We propose an inference technique that integrates out missing measurements analytically and scales linearly with the number of taxa by using a post-order traversal algorithm under a multivariate Brownian diffusion (MBD) model to characterize trait evolution. We further exploit this technique to extend the MBD model to account for sampling error or non-heritable residual variance. We test these methods to examine mammalian life history traits, prokaryotic genomic and phenotypic traits, and HIV infection traits. We find computational efficiency increases that top two orders-of-magnitude over current best practices. While we focus on the utility of this algorithm in phylogenetic comparative methods, our approach generalizes to solve long-standing challenges in computing the likelihood for matrix-normal and multivariate normal distributions with missing data at scale.

The exponential-family random graph models (ERGMs) have emerged as animportant framework for modeling social networks for a wide variety ofrelational types. ERGMs for valued networks are less well-developed than theirunvalued counterparts, and pose particular computational challenges. Networkdata with edge values on the non-negative integers (count-valued networks) isan important such case, with examples ranging from the magnitude of migrationand trade flows between places to the frequency of interactions and encountersbetween individuals. Here, we propose an efficient parallelable subsampledmaximum pseudo-likelihood estimation (MPLE) scheme for count-valued ERGMs, andcompare its performance with existing Contrastive Divergence (CD) and MonteCarlo Maximum Likelihood Estimation (MCMLE) approaches via a simulation studybased on migration flow networks in two U.S. states. Our results suggest thatedge value variance is a key factor in method performance, while network sizemainly influences their relative merits in computational time. Forsmall-variance networks, all methods perform well in point estimations while CDgreatly overestimates uncertainties, and MPLE underestimates them fordependence terms; all methods have fast estimation for small networks, but CDand subsampled multi-core MPLE provides speed advantages as network sizeincreases. For large-variance networks, both MPLE and MCMLE offer high-qualityestimates of coefficients and their uncertainty, but MPLE is significantlyfaster than MCMLE; MPLE is also a better seeding method for MCMLE than CD, asthe latter makes MCMLE more prone to convergence failure.

Across research disciplines, cluster randomized trials (CRTs) are commonlyimplemented to evaluate interventions delivered to groups of participants, suchas communities and clinics. Despite advances in the design and analysis ofCRTs, several challenges remain. First, there are many possible ways to specifythe causal effect of interest (e.g., at the individual-level or at thecluster-level). Second, the theoretical and practical performance of commonmethods for CRT analysis remain poorly understood. Here, we present a generalframework to formally define an array of causal effects in terms of summarymeasures of counterfactual outcomes. Next, we provide a comprehensive overviewof CRT estimators, including the t-test, generalized estimating equations(GEE), augmented-GEE, and targeted maximum likelihood estimation (TMLE). Usingfinite sample simulations, we illustrate the practical performance of theseestimators for different causal effects and when, as commonly occurs, there arelimited numbers of clusters of different sizes. Finally, our application todata from the Preterm Birth Initiative (PTBi) study demonstrates the real-worldimpact of varying cluster sizes and targeting effects at the cluster-level orat the individual-level. Specifically, the relative effect of the PTBIintervention was 0.81 at the cluster-level, corresponding to a 19% reduction inoutcome incidence, and was 0.66 at the individual-level, corresponding to a 34%reduction in outcome risk. Given its flexibility to estimate a variety ofuser-specified effects and ability to adaptively adjust for covariates forprecision gains while maintaining Type-I error control, we conclude TMLE is apromising tool for CRT analysis.

This paper develops two orthogonal contributions to scalable sparse regression for competing risks time-to-event data. First, we study and accelerate the broken adaptive ridge method (BAR), a surrogate ℓ 0-based iteratively reweighted ℓ 2-penalization algorithm that achieves sparsity in its limit, in the context of the Fine-Gray (1999) proportional subdistributional hazards (PSH) model. In particular, we derive a new algorithm for BAR regression, named cycBAR, that performs cyclic update of each coordinate using an explicit thresholding formula. The new cycBAR algorithm effectively avoids fitting multiple reweighted ℓ 2-penalizations and thus yields impressive speedups over the original BAR algorithm. Second, we address a pivotal computational issue related to fitting the PSH model. Specifically, the computation costs of the log-pseudo likelihood and its derivatives for PSH model grow at the rate of O(n 2) with the sample size n in current implementations. We propose a novel forward-backward scan algorithm that reduces the computation costs to O(n). The proposed method applies to both unpenalized and penalized estimation for the PSH model and has exhibited drastic speedups over current implementations. Finally, combining the two algorithms can yields > 1, 000 fold speedups over the original BAR algorithm. Illustrations of the impressive scalability of our proposed algorithm for large competing risks data are given using both simulations and a United States Renal Data System data. Supplementary materials for this article are available online.

In this paper, we focus on identifying differentially activated brain regionsusing a light sheet fluorescence microscopy - a recently developed techniquefor whole-brain imaging. Most existing statistical methods solve this problemby partitioning the brain regions into two classes: significantly andnon-significantly activated. However, for the brain imaging problem at thecenter of our study, such binary grouping may provide overly simplisticdiscoveries by filtering out weak but important signals, that are typicallyadulterated by the noise present in the data. To overcome this limitation, weintroduce a new Bayesian approach that allows classifying the brain regionsinto several tiers with varying degrees of relevance. Our approach is based ona combination of shrinkage priors - widely used in regression and multiplehypothesis testing problems - and mixture models - commonly used in model-basedclustering. In contrast to the existing regularizing prior distributions, whichuse either the spike-and-slab prior or continuous scale mixtures, our class ofpriors is based on a discrete mixture of continuous scale mixtures and devisesa cluster-shrinkage version of the Horseshoe prior. As a result, our approachprovides a more general setting for Bayesian sparse estimation, drasticallyreduces the number of shrinkage parameters needed, and creates a framework forsharing information across units of interest. We show that this approach leadsto more biologically meaningful and interpretable results in our brain imagingproblem, since it allows the discrimination between active and inactiveregions, while at the same time ranking the discoveries into clustersrepresenting tiers of similar importance.

Optical imaging of genetically encoded calcium indicators is a powerful toolto record the activity of a large number of neurons simultaneously over a longperiod of time from freely behaving animals. However, determining the exacttime at which a neuron spikes and estimating the underlying firing rate fromcalcium fluorescence data remains challenging, especially for calcium imagingdata obtained from a longitudinal study. We propose a multi-trial time-varying$\ell_0$ penalized method to jointly detect spikes and estimate firing rates byrobustly integrating evolving neural dynamics across trials. Our simulationstudy shows that the proposed method performs well in both spike detection andfiring rate estimation. We demonstrate the usefulness of our method on calciumfluorescence trace data from two studies, with the first study showingdifferential firing rate functions between two behaviors and the second studyshowing evolving firing rate function across trials due to learning.

Relaxed random walk (RRW) models of trait evolution introduce branch-specific rate multipliers to modulate the variance of a standard Brownian diffusion process along a phylogeny and more accurately model overdispersed biological data. Increased taxonomic sampling challenges inference under RRWs as the number of unknown parameters grows with the number of taxa. To solve this problem, we present a scalable method to efficiently fit RRWs and infer this branch-specific variation in a Bayesian framework. We develop a Hamiltonian Monte Carlo (HMC) sampler to approximate the high-dimensional, correlated posterior that exploits a closed-form evaluation of the gradient of the trait data log-likelihood with respect to all branch-rate multipliers simultaneously. Our gradient calculation achieves computational complexity that scales only linearly with the number of taxa under study. We compare the efficiency of our HMC sampler to the previously standard univariable Metropolis-Hastings approach while studying the spatial emergence of the West Nile virus in North America in the early 2000s. Our method achieves at least a 6-fold speed increase over the univariable approach. Additionally, we demonstrate the scalability of our method by applying the RRW to study the correlation between five mammalian life history traits in a phylogenetic tree with $3650$ tips.[Bayesian inference; BEAST; Hamiltonian Monte Carlo; life history; phylodynamics, relaxed random walk.].

The haematopoietic system has a highly regulated and complex structure in which cells are organized to successfully create and maintain new blood cells. It is known that feedback regulation is crucial to tightly control this system, but the specific mechanisms by which control is exerted are not completely understood. In this work, we aim to uncover the underlying mechanisms in haematopoiesis by conducting perturbation experiments, where animal subjects are exposed to an external agent in order to observe the system response and evolution. We have developed a novel Bayesian hierarchical framework for optimal design of perturbation experiments and proper analysis of the data collected. We use a deterministic model that accounts for feedback and feedforward regulation on cell division rates and self-renewal probabilities. A significant obstacle is that the experimental data are not longitudinal, rather each data point corresponds to a different animal. We overcome this difficulty by modelling the unobserved cellular levels as latent variables. We then use principles of Bayesian experimental design to optimally distribute time points at which the haematopoietic cells are quantified. We evaluate our approach using synthetic and real experimental data and show that an optimal design can lead to better estimates of model parameters.

The purpose of this paper is to construct confidence intervals for theregression coefficients in high-dimensional Cox proportional hazards regressionmodels where the number of covariates may be larger than the sample size. Ourdebiased estimator construction is similar to those in Zhang and Zhang (2014)and van de Geer et al. (2014), but the time-dependent covariates and censoredrisk sets introduce considerable additional challenges. Our theoreticalresults, which provide conditions under which our confidence intervals areasymptotically valid, are supported by extensive numerical experiments.

Staggered adoption of policies by different units at different times createspromising opportunities for observational causal inference. Estimation remainschallenging, however, and common regression methods can give misleadingresults. A promising alternative is the synthetic control method (SCM), whichfinds a weighted average of control units that closely balances the treatedunit's pre-treatment outcomes. In this paper, we generalize SCM, originallydesigned to study a single treated unit, to the staggered adoption setting. Wefirst bound the error for the average effect and show that it depends on boththe imbalance for each treated unit separately and the imbalance for theaverage of the treated units. We then propose "partially pooled" SCM weights tominimize a weighted combination of these measures; approaches that focus onlyon balancing one of the two components can lead to bias. We extend thisapproach to incorporate unit-level intercept shifts and auxiliary covariates.We assess the performance of the proposed method via extensive simulations andapply our results to the question of whether teacher collective bargainingleads to higher school spending, finding minimal impacts. We implement theproposed method in the augsynth R package.

The synthetic control method (SCM) is a popular approach for estimating theimpact of a treatment on a single unit in panel data settings. The "syntheticcontrol" is a weighted average of control units that balances the treatedunit's pre-treatment outcomes as closely as possible. A critical feature of theoriginal proposal is to use SCM only when the fit on pre-treatment outcomes isexcellent. We propose Augmented SCM as an extension of SCM to settings wheresuch pre-treatment fit is infeasible. Analogous to bias correction for inexactmatching, Augmented SCM uses an outcome model to estimate the bias due toimperfect pre-treatment fit and then de-biases the original SCM estimate. Ourmain proposal, which uses ridge regression as the outcome model, directlycontrols pre-treatment fit while minimizing extrapolation from the convex hull.This estimator can also be expressed as a solution to a modified syntheticcontrols problem that allows negative weights on some donor units. We bound theestimation error of this approach under different data generating processes,including a linear factor model, and show how regularization helps to avoidover-fitting to noise. We demonstrate gains from Augmented SCM with extensivesimulation studies and apply this framework to estimate the impact of the 2012Kansas tax cuts on economic growth. We implement the proposed method in the newaugsynth R package.

Markov models of character substitution on phylogenies form the foundation of phylogenetic inference frameworks. Early models made the simplifying assumption that the substitution process is homogeneous over time and across sites in the molecular sequence alignment. While standard practice adopts extensions that accommodate heterogeneity of substitution rates across sites, heterogeneity in the process over time in a site-specific manner remains frequently overlooked. This is problematic, as evolutionary processes that act at the molecular level are highly variable, subjecting different sites to different selective constraints over time, impacting their substitution behavior. We propose incorporating time variability through Markov-modulated models (MMMs), which extend covarion-like models and allow the substitution process (including relative character exchange rates as well as the overall substitution rate) at individual sites to vary across lineages. We implement a general MMM framework in BEAST, a popular Bayesian phylogenetic inference software package, allowing researchers to compose a wide range of MMMs through flexible XML specification. Using examples from bacterial, viral, and plastid genome evolution, we show that MMMs impact phylogenetic tree estimation and can substantially improve model fit compared to standard substitution models. Through simulations, we show that marginal likelihood estimation accurately identifies the generative model and does not systematically prefer the more parameter-rich MMMs. To mitigate the increased computational demands associated with MMMs, our implementation exploits recent developments in BEAGLE, a high-performance computational library for phylogenetic inference. [Bayesian inference; BEAGLE; BEAST; covarion, heterotachy; Markov-modulated models; phylogenetics.].

Many scientific and engineering challenges -- ranging from pharmacokineticdrug dosage allocation and personalized medicine to marketing mix (4Ps)recommendations -- require an understanding of the unobserved heterogeneity inorder to develop the best decision making-processes. In this paper, we developa hypothesis test and the corresponding p-value for testing for thesignificance of the homogeneous structure in linear mixed models. A robustmatching moment construction is used for creating a test that adapts to thesize of the model sparsity. When unobserved heterogeneity at a cluster level isconstant, we show that our test is both consistent and unbiased even when thedimension of the model is extremely high. Our theoretical results rely on a newfamily of adaptive sparse estimators of the fixed effects that do not requireconsistent estimation of the random effects. Moreover, our inference results donot require consistent model selection. We showcase that moment matching can beextended to nonlinear mixed effects models and to generalized linear mixedeffects models. In numerical and real data experiments, we find that thedeveloped method is extremely accurate, that it adapts to the size of theunderlying model and is decidedly powerful in the presence of irrelevantcovariates.

Non-ignorable dropout is common in studies with long follow-up time, and it can bias study results unless handled carefully in the study design and the statistical analysis. A double-sampling design allocates additional resources to pursue a subsample of the dropouts and find out their outcomes, which can address potential biases due to non-ignorable dropout. It is desirable to construct semiparametric estimators for the double-sampling design because of their robustness properties. However, obtaining such semiparametric estimators remains a challenge due to the requirement of the analytic form of the efficient influence function (EIF), the derivation of which can be ad hoc and difficult for the double-sampling design. Recent work has shown how the derivation of EIF can be made deductive and computerizable using the functional derivative representation of the EIF in nonparametric models. This approach, however, requires deriving the mixture of a continuous distribution and a point mass, which can itself be challenging for complicated problems such as the double-sampling design. We propose semiparametric estimators for the survival probability in double-sampling designs by generalizing the deductive and computerizable estimation approach. In particular, we propose to build the semiparametric estimators based on a discretized support structure, which approximates the possibly continuous observed data distribution and circumvents the derivation of the mixture distribution. Our approach is deductive in the sense that it is expected to produce semiparametric locally efficient estimators within finite steps without knowledge of the EIF. We apply the proposed estimators to estimating the mortality rate in a double-sampling design component of the President’s Emergency Plan for AIDS Relief (PEPFAR) program. We evaluate the impact of double-sampling selection criteria on the mortality rate estimates. Simulation studies are conducted to evaluate the robustness of the proposed estimators.

Modeling correlation (and covariance) matrices can be challenging due to the positive-definiteness constraint and potential high-dimensionality. Our approach is to decompose the covariance matrix into the correlation and variance matrices and propose a novel Bayesian framework based on modeling the correlations as products of unit vectors. By specifying a wide range of distributions on a sphere (e.g. the squared-Dirichlet distribution), the proposed approach induces flexible prior distributions for covariance matrices (that go beyond the commonly used inverse-Wishart prior). For modeling real-life spatio-temporal processes with complex dependence structures, we extend our method to dynamic cases and introduce unit-vector Gaussian process priors in order to capture the evolution of correlation among components of a multivariate time series. To handle the intractability of the resulting posterior, we introduce the adaptive Δ-Spherical Hamiltonian Monte Carlo. We demonstrate the validity and flexibility of our proposed framework in a simulation study of periodic processes and an analysis of rat's local field potential activity in a complex sequence memory task.

Coalescent theory combined with statistical modeling allows us to estimate effective population size fluctuations from molecular sequences of individuals sampled from a population of interest. When sequences are sampled serially through time and the distribution of the sampling times depends on the effective population size, explicit statistical modeling of sampling times improves population size estimation. Previous work assumed that the genealogy relating sampled sequences is known and modeled sampling times as an inhomogeneous Poisson process with log-intensity equal to a linear function of the log-transformed effective population size. We improve this approach in two ways. First, we extend the method to allow for joint Bayesian estimation of the genealogy, effective population size trajectory, and other model parameters. Next, we improve the sampling time model by incorporating additional sources of information in the form of time-varying covariates. We validate our new modeling framework using a simulation study and apply our new methodology to analyses of population dynamics of seasonal influenza and to the recent Ebola virus outbreak in West Africa.

Phylodynamics is an area of population genetics that uses genetic sequence data to estimate past population dynamics. Modern state-of-the-art Bayesian nonparametric methods for recovering population size trajectories of unknown form use either change-point models or Gaussian process priors. Change-point models suffer from computational issues when the number of change-points is unknown and needs to be estimated. Gaussian process-based methods lack local adaptivity and cannot accurately recover trajectories that exhibit features such as abrupt changes in trend or varying levels of smoothness. We propose a novel, locally adaptive approach to Bayesian nonparametric phylodynamic inference that has the flexibility to accommodate a large class of functional behaviors. Local adaptivity results from modeling the log-transformed effective population size a priori as a horseshoe Markov random field, a recently proposed statistical model that blends together the best properties of the change-point and Gaussian process modeling paradigms. We use simulated data to assess model performance, and find that our proposed method results in reduced bias and increased precision when compared to contemporary methods. We also use our models to reconstruct past changes in genetic diversity of human hepatitis C virus in Egypt and to estimate population size changes of ancient and modern steppe bison. These analyses show that our new method captures features of the population size trajectories that were missed by the state-of-the-art methods.

The human body generates a diverse set of high affinity antibodies, the soluble form of B cell receptors (BCRs), that bind to and neutralize invading pathogens. The natural development of BCRs must be understood in order to design vaccines for highly mutable pathogens such as influenza and HIV. BCR diversity is induced by naturally occurring combinatorial "V(D)J" rearrangement, mutation, and selection processes. Most current methods for BCR sequence analysis focus on separately modeling the above processes. Statistical phylogenetic methods are often used to model the mutational dynamics of BCR sequence data, but these techniques do not consider all the complexities associated with B cell diversification such as the V(D)J rearrangement process. In particular, standard phylogenetic approaches assume the DNA bases of the progenitor (or "naive") sequence arise independently and according to the same distribution, ignoring the complexities of V(D)J rearrangement. In this paper, we introduce a novel approach to Bayesian phylogenetic inference for BCR sequences that is based on a phylogenetic hidden Markov model (phylo-HMM). This technique not only integrates a naive rearrangement model with a phylogenetic model for BCR sequence evolution but also naturally accounts for uncertainty in all unobserved variables, including the phylogenetic tree, via posterior distribution sampling.

The use of high-dimensional data for targeted therapeutic interventionsrequires new ways to characterize the heterogeneity observed across subgroupsof a specific population. In particular, models for partially exchangeable dataare needed for inference on nested datasets, where the observations are assumedto be organized in different units and some sharing of information is requiredto learn distinctive features of the units. In this manuscript, we propose anested Common Atoms Model (CAM) that is particularly suited for the analysis ofnested datasets where the distributions of the units are expected to differonly over a small fraction of the observations sampled from each unit. Theproposed CAM allows a two-layered clustering at the distributional andobservational level and is amenable to scalable posterior inference through theuse of a computationally efficient nested slice-sampler algorithm. We furtherdiscuss how to extend the proposed modeling framework to handle discretemeasurements, and we conduct posterior inference on a real microbiome datasetfrom a diet swap study to investigate how the alterations in intestinalmicrobiota composition are associated with different eating habits. We furtherinvestigate the performance of our model in capturing true distributionalstructures in the population by means of a simulation study.

We propose an evolutionary state space model (E-SSM) for analyzing high dimensional brain signals whose statistical properties evolve over the course of a non-spatial memory experiment. Under E-SSM, brain signals are modeled as mixtures of components (e.g., AR(2) process) with oscillatory activity at pre-defined frequency bands. To account for the potential non-stationarity of these components (since the brain responses could vary throughout the entire experiment), the parameters are allowed to vary over epochs. Compared with classical approaches such as independent component analysis and filtering, the proposed method accounts for the entire temporal correlation of the components and accommodates non-stationarity. For inference purpose, we propose a novel computational algorithm based upon using Kalman smoother, maximum likelihood and blocked resampling. The E-SSM model is applied to simulation studies and an application to a multi-epoch local field potentials (LFP) signal data collected from a non-spatial (olfactory) sequence memory task study. The results confirm that our method captures the evolution of the power for different components across different phases in the experiment and identifies clusters of electrodes that behave similarly with respect to the decomposition of different sources. These findings suggest that the activity of different electrodes does change over the course of an experiment in practice; treating these epoch recordings as realizations of an identical process could lead to misleading results. In summary, the proposed method underscores the importance of capturing the evolution in brain responses over the study period.

Constructing confidence intervals for the coefficients of high-dimensional sparse linear models remains a challenge, mainly because of the complicated limiting distributions of the widely used estimators, such as the lasso. Several methods have been developed for constructing such intervals. Bootstrap lasso+ols is notable for its technical simplicity, good interpretability, and performance that is comparable with that of other more complicated methods. However, bootstrap lasso+ols depends on the beta-min assumption, a theoretic criterion that is often violated in practice. Thus, we introduce a new method, called bootstrap lasso+partial ridge, to relax this assumption. Lasso+partial ridge is a two-stage estimator. First, the lasso is used to select features. Then, the partial ridge is used to refit the coefficients. Simulation results show that bootstrap lasso+partial ridge outperforms bootstrap lasso+ols when there exist small, but nonzero coefficients, a common situation that violates the beta-min assumption. For such coefficients, the confidence intervals constructed using bootstrap lasso+partial ridge have, on average, 50% larger coverage probabilities than those of bootstrap lasso+ols. Bootstrap lasso+partial ridge also has, on average, 35% shorter confidence interval lengths than those of the de-sparsified lasso methods, regardless of whether the linear models are misspecified. Additionally, we provide theoretical guarantees for bootstrap lasso+partial ridge under appropriate conditions, and implement it in the R package “HDCI”.

Feature selection is demanded in many modern scientific research problemsthat use high-dimensional data. A typical example is to find the most usefulgenes that are related to a certain disease (eg, cancer) from high-dimensionalgene expressions. The expressions of genes have grouping structures, forexample, a group of co-regulated genes that have similar biological functionstend to have similar expressions. Many statistical methods have been proposedto take the grouping structure into consideration in feature selection,including group LASSO, supervised group LASSO, and regression on grouprepresentatives. In this paper, we propose a fully Bayesian Robit regressionmethod with heavy-tailed (sparsity) priors (shortened by FBRHT) for selectingfeatures with grouping structure. The main features of FBRHT include that itdiscards more aggressively unrelated features than LASSO, and it can makefeature selection within groups automatically without a pre-specified groupingstructure. In this paper, we use simulated and real datasets to demonstratethat the predictive power of the sparse feature subsets selected by FBRHT arecomparable with other much larger feature subsets selected by LASSO, groupLASSO, supervised group LASSO, penalized logistic regression and random forest,and that the succinct feature subsets selected by FBRHT have significantlybetter predictive power than the feature subsets of the same size taken fromthe top features selected by the aforementioned methods.

Reconstructing pathogen dynamics from genetic data as they become available during an outbreak or epidemic represents an important statistical scenario in which observations arrive sequentially in time and one is interested in performing inference in an "online" fashion. Widely used Bayesian phylogenetic inference packages are not set up for this purpose, generally requiring one to recompute trees and evolutionary model parameters de novo when new data arrive. To accommodate increasing data flow in a Bayesian phylogenetic framework, we introduce a methodology to efficiently update the posterior distribution with newly available genetic data. Our procedure is implemented in the BEAST 1.10 software package, and relies on a distance-based measure to insert new taxa into the current estimate of the phylogeny and imputes plausible values for new model parameters to accommodate growing dimensionality. This augmentation creates informed starting values and re-uses optimally tuned transition kernels for posterior exploration of growing data sets, reducing the time necessary to converge to target posterior distributions. We apply our framework to data from the recent West African Ebola virus epidemic and demonstrate a considerable reduction in time required to obtain posterior estimates at different time points of the outbreak. Beyond epidemic monitoring, this framework easily finds other applications within the phylogenetics community, where changes in the data-in terms of alignment changes, sequence addition or removal-present common scenarios that can benefit from online inference.

This study introduces a new approach to power analysis in the context ofestimating a local average treatment effect (LATE), where the study subjectsexhibit noncompliance with treatment assignment. As a result of distributionalcomplications in the LATE context, compared to the simple ATE context, there iscurrently no standard method of power analysis for the LATE. Moreover, existingmethods and commonly used substitutes - which include instrumental variable(IV), intent-to-treat (ITT), and scaled ATE power analyses - require specifyinggenerally unknown variance terms and/or rely upon strong and unrealisticassumptions, thus providing unreliable guidance on the power of tests of theLATE. This study develops a new approach that uses standardized effect sizes toplace bounds on the power for the most commonly used estimator of the LATE, theWald IV estimator, whereby variance terms and distributional parameters neednot be specified nor assumed. Instead, in addition to the effect size, samplesize, and error tolerance parameters, the only other parameter that must bespecified by the researcher is the compliance rate. Additional conditions canalso be introduced to further narrow the bounds on the power calculation. Theresult is a generalized approach to power analysis in the LATE context that issimple to implement.

Experiments often include multiple treatments, with the primary goal to compare the causal effects of those treatments. This study focuses on comparing the causal anatomies of multiple treatments through the use of causal mediation analysis. It proposes a novel set of comparative causal mediation (CCM) estimands that compare the mediation effects of different treatments via a common mediator. Further, it derives the properties of a set of estimators for the CCM estimands and shows these estimators to be consistent (or conservative) under assumptions that do not require the absence of unobserved confounding of the mediator–outcome relationship, which is a strong and nonrefutable assumption that must typically be made for consistent estimation of individual causal mediation effects. To illustrate the method, the study presents an original application investigating whether and how the international legal status of a foreign policy commitment can increase the domestic political “audience costs” that democratic governments suffer for violating such a commitment. The results provide novel evidence that international legalization can enhance audience costs via multiple causal channels, including by amplifying the perceived immorality of violating the commitment.

Biclustering is a class of techniques that simultaneously clusters the rowsand columns of a matrix to sort heterogeneous data into homogeneous blocks.Although many algorithms have been proposed to find biclusters, existingmethods suffer from the pre-specification of the number of biclusters or placeconstraints on the model structure. To address these issues, we develop anovel, non-parametric probabilistic biclustering method based on Dirichletprocesses to identify biclusters with strong co-occurrence in both rows andcolumns. The proposed method utilizes dual Dirichlet process mixture models tolearn row and column clusters, with the number of resulting clusters determinedby the data rather than pre-specified. Probabilistic biclusters are identifiedby modeling the mutual dependence between the row and column clusters. We applyour method to two different applications, text mining and gene expressionanalysis, and demonstrate that our method improves bicluster extraction in manysettings compared to existing approaches.

Stochastic epidemic models (SEMs) fit to incidence data are critical toelucidating outbreak dynamics, shaping response strategies, and preparing forfuture epidemics. SEMs typically represent counts of individuals in discreteinfection states using Markov jump processes (MJPs), but are computationallychallenging as imperfect surveillance, lack of subject-level information, andtemporal coarseness of the data obscure the true epidemic. Analytic integrationover the latent epidemic process is impossible, and integration via Markovchain Monte Carlo (MCMC) is cumbersome due to the dimensionality anddiscreteness of the latent state space. Simulation-based computationalapproaches can address the intractability of the MJP likelihood, but arenumerically fragile and prohibitively expensive for complex models. A linearnoise approximation (LNA) that approximates the MJP transition density with aGaussian density has been explored for analyzing prevalence data inlarge-population settings, but requires modification for analyzing incidencecounts without assuming that the data are normally distributed. We demonstratehow to reparameterize SEMs to appropriately analyze incidence data, and foldthe LNA into a data augmentation MCMC framework that outperforms deterministicmethods, statistically, and simulation-based methods, computationally. Ourframework is computationally robust when the model dynamics are complex andapplies to a broad class of SEMs. We evaluate our method in simulations thatreflect Ebola, influenza, and SARS-CoV-2 dynamics, and apply our method tonational surveillance counts from the 2013--2015 West Africa Ebola outbreak.

We develop a multiscale spatial kernel convolution technique with higherorder functions to capture fine scale local features and lower order terms tocapture large scale features. To achieve parsimony, the coefficients in themultiscale kernel convolution model is assigned a new class of "Tree shrinkageprior" distributions. Tree shrinkage priors exert increasing shrinkage on thecoefficients as resolution grows so as to adapt to the necessary degree ofresolution at any sub-domain. Our proposed model has a number of significantfeatures over the existing multi-scale spatial models for big data. In contrastto the existing multiscale approaches, the proposed approach auto-tunes thedegree of resolution necessary to model a subregion in the domain, achievesscalability by suitable parallelization of local updating of parameters and isbuttressed by theoretical support. Excellent empirical performances areillustrated using several simulation experiments and a geostatistical analysisof the sea surface temperature data from the pacific ocean.

Markov chain Monte Carlo (MCMC) samplers are numerical methods for drawing samples from a given target probability distribution. We discuss one particular MCMC sampler, the MALA-within-Gibbs sampler, from the theoretical and practical perspectives. We first show that the acceptance ratio and step size of this sampler are independent of the overall problem dimension when (i) the target distribution has sparse conditional structure, and (ii) this structure is reflected in the partial updating strategy of MALA-within-Gibbs. If, in addition, the target density is blockwise log-concave, then the sampler's convergence rate is independent of dimension. From a practical perspective, we expect that MALA-within-Gibbs is useful for solving high-dimensional Bayesian inference problems where the posterior exhibits sparse conditional structure at least approximately. In this context, a partitioning of the state that correctly reflects the sparse conditional structure must be found, and we illustrate this process in two numerical examples. We also discuss trade-offs between the block size used for partial updating and computational requirements that may increase with the number of blocks.

Random forests are powerful non-parametric regression method but are severelylimited in their usage in the presence of randomly censored observations, andnaively applied can exhibit poor predictive performance due to the incurredbiases. Based on a local adaptive representation of random forests, we developits regression adjustment for randomly censored regression quantile models.Regression adjustment is based on new estimating equations that adapt tocensoring and lead to quantile score whenever the data do not exhibitcensoring. The proposed procedure named censored quantile regression forest,allows us to estimate quantiles of time-to-event without any parametricmodeling assumption. We establish its consistency under mild modelspecifications. Numerical studies showcase a clear advantage of the proposedprocedure.

Single-cell lineage tracking strategies enabled by recent experimental technologies have produced significant insights into cell fate decisions, but lack the quantitative framework necessary for rigorous statistical analysis of mechanistic models describing cell division and differentiation. In this paper, we develop such a framework with corresponding moment-based parameter estimation techniques for continuous-time, multi-type branching processes. Such processes provide a probabilistic model of how cells divide and differentiate, and we apply our method to study hematopoiesis, the mechanism of blood cell production. We derive closed-form expressions for higher moments in a general class of such models. These analytical results allow us to efficiently estimate parameters of much richer statistical models of hematopoiesis than those used in previous statistical studies. To our knowledge, the method provides the first rate inference procedure for fitting such models to time series data generated from cellular barcoding experiments. After validating the methodology in simulation studies, we apply our estimator to hematopoietic lineage tracking data from rhesus macaques. Our analysis provides a more complete understanding of cell fate decisions during hematopoiesis in nonhuman primates, which may be more relevant to human biology and clinical strategies than previous findings from murine studies. For example, in addition to previously estimated hematopoietic stem cell self-renewal rate, we are able to estimate fate decision probabilities and to compare structurally distinct models of hematopoiesis using cross validation. These estimates of fate decision probabilities and our model selection results should help biologists compare competing hypotheses about how progenitor cells differentiate. The methodology is transferrable to a large class of stochastic compartmental and multi-type branching models, commonly used in studies of cancer progression, epidemiology and many other fields.

In randomized trials, appropriately adjusting for baseline variables andshort-term outcomes can lead to increased precision and reduced sample size. Weexamine the impact of such adjustment in group sequential designs, i.e.,designs with preplanned interim analyses where enrollment may be stopped earlyfor efficacy or futility. We address the following questions: how muchprecision gain can be obtained by appropriately adjusting for baselinevariables and a short-term outcome? How is this precision gain impacted byfactors such as the proportion of pipeline participants (those who enrolled buthaven't yet had their primary outcomes measured) and treatment effectheterogeneity? What is the resulting impact on power and average sample size ina group sequential design? We derive an asymptotic formula that decomposes theoverall precision gain from adjusting for baseline variables and a short-termoutcome into contributions from factors mentioned above, for efficientestimators in the model that only assumes randomization and independentcensoring. We use our formula to approximate the precision gain from a targetedminimum loss-based estimator applied to data from a completed trial of a newsurgical intervention for stroke. Our formula implies that (for an efficientestimator) adjusting for a prognostic baseline variable leads to at least asmuch asymptotic precision gain as adjusting for an equally prognosticshort-term outcome. In many cases, such as our stroke trial application, theformer leads to substantially greater precision gains than the latter. In oursimulation study, we show how precision gains from adjustment can be convertedinto sample size reductions (even when there is no treatment effect).

Large-scale multiple testing with correlated and heavy-tailed data arises in a wide range of research areas from genomics, medical imaging to finance. Conventional methods for estimating the false discovery proportion (FDP) often ignore the effect of heavy-tailedness and the dependence structure among test statistics, and thus may lead to inefficient or even inconsistent estimation. Also, the commonly imposed joint normality assumption is arguably too stringent for many applications. To address these challenges, in this paper we propose a Factor-Adjusted Robust Multiple Testing (FarmTest) procedure for large-scale simultaneous inference with control of the false discovery proportion. We demonstrate that robust factor adjustments are extremely important in both controlling the FDP and improving the power. We identify general conditions under which the proposed method produces consistent estimate of the FDP. As a byproduct that is of independent interest, we establish an exponential-type deviation inequality for a robust U-type covariance estimator under the spectral norm. Extensive numerical experiments demonstrate the advantage of the proposed method over several state-of-the-art methods especially when the data are generated from heavy-tailed distributions. The proposed procedures are implemented in the R-package FarmTest.

We develop a model-based method for evaluating heterogeneity among several p × p covariance matrices in the large p, small n setting. This is done by assuming a spiked covariance model for each group and sharing information about the space spanned by the group-level eigenvectors. We use an empirical Bayes method to identify a low-dimensional subspace which explains variation across all groups and use an MCMC algorithm to estimate the posterior uncertainty of eigenvectors and eigenvalues on this subspace. The implementation and utility of our model is illustrated with analyses of high-dimensional multivariate gene expression.

We offer a survey of recent results on covariance estimation for heavy-taileddistributions. By unifying ideas scattered in the literature, we proposeuser-friendly methods that facilitate practical implementation. Specifically,we introduce element-wise and spectrum-wise truncation operators, as well astheir $M$-estimator counterparts, to robustify the sample covariance matrix.Different from the classical notion of robustness that is characterized by thebreakdown property, we focus on the tail robustness which is evidenced by theconnection between nonasymptotic deviation and confidence level. The keyobservation is that the estimators needs to adapt to the sample size,dimensionality of the data and the noise level to achieve optimal tradeoffbetween bias and robustness. Furthermore, to facilitate their practical use, wepropose data-driven procedures that automatically calibrate the tuningparameters. We demonstrate their applications to a series of structured modelsin high dimensions, including the bandable and low-rank covariance matrices andsparse precision matrices. Numerical studies lend strong support to theproposed methods.

Estimating causal effects for survival outcomes in the high-dimensionalsetting is an extremely important topic for many biomedical applications aswell as areas of social sciences. We propose a new orthogonal score method fortreatment effect estimation and inference that results in asymptotically validconfidence intervals assuming only good estimation properties of the hazardoutcome model and the conditional probability of treatment. This guaranteeallows us to provide valid inference for the conditional treatment effect underthe high-dimensional additive hazards model under considerably more generalitythan existing approaches. In addition, we develop a new Hazards Difference(HDi), estimator. We showcase that our approach has double-robustnessproperties in high dimensions: with cross-fitting, the HDi estimate isconsistent under a wide variety of treatment assignment models; the HDiestimate is also consistent when the hazards model is misspecified and insteadthe true data generating mechanism follows a partially linear additive hazardsmodel. We further develop a novel sparsity doubly robust result, where eitherthe outcome or the treatment model can be a fully dense high-dimensional model.We apply our methods to study the treatment effect of radical prostatectomyversus conservative management for prostate cancer patients using theSEER-Medicare Linked Data.

The central aim in this paper is to address variable selection questions in nonlinear and nonparametric regression. Motivated by statistical genetics, where nonlinear interactions are of particular interest, we introduce a novel and interpretable way to summarize the relative importance of predictor variables. Methodologically, we develop the "RelATive cEntrality" (RATE) measure to prioritize candidate genetic variants that are not just marginally important, but whose associations also stem from significant covarying relationships with other variants in the data. We illustrate RATE through Bayesian Gaussian process regression, but the methodological innovations apply to other "black box" methods. It is known that nonlinear models often exhibit greater predictive accuracy than linear models, particularly for phenotypes generated by complex genetic architectures. With detailed simulations and two real data association mapping studies, we show that applying RATE enables an explanation for this improved performance.

Advances in wearables and digital technology now make it possible to deliverbehavioral mobile health interventions to individuals in their everyday life.The micro-randomized trial (MRT) is increasingly used to provide data to informthe construction of these interventions. In an MRT, each individual isrepeatedly randomized among multiple intervention options, often hundreds oreven thousands of times, over the course of the trial. This work is motivatedby multiple MRTs that have been conducted, or are currently in the field, inwhich the primary outcome is a longitudinal binary outcome. The primary aim ofsuch MRTs is to examine whether a particular time-varying intervention has aneffect on the longitudinal binary outcome, often marginally over all but asmall subset of the individual's data. We propose the definition of causalexcursion effect that can be used in such primary aim analysis for MRTs withbinary outcomes. Under rather restrictive assumptions one can, based onexisting literature, derive a semiparametric, locally efficient estimator ofthe causal effect. We, starting from this estimator, develop an estimator thatcan be used as the basis of a primary aim analysis under more plausibleassumptions. Simulation studies are conducted to compare the estimators. Weillustrate the developed methods using data from the MRT, BariFit. In BariFit,the goal is to support weight maintenance for individuals who receivedbariatric surgery.

Predictability of behavior is an important characteristic in many fields including biology, medicine, marketing, and education. When a sequence of actions performed by an individual can be modeled as a stationary time-homogeneous Markov chain the predictability of the individual’s behavior can be quantified by the entropy rate of the process. This article compares three estimators of the entropy rate of finite Markov processes. The first two methods directly estimate the entropy rate through estimates of the transition matrix and stationary distribution of the process. The third method is related to the sliding-window Lempel–Ziv compression algorithm. The methods are compared via a simulation study and in the context of a study of interactions between mothers and their children.

Many popular methods for building confidence intervals on causal effectsunder high-dimensional confounding require strong "ultra-sparsity" assumptionsthat may be difficult to validate in practice. To alleviate this difficulty, wehere study a new method for average treatment effect estimation that yieldsasymptotically exact confidence intervals assuming that either the conditionalresponse surface or the conditional probability of treatment allows for anultra-sparse representation (but not necessarily both). This guarantee allowsus to provide valid inference for average treatment effect in high dimensionsunder considerably more generality than available baselines. In addition, weshowcase that our results are semi-parametrically efficient.

We consider an empirical likelihood framework for inference for a statisticalmodel based on an informative sampling design. Covariate information isincorporated both through the weights and the estimating equations. Theestimator is based on conditional weights. We show that under usual conditions,with population size increasing unbounded, the estimates are stronglyconsistent, asymptotically unbiased and normally distributed. Our frameworkprovides additional justification for inverse probability weighted scoreestimators in terms of conditional empirical likelihood. In doing so, itbridges the gap between design-based and model-based modes of inference insurvey sampling settings. We illustrate these ideas with an application to anelectoral survey.

Understanding of the effect of a particular treatment or a policy pertains tomany areas of interest -- ranging from political economics, marketing tohealth-care and personalized treatment studies. In this paper, we develop anon-parametric, model-free test for detecting the effects of treatment overtime that extends widely used Synthetic Control tests. The test is built oncounterfactual predictions arising from many learning algorithms. In theNeyman-Rubin potential outcome framework with possible carry-over effects, weshow that the proposed test is asymptotically consistent for stationary, betamixing processes. We do not assume that class of learners captures the correctmodel necessarily. We also discuss estimates of the average treatment effect,and we provide regret bounds on the predictive performance. To the best of ourknowledge, this is the first set of results that allow for example any RandomForest to be useful for provably valid statistical inference in the SyntheticControl setting. In experiments, we show that our Synthetic Learner issubstantially more powerful than classical methods based on Synthetic Controlor Difference-in-Differences, especially in the presence of non-linear outcomemodels.

Despite the availability of numerous statistical and machine learning toolsfor joint feature modeling, many scientists investigate features marginally,i.e., one feature at a time. This is partly due to training and convention butalso roots in scientists' strong interests in simple visualization andinterpretability. As such, marginal feature ranking for some predictive tasks,e.g., prediction of cancer driver genes, is widely practiced in the process ofscientific discoveries. In this work, we focus on marginal ranking for binaryprediction, the arguably most common predictive tasks. We argue that the mostwidely used marginal ranking criteria, including the Pearson correlation, thetwo-sample t test, and two-sample Wilcoxon rank-sum test, do not fully takefeature distributions and prediction objectives into account. To address thisgap in practice, we propose two ranking criteria corresponding to twoprediction objectives: the classical criterion (CC) and the Neyman-Pearsoncriterion (NPC), both of which use model-free nonparametric implementations toaccommodate diverse feature distributions. Theoretically, we show that underregularity conditions both criteria achieve sample-level ranking consistentwith their population-level counterpart with high probability. Moreover, NPC isrobust to sampling bias when the two class proportions in a sample deviate fromthose in the population. This property endows NPC good potential in biomedicalresearch where sampling bias is common. We demonstrate the use and relativeadvantages of CC and NPC in simulation and real data studies. Our model-freeobjective-based ranking idea is extendable to ranking feature subsets andgeneralizable to other prediction tasks and learning objectives.

There is growing interest in estimating and analyzing heterogeneous treatment effects in experimental and observational studies. We describe a number of metaalgorithms that can take advantage of any supervised learning or regression method in machine learning and statistics to estimate the conditional average treatment effect (CATE) function. Metaalgorithms build on base algorithms-such as random forests (RFs), Bayesian additive regression trees (BARTs), or neural networks-to estimate the CATE, a function that the base algorithms are not designed to estimate directly. We introduce a metaalgorithm, the X-learner, that is provably efficient when the number of units in one treatment group is much larger than in the other and can exploit structural properties of the CATE function. For example, if the CATE function is linear and the response functions in treatment and control are Lipschitz-continuous, the X-learner can still achieve the parametric rate under regularity conditions. We then introduce versions of the X-learner that use RF and BART as base learners. In extensive simulation studies, the X-learner performs favorably, although none of the metalearners is uniformly the best. In two persuasion field experiments from political science, we demonstrate how our X-learner can be used to target treatment regimes and to shed light on underlying mechanisms. A software package is provided that implements our methods.

We investigate how ideas from covariance localization in numerical weather prediction can be used in Markov chain Monte Carlo (MCMC) sampling of high-dimensional posterior distributions arising in Bayesian inverse problems. To localize an inverse problem is to enforce an anticipated “local” structure by (i) neglecting small off-diagonal elements of the prior precision and covariance matrices; and (ii) restricting the influence of observations to their neighborhood. For linear problems we can specify the conditions under which posterior moments of the localized problem are close to those of the original problem. We explain physical interpretations of our assumptions about local structure and discuss the notion of high dimensionality in local problems, which is different from the usual notion of high dimensionality in function space MCMC. The Gibbs sampler is a natural choice of MCMC algorithm for localized inverse problems and we demonstrate that its convergence rate is independent of dimension for localized linear problems. Nonlinear problems can also be tackled efficiently by localization and, as a simple illustration of these ideas, we present a localized Metropolis-within-Gibbs sampler. Several linear and nonlinear numerical examples illustrate localization in the context of MCMC samplers for inverse problems.

Mobile health is a rapidly developing field in which behavioral treatmentsare delivered to individuals via wearables or smartphones to facilitatehealth-related behavior change. Micro-randomized trials (MRT) are anexperimental design for developing mobile health interventions. In an MRT thetreatments are randomized numerous times for each individual over course of thetrial. Along with assessing treatment effects, behavioral scientists aim tounderstand between-person heterogeneity in the treatment effect. A naturalapproach is the familiar linear mixed model. However, directly applying linearmixed models is problematic because potential moderators of the treatmenteffect are frequently endogenous---that is, may depend on prior treatment. Wediscuss model interpretation and biases that arise in the absence of additionalassumptions when endogenous covariates are included in a linear mixed model. Inparticular, when there are endogenous covariates, the coefficients no longerhave the customary marginal interpretation. However, these coefficients stillhave a conditional-on-the-random-effect interpretation. We provide anadditional assumption that, if true, allows scientists to use standard softwareto fit linear mixed model with endogenous covariates, and person-specificpredictions of effects can be provided. As an illustration, we assess theeffect of activity suggestion in the HeartSteps MRT and analyze thebetween-person treatment effect heterogeneity.

Phylodynamics is a set of population genetics tools that aim atreconstructing demographic history of a population based on molecular sequencesof individuals sampled from the population of interest. One important task inphylodynamics is to estimate changes in (effective) population size. Whenapplied to infectious disease sequences such estimation of population sizetrajectories can provide information about changes in the number of infections.To model changes in the number of infected individuals, current phylodynamicmethods use non-parametric approaches, parametric approaches, and stochasticmodeling in conjunction with likelihood-free Bayesian methods. The first classof methods yields results that are hard-to-interpret epidemiologically. Thesecond class of methods provides estimates of important epidemiologicalparameters, such as infection and removal/recovery rates, but ignores variationin the dynamics of infectious disease spread. The third class of methods is themost advantageous statistically, but relies on computationally intensiveparticle filtering techniques that limits its applications. We propose aBayesian model that combines phylodynamic inference and stochastic epidemicmodels, and achieves computational tractability by using a linear noiseapproximation (LNA) --- a technique that allows us to approximate probabilitydensities of stochastic epidemic model trajectories. LNA opens the door forusing modern Markov chain Monte Carlo tools to approximate the joint posteriordistribution of the disease transmission parameters and of high dimensionalvectors describing unobserved changes in the stochastic epidemic modelcompartment sizes (e.g., numbers of infectious and susceptible individuals). Weapply our estimation technique to Ebola genealogies estimated using viralgenetic data from the 2014 epidemic in Sierra Leone and Liberia.

We provide a high-dimensional semi-supervised inference framework focused onthe mean and variance of the response. Our data are comprised of an extensiveset of observations regarding the covariate vectors and a much smaller set oflabeled observations where we observe both the response as well as thecovariates. We allow the size of the covariates to be much larger than thesample size and impose weak conditions on a statistical form of the data. Weprovide new estimators of the mean and variance of the response that extendsome of the recent results presented in low-dimensional models. In particular,at times we will not necessitate consistent estimation of the functional formof the data. Together with estimation of the population mean and variance, weprovide their asymptotic distribution and confidence intervals where weshowcase gains in efficiency compared to the sample mean and variance. Ourprocedure, with minor modifications, is then presented to make importantcontributions regarding inference about average treatment effects. We alsoinvestigate the robustness of estimation and coverage and showcase widespreadapplicability and generality of the proposed method.

Bayesian networks, with structure given by a directed acyclic graph (DAG), are a popular class of graphical models. However, learning Bayesian networks from discrete or categorical data is particularly challenging, due to the large parameter space and the difficulty in searching for a sparse structure. In this article, we develop a maximum penalized likelihood method to tackle this problem. Instead of the commonly used multinomial distribution, we model the conditional distribution of a node given its parents by multi-logit regression, in which an edge is parameterized by a set of coefficient vectors with dummy variables encoding the levels of a node. To obtain a sparse DAG, a group norm penalty is employed, and a blockwise coordinate descent algorithm is developed to maximize the penalized likelihood subject to the acyclicity constraint of a DAG. When interventional data are available, our method constructs a causal network, in which a directed edge represents a causal relation. We apply our method to various simulated and real data sets. The results show that our method is very competitive, compared to many existing methods, in DAG estimation from both interventional and high-dimensional observational data.