Papers

Many physical, biological, and epidemiological processes are governed by ordinary differential equations (ODEs) that are nonlinear in the state variable, including logistic population growth, chemical reaction kinetics, and epidemiological compartment models. We develop a differential equation-constrained local polynomial regression (DE-constrained LPR) framework for the general first-order ODE constraint g'(x) = F(x, g(x)), where F may be any Lipschitz continuous function, extending prior work restricted to exponential and linear ODE structures. Because F is generally nonlinear in g, the Taylor coefficients of the DE1-k estimator cannot be written in closed form; instead they are obtained by successive symbolic differentiation of F, and the estimator is computed by nonlinear least squares, requiring only a single local parameter at each evaluation point regardless of polynomial degree k. We derive the asymptotic conditional bias and variance of the DE1-k estimator, propose an AIMSE-optimal bandwidth that exploits the ODE structure to avoid direct estimation of high-order derivatives, and evaluate the method in a simulation study based on logistic growth, benchmarking against the parameter cascading method of Ramsay et al. (2007) (PCODE) and classical local linear regression. The DE-constrained estimator consistently outperforms local linear regression and is competitive with PCODE even though it estimates no structural parameter of the ODE; a sensitivity analysis across growth rates shows DE-constrained estimation becomes more accurate and more robust than PCODE as the curve steepens and PCODE's parameter estimation grows less stable. These results position DE-constrained LPR as a practical nonparametric alternative to parametric ODE-fitting methods when structural parameters are difficult to identify reliably.

Parameter identification from observations of dynamical systems is a fundamental problem in population biology. Mechanistic models of ecological systems rely on optimization methods that require accurate initial guesses to guarantee convergence. In ecological applications, datasets contain observation noise and are collected at sparse time points. This sparsity creates irregular likelihoods that cause standard optimization methods to struggle, while the ordinary differential equation solvers can become stiff or unstable in certain regions of the parameter space. These instabilities cause long running times or runtime errors. Here we present a computational framework for parameter identification that addresses these numerical instabilities by employing Natural Gradient Ascent, and we apply it to the classical Lotka-Volterra predator-prey model. We exploit the non-dimensionalization of the ordinary differential equations to treat scaling factors as nuisance parameters, reducing the dimensionality of the optimization problem. To prevent the solver step from becoming small, we implement an adaptive solver that switches between two independent second-order equations derived from the two components of the model. This approach allows Natural Gradient Ascent to converge in fewer iterations and with more stability than standard gradient ascent or BFGS methods. This framework provides a reliable method for parameter estimation in ecology when data is limited. The method can be generalized to other dynamical systems as long as the different components of the system do not become numerically problematic at the same time.

Selective testing based on values of auxiliary variables is an increasingly popular design strategy in observational studies and is ubiquitous in electronic health record data. Ignoring this underlying sampling mechanism can lead to biased estimation and erroneous scientific conclusions. Yet, rigorous analytic methods for accounting for two-phase sampling designs in observational settings remain under-utilized. Motivated by the Researching COVID to Enhance Recovery (RECOVER) Adult and Pediatric observational cohort studies, we describe common pitfalls and an approach for analysis of data collected via auxiliary variable dependent sampling.

Network experiments are used throughout the social and medical sciences to investigate causal effects under the presence of interference. While a large body of work has developed improved statistical procedures, the fundamental limits of statistical estimation in these settings is less well understood. In this paper, we develop and investigate a design-based theory of minimax risk for network experiments under an arbitrary neighborhood interference model. Our notion of minimax risk describes the optimal precision among all statistical procedures for investigating a particular causal effect on the observed interference network. We show that the minimax risk is a function of the corresponding conflict graph, which captures inherent unobservability of estimand-relevant potential outcomes given the observed interference network. Our main contribution is a series of upper and lower bounds on the minimax rate in terms of local and global connectivity properties of the conflict graph. To illustrate their utility, we apply these general results to obtain minimax analyses for two commonly studied effects: the direct treatment effect and global average treatment effect.

This paper develops procedures for nonparametric goodness-of-fit testing under covariate shift, where labelled data are drawn from a source population but goodness-of-fit is evaluated for a target population. The distribution mismatch is quantified by either a bounded moment condition or a sub-exponential tail condition on the target-to-source density ratio. Our method combines truncated importance-weighting kernel ridge regression with a multiplier bootstrap to construct confidence sets for the regression function. The truncation stabilizes the importance- weighting kernel ridge regression as well as the bootstrap calibration, making our approach applicable even when the density ratio has heavy tails. We prove nonasymptotic validity and sharpness of the resulting confidence sets under suitable operator compatibility conditions, and establish explicit error rates for coverage probability under specific conditions on the target- to-source density ratio and on the spectral decay of the kernel integral operator. Numerical experiments corroborate our theoretical findings.

In pan-omics pan-cancer studies, it is critical to identify latent sources of variation that are shared across particular subsets. This task often requires bidimensionally linked data matrices to be decomposed into a sum of block-sparse, low-rank modules. Existing approaches often rely on pre-specified module numbers, ranks, or post-hoc thresholding and can be sensitive to model specification when the underlying sharing structure is complex. To address these issues, we propose GL-BIDIFAC+, a group-regularized matrix factorization framework for discovering partially shared modules. It requires only an upper bound on the latent dimension and encourages module selection through group regularization with theoretically-motivated tuning parameter selection and local support recovery analysis, providing both scalability and principled guidance for module discovery. It also admits a probabilistic interpretation that enables model-based imputation of missing data. Simulation studies demonstrate accurate module recovery and favorable computational performance relative to existing approaches. We further apply GL-BIDIFAC+ to analyze the Cancer Genome Atlas data, where well-established molecular structure provides interpretable biological references. Our analysis distinguishes broad pan-cancer variation, cancer-specific subtype structure, and variation shared across cancers with related tissue origins or histologic features.

Advances in automation and high-throughput experimentation have enabled larger and more complex studies involving many factors and tests, creating a growing demand for computationally effective design construction methods. Efficient experimental design remains a key challenge in this context, creating a need for frameworks that can generate large experiments while preserving orthogonality and minimal aliasing. Unlike existing approaches which struggle with scalability, this work introduces an algorithmic framework for constructing large Orthogonal Minimally Aliased Response Surface (OMARS) designs by enumerating and combining weighing designs, three-level matrices with orthogonal columns and a fixed number of non-zero entries per column. Complete enumerations of weighing designs are achieved for designs with up to 24 tests, covering multiple numbers of factors and weights corresponding to two or three zeros per factor. In addition, a validated partial enumeration procedure and a combination method extend the catalog to substantially larger designs. The combination method enables the construction of OMARS designs for any test size that is a multiple of selected base sizes. This paper thus provides the methodology for generating large catalogs of high-quality OMARS designs, well-suited for high-dimensional screening and response-surface modelling in complex industrial and scientific experiments.

We propose the stochastic block Ornstein-Uhlenbeck (SBOU) process, a continuous-time multivariate model in which the drift matrix encodes a latent group structure among its components. Our main contribution is a community-detection algorithm whose misclassification proportion converges to zero in a regime combining infill, long-span, and high-dimensional asymptotics. To our knowledge, this is the first consistency result of this kind for latent group recovery in a discretely observed continuous-time multivariate model. As a key intermediate result, we establish consistency of the discretely observed maximum likelihood estimator of the drift matrix in the same regime, thereby extending the high-dimensional Lévy-driven Ornstein-Uhlenbeck literature. For practical implementation, we develop a feasible model-selection procedure for estimating the support of the drift matrix, which enables data-driven selection of the number of latent groups. The SBOU framework can be viewed as a continuous-time generalisation of the discrete-time stochastic-block VAR model, allowing for both positive and negative dynamic interactions between groups as opposed to only positive. We illustrate the methodology on the RE-Europe wind-capacity dataset and recover a country-level grouping consistent with the geographic benchmark.

Estimating the dynamic effects of economic shocks in short and very short samples is impeded by a lack of degrees of freedom. We offer a solution based on a Bayesian hierarchical framework for estimating local projection (LP) impulse response functions across a panel of related time series. The framework explicitly accommodates unbalanced panels in which some series are substantially shorter than others, allowing the short series to borrow information from longer ones at horizons where the short series carry little or no own data. Since series might exhibit heterogeneous dynamics, we develop a sparse finite mixture pool that clusters units by similarity of their impulse response profiles. We show in simulations that our approach substantially improves LP estimation accuracy relative to the standard approach if the time series are short while producing similar LPs for longer time series. Using a US price dataset, augmented with survey responses, we find that supply-chain and oil shocks trigger heterogeneous reactions of different price measures, with headline price indices responding more sharply than their core counterparts and goods prices changing more than services prices.

Formula One outcomes reflect the joint contributions of drivers and constructors, but these contributions are unobserved and vary over time. We propose a Bayesian state-space model that disentangles dynamic driver and constructor abilities using two observed outcomes: fastest qualifying lap times and race rankings. Both outcomes depend jointly on latent driver and constructor states that evolve at the Grand Prix level, while the race equation additionally accounts for starting-grid position. The decomposition is supported by constraints that center the driver and constructor abilities at zero, together with variation in driver-constructor assignments over time. Bayesian inference is performed using the No-U-Turn sampler under weakly informative priors that treat driver and constructor abilities symmetrically. Applying the model to the Formula One hybrid era from 2014 to 2021, we find substantial heterogeneity in both driver and constructor abilities. Driver abilities are generally more stable over time, whereas constructor abilities exhibit greater variation and, for many driver–constructor combinations, contribute more strongly to observed performance.

The mixture of experts (MoE) model provides a flexible approach for modeling heterogeneous regression relationships by allowing covariate-dependent mixing through a gating network, but most existing MoE models rely on parametric assumptions for expert error distributions, typically Gaussian, which can lead to inefficiency and sensitivity to outliers or heavy-tailed behavior when misspecified. We propose a semiparametric MoE model in which each expert error distribution is represented as a nonparametric Gaussian scale mixture estimated via nonparametric maximum likelihood, relaxing parametric assumptions within the Gaussian scale-mixture class while preserving the interpretability and structure of the MoE framework. The resulting model adapts to complex error structures, improves robustness under contamination and heavy tails, and remains competitive under well-specified Gaussian settings, providing a practical and theoretically grounded alternative to parametric MoE formulations.

On two-sided content platforms, symmetric two-sided isolation (assigning matched fractions of creators and viewers to isolated treatment and control submarkets) is widely used for creator-side and cold-start experiments because it removes cross-arm marketplace interference. Isolation, however, thins each viewer's candidate catalog, and intuition suggests the resulting engagement cost should fade as the platform grows: a small fraction of a vast catalog is still vast. We show that, in an order-statistics model of engagement, whether this intuition holds depends on the upper tail of match quality. Extreme-value theory yields tail-class loss laws with a sharp dichotomy: for light or bounded tails the loss vanishes as the candidate pool grows, whereas under heavy tails it converges to a size-independent constant, so expanding the candidate pool, even by orders of magnitude, does not asymptotically eliminate the cost. Evidence from two production experiments on a platform with millions of active creators is consistent with this picture: a pure A/A traffic sweep reveals a measurable, depth-graded engagement cost; a one-sided catalog ablation independently shows that per-viewer thinning contributes to the loss; and a tail index calibrated on the small exploration pool predicts an effect consistent with the one observed in the far larger full-catalog ablation. Isolation thus carries a price that experimenters should budget for, like any other cost. We give practitioners a preflight procedure that estimates it before launch, sizes traffic accordingly, and recommends a fallback design when the predicted cost exceeds a chosen tolerance.

On a non-collapsible survival mechanism, longitudinal-matching and target-trial-emulation methods are not competing estimators of one truth but answers to different causal questions, so a benchmark that scores them against a single "true hazard ratio" fabricates bias. We provide the direct comparison of relative efficiency, variance estimation, and model sensitivity that reviews find lacking. On a deliberately non-collapsible continuous-time Cox mechanism with known truth, the dominant families (sequential Cox, sequential stratification, risk-set matching, and inverse-probability-of-treatment-weighted (IPTW) marginal structural models) target numerically distinct causal estimands (marginal, conditional, two average-treatment-effect-on-the-treated, and intention-to-treat versus per-protocol). First, we quantify the phantom bias a shared marginal truth fabricates: 0.32-0.33 log-cumulative-hazard-ratio units for the matching estimators and 0.15 for the conditional method; the associational naive time-dependent Cox sits 0.76 away, a total discrepancy compounding the estimand gap with confounding. Second, a rank reversal: the recommended method flips with the target estimand, and a low-variance off-target estimator can still win on mean-squared error. Third, a cross-family variance result: the cluster-robust sandwich is closer to nominal for the trial-stacking estimator (0.90) but under-covers the matching estimators (0.77-0.82), which a prespecified n=500 bootstrap sub-study brings to 0.95-0.96. Fourth, model sensitivity: omitting a confounder induces 0.45-0.50 log-hazard-ratio bias and undercoverage, and intention-to-treat and per-protocol effects diverge as switching increases; a heart-transplant analysis illustrates these. On a second mechanism three of four findings replicate, the rank reversal attenuating and model sensitivity proving calibration-dependent.

Bayesian Go/No-Go decisions with co-primary endpoints require specifying prior distributions under the Normal-Inverse-Wishart framework; however guidance on how prior hyperparameters influence trial decisions remains limited. We propose a calibrated prior specification framework for bivariate Go/No-Go decisions. Skeptical and enthusiastic priors are calibrated so that each assigns a target probability to a clinically relevant decision region. We prove that for any prior precision $\kappa > 0$, a unique scale parameter $\lambda_0$ achieves the target calibration. Operating characteristics are evaluated across different $\kappa $ via simulation and applied to a phase~3 telitacicept lupus trial.The simulation result indicates $\kappa$ is the primary driver of prior discrimination. At $\kappa = 1$, the go rate difference between priors was 0.07; at $\kappa = 10$ it reached 0.56, with false positive rates below 0.01. Operating characteristics were robust to the degrees of freedom parameter $\nu_0$ and prior correlation $\rho_0$, supporting a default of $\nu_0 = 2$. In the lupus application, prior sensitivity was negligible at $\kappa = 1$ but at $\kappa = 10$ the enthusiastic go rate was three times the skeptical rate at small sample sizes. The framework reduces prior specification to two choices: the prior center and the prior precision $\kappa$. The identification of $\kappa$ as the dominant parameter, together with the cautious choice of $\kappa$ before the trial, motivates adaptive approaches to prior precision.

We comment on the optimally-transported generalized method of moments (OTGMM) estimator proposed by Schennach & Starck (2026a) and give counterexamples to Theorems 2-6 under their stated assumptions. First, the assumptions used in the small-error analysis are insufficient for consistency in Theorem 2 and asymptotic normality in Theorem 3. Next, we consider the large-error analysis, in which Theorem 4 states that the OTGMM estimator is equivalent to a GMM estimator with modified moments. We show that in a scalar model, Theorem 4 selects a value that differs from the unique OTGMM minimizer and violates the OTGMM sample moment restriction. In an overidentified model satisfying the assumptions used in Theorems 5 and 6, the first component of the Lagrange multiplier has different probability limits under the OTGMM estimator and the GMM estimator with modified moments. Under misspecification, the population value selected by OTGMM depends on the transport metric and on which variables may be adjusted. We give a sufficient condition under which solutions of the modified moment equations also solve the original constrained problem at a given parameter value, and separate conditions for consistency and asymptotic normality of the OTGMM estimator. We also show that Assumption 16 does not imply the matrix bound used in the supplemental proofs and replace Assumption 16 with a matrix condition that yields the bound.

Physics-informed and hybrid models integrate data and physics, and read learned coefficients as physical. Yet a small prediction error guarantees neither that the postulated operator is right, nor that the design identifies the coefficient, nor that training reached its minimum. On one parabolic benchmark where a sensor aggregates several diffusion rates and the truth is known, a correct operator recovers the coefficient to machine precision. A wrong one leaves it $30\%$ off while the in-domain error is 1.36 times the noise, and drifts $27\%$ across windows. At $7\%$ noise the error falls below the noise; the coefficient stays $30\%$ off. Without the truth, a post-training coefficient diagnostic questions the trained model with three classical statistics, without oracle; the first two need no refit. Does it explain the data? Can the design resolve the named physics? Does the coefficient depend on the window? The specification test keeps its pre-registered false-alarm ceiling and rejects every replicate of both wrong operators tested; a network frozen at its minimum returns the same verdicts. Where the design is blind, a wrong operator confined to the unseen rate leaves the test mute; the design reading flags it beforehand. A network stopped short of its minimum mimics a wrong operator; an analytic second fit separates the two. The test has a notch against a missing reaction term, miscalibrates under correlated noise and loses power beyond $7\%$ noise. Each warning names its remedy (operator, experiment or training); silence certifies nothing.

Longitudinal prediction of Alzheimer's disease biomarkers increasingly informs clinical decisions, and a forecast is only useful if it also reports how much to trust it. Conformal prediction supplies this by wrapping any forecaster in a prediction band with a finite-sample coverage guarantee under exchangeability. However, standard population-level conformal prediction guarantees only marginal coverage and may mask substantial under-coverage within clinically important subgroups. We introduce a general mechanism-driven framework for auditing and repairing such subgroup under-coverage. Across two cohorts (ADNI, OASIS-3), two base forecasters, and nine attributes spanning genetic risk, demographics, and clinical severity, we find that population-level bands under-cover high-risk subgroups in 57 of 68 audited combinations, despite achieving nominal marginal coverage. We trace these failures to two mechanisms: (A) rarity, where a group-conditional band calibrated on only $n$ patients covers at most $k/(n+1)$; and (B) tail-heaviness, where a population-wide band is too narrow for a heavy-tailed subgroup and additional data cannot close the gap. Under-coverage falls disproportionately on patients with high genetic risk and disease severity (6.1 pp mean deficit, 95% CI [3.3, 8.9]), while demographic groups remain at the target level on average (0.0 pp, CI [$-1.9$, 1.7]). We pair each mechanism with a corresponding conformal correction: cross-conformal pooling for rarity, per-subgroup calibration for tail-heaviness, and a coverage-safe marginal floor when both arise. Together, these corrections restore target coverage for nearly every high-risk subgroup across both cohorts and forecasters.

Objectives: While causal analysis of travel behavior is an emerging field, estimating heterogeneity in mode choice through causal modeling remains unexplored. This study demonstrates the application of a novel causal method, causal forest, to quantify the heterogeneity in travel mode choice shifts caused by the COVID-19 pandemic. Methods: We applied causal forests, a non-parametric causal machine learning method, to 802,935 trip records from the 2017 and 2022 waves of the National Household Travel Survey. The 2017 wave serves as the pre-pandemic control group, while the 2022 wave represents the treatment condition. Within the potential outcomes framework, we estimate average treatment effects (ATE), heterogeneous treatment effects (HTE), and conditional average treatment effects (CATE) across diverse socio-demographic groups and trip characteristics. Findings: Our results reveal an estimated ATE of a 1.86 percentage point (pp) increase in car-mode share, contrasted with decreases of 0.38 pp and 1.57 pp in public transit and walking, respectively. The largest increases in car use appeared for short-distance trips (one mile or less), households with annual incomes exceeding USD 200,000, and female travelers. Novelty: This is one of the first applications of causal forests to travel mode choice, and the first to use causal machine learning to estimate the pandemic's causal effect on mode choice analysis. Practical Applications: This study discusses methodological advantages, inherent assumptions, and limitations of causal forests within the context of transportation planning. This methodology is applied to COVID-19 travel data to illustrate how causal heterogeneity analysis can offer a deeper understanding of changes in mode choice. These insights are valuable for planners and policymakers in making policies related to mode shifts under an intervention.

Randomized controlled trials with similar protocols may yield conflicting findings when the distribution of relevant effect modifiers differs across study populations. Yet no formal statistical framework exists for defining and assessing whether conflicting trials are reconcilable, despite the importance of this question for evidence synthesis and regulatory decision making. To address this gap, we develop a causal inference framework for evaluating conditional and marginal reconcilability on additive and multiplicative scales in the presence of unmeasured effect modifiers. Within this framework, we use proxy variables for hypothesized unmeasured effect modifiers to develop regression-based tests of conditional reconcilability under parametric structural models. To assess marginal reconcilability, we extend existing transportability methods and develop an equivalence testing framework. We also introduce a reconciliation proportion to quantify the degree of marginal reconciliation. We illustrate these methods using the conflicting Meis and PROLONG trials of 17-alpha-hydroxyprogesterone caproate for preventing recurrent preterm birth. The analyses provided limited evidence that unmeasured effect modifiers such as cervical length, as captured by the selected proxies, were sufficient to marginally reconcile the trials. These findings demonstrate how proximal reconciliation methods may help regulators, researchers, and clinicians evaluate whether differences in study populations explain conflicting trial findings.

We propose a resolution-adaptive Bayesian wavelet denoising method for noisy one-dimensional signals. Its central innovation is a spike-and-slab prior whose continuous slab mixes a compactly supported Wendland-type polynomial density with the more dispersed semicircle density. A low-dimensional empirical-Bayes trend produces data-adaptive mixture weights by resolution, while a data-adaptive support scale controls the common bounded interval. Thus, the method combines sparsity, explicit support control, and interpretable resolution-dependent shrinkage. Under squared-error loss, we derive the posterior-mean estimator and establish symmetry, boundedness, continuity, and limiting properties. We define fixed-hyperparameter bias, variance, and risk and develop an empirical-Bayes fitting procedure. Under a Laplace working likelihood, the Wendland contribution has finite-sum expressions, while the semicircle contribution is evaluated by stable one-dimensional integration. Simulations with the Bumps, Blocks, Doppler, and HeaviSine signals compare Gaussian- and Laplace-likelihood versions with universal thresholding, false-discovery-rate (FDR) thresholding, cross-validation (CV), Stein's unbiased risk estimate (SURE), the Bayesian adaptive multiresolution shrinker (BAMS), and a nonlocal-prior (NLP) method. In the primary Gaussian-error study, WS–Gaussian was the best non-NLP method in 24 of 36 cells, including 11 of 12 low-SNR cells, with a much more favorable computational profile than WS–Laplace. A real seismic acceleration record from the 2008 Chino Hills earthquake illustrates attenuation of rapid fluctuations and preservation of the dominant event. A semi-synthetic study using the processed trace as surrogate truth showed improvement over the noisy observation at lower and moderate SNRs, but not at the highest SNR.

Anytime-valid inference enables analysts to continuously monitor their data and stop experiments early. However, the majority of these methods incur a certain conservativeness by remaining valid on infinite time horizons. In practice, a bound on the horizon may be imposed due to budgetary, practical, or ethical constraints. In this paper, we ask the question: "Is it possible to obtain sharper large-sample anytime-valid inference by forgoing validity beyond some finite time horizon?". We provide a positive answer to this question by proposing a family of statistical objects that we call "confidence horizons". These objects can be viewed as large-sample confidence sequences on bounded time horizons, or alternatively as group sequential repeated confidence intervals with a maximal number of interim peeking times. We make explicit connections to the group sequential boundaries of Pocock [1977], O'Brien–Fleming [1979], and Wang–Tsiatis [1987]. We derive closed-form distribution functions of certain statistics which can be used to calculate the asymptotic quantiles of confidence horizons exactly, sidestepping the repeated integration typically employed in group sequential methods. We illustrate the use of confidence horizons for treatment effect estimation in sequentially randomized experiments under adaptive Neyman allocation.

Empirical work often removes fixed effects, latent factors, or high-dimensional controls before estimating structural relationships. These transformations reduce confounding but may also remove identifying variation. We study linear panel IV after one equation-compatible nuisance projection under two-way dependence. The projected Jacobian determines which structural directions remain visible; the projected-score law determines their precision; and, on Gaussian fixed-rank strata, they combine in a Projected Information Matrix. We derive weak-identification limits with dimension-specific information accumulation, feasible factor-transfer conditions, identification-robust tests, bootstrap procedures for non-Gaussian interaction limits, and inference for the projected spectrum, rank, subspaces, and information matrix. Simulations show that a raw first-stage statistic above 500 can support the wrong sign while projected diagnostics reveal weak valid information. In an international monetary application, common projection substantially attenuates apparent foreign-output persistence, while Gaussian-reference Anderson–Rubin sets remain unbounded. Identification should therefore be assessed after nuisance removal.

Dimensionality reduction and effect estimation are central tasks in the analysis of high-dimensional data such as in genomics. Generalized bilinear models (GBMs) provide a versatile framework for these tasks, however, the statistical and computational efficiency of GBMs degrades rapidly as the number of sample covariates grows. To address this limitation, we introduce reduced-rank generalized bilinear models (RR-GBMs), which employ a reduced-rank sample coefficient matrix to model the effect of a large number of covariates without requiring an excessive number of parameters. In simulation studies, we find that when the true sample coefficient matrix is reduced rank or close to reduced rank, RR-GBM outperforms the standard full-rank GBM both statistically and computationally, providing more accurate estimates with lower computational burden. We develop a data thinning approach for model selection in the RR-GBM framework, facilitating rank selection. Furthermore, RR-GBM enables a new approach to visualizing the relationships among covariates and among features. We demonstrate the method in an application to Perturb-seq data for pancreatic cancer.

Safety surveillance increasingly involves repeated monitoring of many exposure-outcome signals in observational healthcare data, where sparse information, dependence across related signals, and systematic error can complicate inference. Existing frameworks typically focus on either correcting residual bias using negative controls or borrowing information across exposure-outcome pairs, but not both. We propose a multi-signal Bayesian sequential surveillance framework that integrates empirical bias correction with low-rank latent factor modeling. At each analysis time, a hierarchical Bayesian model learns exposure-specific bias distributions from negative control outcomes assumed to have null latent effects. Conditional on these distributions, low-rank latent factors are estimated across exposures and outcomes of interest to share information across correlated signals. As new data accrue, posterior inference is updated sequentially, yielding bias-corrected posterior summaries of effect sizes across multiple monitored signals. We illustrate the method in a postmarket vaccine safety surveillance study using a large US insurance claims database.

Missing values present a common challenge in statistical modeling, so handling them properly is an important research direction. Among the various mechanisms that can generate missing values, the most common is the missing-at-random (MAR) mechanism, in which the probability of missingness depends only on observed data and not on unobserved data. This paper addresses the problem of estimating a multivariate linear regression model with multiple random covariates in the presence of MAR values in both the response and covariate spaces using a maximum likelihood (ML) framework. The proposed methodology models the joint distribution of responses and covariates through a conditional-marginal factorization of a multivariate Gaussian distribution. This formulation can be interpreted as a reparameterization of the multivariate normal distribution when the variables can be naturally partitioned into responses and covariates. Parameter estimation is performed using the expectation-maximization (EM) algorithm, which facilitates the imputation of missing values while preserving the distinct roles of responses and covariates. We extend this framework to the model-based clustering setting by considering a mixture of multivariate linear regressions with multiple random covariates. This extension enables soft clustering under incomplete data and accommodates MAR values in both the multivariate responses and covariates. Hence, it represents one of the most general model-based clustering solutions for regression data currently available in the literature. The effectiveness of the methodology is demonstrated through a simulation study, and the advantages of the proposed reparameterization are illustrated using the Automobile dataset, which contains missing values.

Many prediction problems seek to infer an unobserved, instance-specific parameter that governs the distribution of an observable response, even though the latent parameter is unavailable for both historical and future instances. We develop LatentCP, a prior-free conformal framework that constructs uncertainty sets for latent distributional parameters using only observed context–response pairs and a specified forward model. The method first constructs a conformal prediction set in the observable response space and then retains candidate latent parameters according to the probability their induced response distributions assign to that set. This inversion provides finite-sample marginal coverage without requiring latent calibration labels, a unique inverse mapping, or knowledge of the latent mixing distribution. Because latent-set efficiency depends nonmonotonically on the response-space miscoverage rate, we further introduce a multilevel procedure that aggregates normalized incompatibility scores across several response sets and selects the aggregation distribution using an independent tuning sample. Across synthetic experiments, LatentCP maintains nominal latent coverage under weak forward identification, observational nonidentifiability, and latent heterogeneity and multimodality, where empirical-Bayes, likelihood-based, and proxy-label conformal methods can substantially under-cover. On a California wildfire real dataset, it produces spatially adaptive uncertainty sets for latent fire intensity. Independent tuning improves efficiency, while multilevel aggregation provides additional gains when different response levels contain complementary information, without sacrificing validity.

The drift diffusion model (DDM) is a cornerstone of cognitive decision-making research. Although numerous estimation methods exist, researchers continue to seek inference approaches that are both fast and flexible across diverse study designs. Amortized Bayesian inference (ABI) can provide nearly instantaneous inference for complex stochastic models like the DDM, but neural networks trained for one study design cannot generalize to others. In this paper, we propose a divide-and-conquer framework that address this limitation. The core idea is that the DDM's independence assumption allows the full dataset to be decomposed into pairwise shards, each sharing a common structure that a single neural network can learn. Inference is performed on each shard separately and the resulting posteriors are combined via consensus MCMC to approximate the full posterior. Using simulated datasets, we evaluate the accuracy and uncertainty of this method. Our results show that the proposed divide-and-conquer approach achieves accuracy and uncertainty comparable to MCMC while reducing computational cost by several orders of magnitude. This work not only advances DDM estimation but also demonstrates a general strategy for improving the scalability and generalizability of ABI methods across diverse applications.

We develop a fully nonlinear structural vector autoregressive framework in which the contemporaneous structural mapping may be nonlinear and non-additive. Identification is achieved by exploiting variation in the conditional distributions of the mutually independent structural shocks induced by an observed exogenous variable. Specifically, a general contrastive learning framework that makes use of this variation together with the assumed exponential-family structure is employed to recover the shocks. Existing independent innovation analysis results identify such shocks only up to arbitrary componentwise invertible transformations, which is generally insufficient for structural econometric analysis. We strengthen this result by imposing a structured exponential-family specification for the conditional shock distributions. With the imposed sufficient statistics, the remaining ambiguity is reduced to a one-parameter transformed-scale map for each shock. We then show that, under a logistic specification used for the natural parameters, the identification is further strengthened up to permutation and componentwise sign changes. Once the shocks have been recovered, the fully nonlinear structural vector autoregression can be estimated using feed-forward neural networks, motivated by their universal approximation capabilities. The empirical application studies asymmetries in the responses of U.S. industrial production to the real oil price shock. We find modest asymmetries with respect to the sign of the shock and state of the economy. The accompanying R package iiasvar implements the introduced methods.

The effect size and nuisance parameters needed to appropriately size a clinical trial are generally not adequately understood before a trial begins. Through pre-planned early stopping rules, a conservatively planned group sequential design can effectively adapt the sample size for a clinical trial for a new treatment that is ineffective, very effective or minimally effective. One potential issue with such designs is that a substantial number of patients can be enrolled after a data cutoff for an interim analysis while data are being entered, cleaned, analyzed and discussed. A strategy of re-estimating the sample size at an interim analysis based on conditional power has been proposed to reduce somewhat this enrollment overrun issue. A purported advantage sometimes claimed is a smaller up-front planned sample size than a conservatively planned group sequential design. We demonstrate derivation of 2-stage group sequential and conditional power designs using the gsDesign R package and suggest comparing designs with comparable power using expected sample size calculations. While there are cases where conditional power designs may have small advantages, it is quite easy to derive very inefficient conditional power designs. This, along with the fact that a conditional power design may reveal something about the interim treatment effect, will often leave a group sequential design as the design of choice.

Scalable spatiotemporal modeling remains challenging because conventional methods rely on covariance models that tightly couple spatial representation and temporal inference, often leading to high computational costs. To address this difficulty, we developed coarse-to-fine spatiotemporal modeling (CF-STM), a framework that extends coarse-to-fine spatial modeling (CF-SM) to spatiotemporal settings. CF-STM represents latent spatial processes through multiscale locally weighted models, with temporal dependence modeled separately through state-space models defined at local centers. This covariance-free formulation achieves scalable computation while allowing temporal models to be flexibly specified without altering the spatial representation. Monte Carlo experiments demonstrate predictive performance comparable to that of alternative scalable space-time models at a substantially lower computational cost. An application to long-term residential land price data in the Tokyo metropolitan area shows that CF-STM flexibly captures complex spatiotemporal patterns while enabling interpretable inferences. CF-STM is implemented in an R package spCF (https://cran.r-project.org/web/packages/spCF/).