Papers

A common question about networked time series is whether outcomes changed because shocks transmit more strongly or because the pattern of connections changed. Standard practice inserts a recorded network into an outcome regression and reads movements of the fitted coefficient as changes in transmission strength. When the network is latent, time varying, and measured with error, this reading fails: changes in strength and changes in composition can produce the same outcome distribution at a single date, and the population coefficient moves under composition changes alone. The question becomes answerable when outcomes are analyzed jointly with repeated noisy measurements of the network, such as paired reports of bilateral trade flows. For the joint model we establish necessary and sufficient conditions for local identification, estimators of the strength and composition paths, a simultaneous confidence band for the strength path, confidence sets that remain exact under weak identification, breakdown bounds under common reporting bias, and an exactly sized test that detects changes on the observed path and attributes them to strength or to composition. Simulations assess each procedure at its stated boundary. On a mirror-reported trade panel of eighteen economies over 1995 to 2020, the diagnostics flag exactly the crisis years and the composition coordinate attached to European Union membership declines by roughly two thirds. The estimand is predictive dependence, not a causal effect.

Assessing fine-scale spatio-temporal air pollution contrasts in urban contexts is a major challenge for environmental epidemiology. We propose a Bayesian spatio-temporal model to predict daily concentrations of NO$_2$, PM$_{10}$, and PM$_{2.5}$ in Rome (Italy), over 2011-2022, on a fine grid scale (1 km), using data from 8 to 13 monitoring stations, depending on pollutant. The model includes meteorological and temporal (working day/weekend) fixed effects together with a lag-1 autoregressive spatio-temporal random effect aimed at capturing spatial and daily dependence. Predictive performance was assessed by leave-one-site-out cross-validation. Estimated exposures were then linked to geolocated cause-specific mortality data for Rome (2012-2019), and a case-crossover time-stratified approach was adopted to investigate acute effects. Cross-validation showed good overall predictive performance, with larger errors at high-traffic sites. Estimated exposure (mean lag 0-5) was positively associated with natural-cause mortality, with percent increase in risk of 1.3 (95% CI: 0.6-2.0) for PM$_{10}$, 2.1 for NO$_2$ (1.4-2.9), and 2.4 for PM$_{2.5}$ (1.4-3.3) per 10 $\mu$g/m$^3$ increase. Corresponding estimates when using the city-specific daily average exposure, instead of our 1 km resolution model, were of comparable magnitude. The proposed Bayesian spatio-temporal framework provides reliable fine-scale exposure estimates for epidemiological use, with results consistent across independent exposure estimates, supporting its application in urban air-pollution health studies.

Global forecasting models pool many series and learn one shared function. Gradient-boosted trees are their most common form. We measure a failure of this design that has not, to our knowledge, been documented. Train a global tree on the individual series of a hierarchy, then ask it for the hierarchical aggregate. The aggregate sits far outside the model's training range, and the forecast collapses. The model under-predicts the total by 30-50x in our production deployment, and by up to 496x in a public M5 reconstruction. The mechanism is known: beyond its training range, a tree predicts a constant. It surfaces at the aggregate because the total dwarfs every training series. The cure is not new. Per-series scaling, the preprocessing step that Montero-Manso and Hyndman (2021) recommend, prevents the collapse. So do a weighted aggregate-level training row and seasonal differencing. Our contribution is the characterization. The collapse reproduces on five panels: a production business-to-business marketplace, a synthetic hierarchy, M5, Australian Tourism, and a public business-buyer panel. It holds on three tree libraries, is invariant across training seeds, and is statistically significant. Its onset is immediate and tracks a simple support bound: a scale gap of only 1.15x already costs a third of the total. No standard configuration change prevents it: pooling every hierarchy level into training fails at scale, and the one knob that fits linear models in the leaves softens it without curing it. Rolling the forecasts forward recursively separates the cures: the aggregate-row cure re-collapses, per-series scaling degrades but stays low, and only seasonal differencing keeps its one-step accuracy unchanged. We close with a three-step procedure for diagnosing and preventing the failure in deployed systems.

Autism spectrum disorder is a heterogeneous neurodevelopmental condition whose functional brain organization varies across individuals and imaging centers. Resting-state functional connectivity provides an opportunity to study this variation, but analysis is complicated by high dimensionality, strong dependence among connections, and substantial site heterogeneity. We introduce a Hybrid Spatial Statistical Learning Framework for individualized probability estimation in multisite autism neuroimaging. The framework combines edge-level connectivity, graph-based network organization, and participant characteristics as complementary representations of a common probability target. Fine-scale connectivity retains discriminatory information, while broader representations stabilize probability estimates under site shift. All preprocessing and model development are performed within complete site-held-out validation to assess transportability to unseen acquisition environments. Simulations show that the hybrid preserves the discrimination of the strongest edge model while improving probability accuracy as between-site heterogeneity increases. In 860 participants from 20 ABIDE sites, the hybrid retained ROC AUC while reducing Brier score and log loss. These results show that multiscale integration can improve the reliability and transportability of individualized probabilities from complex, spatially dependent biomedical data.

Functional decline in older adults is typically recognized only after falls or observable gait impairment, closing the window for prevention. We investigated whether temporal trajectories of routine biomarkers, already recorded but rarely analyzed longitudinally, can identify patients in the pre-clinical phase of mobility decline. Using the All of Us Research Program (N = 297,861; 11.1% cases), we derived trajectory features (slope, variability, delta, mean) for twelve biomarkers over a three-year pre-index window. LightGBM models incorporating trajectories significantly outperformed static laboratory summaries (AUROC 0.797 vs. 0.755; DeLong p < 0.001; AUPRC 0.380 vs. 0.304). A 1:1 age- and sex-matched analysis confirmed an independent trajectory signal (AUROC 0.727 vs. demographics-only 0.680). A horizon analysis demonstrated sustained prediction 3-12 months before decline onset (AUROC 0.768-0.740). Because the model uses only measurements already ordered in routine care, it supports passive, zero-burden EHR integration for early detection of pre-clinical functional decline.

We introduce a novel outlier detection method that utilizes a concrete measure of the statistical impact of an outlier rather than a vague notion of "unusualness". This allows practitioners to only flag points that materially affect the inferences made from the data while also ensuring these points are anomalous. This is of particular interest in experimentation, where outliers affect the validity of the causal inferences. In particular, we consider a potential outlier's impact on the mean and the False Positive Rate (FPR) of a hypothesis test.

Emile Durkheim's Suicide: A Study in Sociology (1897) predates much of the statistical machinery that quantitative modellers now take for granted. Yet, working with sparse and imperfect observational data, Durkheim repeatedly arrives at practices that remain remarkably relevant to modern modelling. This paper revisits Suicide from the perspective of quantitative risk modelling, not for its substantive conclusions, but for the reasoning by which Durkheim reached them. His approach illustrates how careful definition, common sense, logical investigation, scepticism toward convenient explanations, and close attention to what the data can and cannot support—all of which must precede, and can often substitute for, statistical sophistication. The broader lesson is simple: good modelling begins not with technique, but with understanding the problem, interrogating the evidence, and reasoning carefully about what it can actually tell us.

We study the news cascade triggered by a single high-impact event as a pair of coupled stochastic processes: the daily article count $N(t)$ and the reliable-vs-unreliable composition $\pi_{\mathrm{pos}}(t)$. Using the 2017 Las Vegas shooting as a primary case study (CC-News + NewsGuard, 91 days, 2,702 articles), we fit and compare four counting-process models for $N(t)$ and four distributional models for $\pi_{\mathrm{pos}}(t)$ on a common likelihood scale. For counts, a hybrid Inhomogeneous-Poisson plus Hawkes model wins decisively on AIC, BIC, and log-likelihood, as the IHP accounts for the initial exogenous shock, while the Hawkes accounts for the subsequent self-excitation. A profile-likelihood sweep of the Hawkes decay rate collapses the optimal Hawkes kernel to an AR(1) one-step lag $\lambda(t)=\mu+nN(t-1)$. This Markovian structure on the counts directly motivates the Part-B reliability models: we treat the reliability label as a first-order Markov chain on a binary state space and test homogeneous, regime-switching, and mean-field parameterizations against a bivariate Hawkes baseline model using a shared conditional binomial likelihood. The results show that the bivariate Hawkes model and the non-homogeneous Markov models are nearly indistinguishable on fit, as they exhaust the signal and start fitting to the noise. A brief replication on the 2017 Hurricane Harvey data shows that the same methodology applies but produces different results due to Hurricane Harvey's gradual endogenous build-up.

Flood-monitoring sensors are often placed according to coverage, access, or expected inundation. However, the value of a measurement depends on the prediction or decision it is intended to inform. Using tRIBS-Urban simulations and a neural-network surrogate of the August 2014 metropolitan Detroit flood, we examine how this learning target changes single-sensor placement. Across 2,576 candidate locations, we compare parameter-oriented optimal experimental design (PO-OED), which values expected information gain (EIG) about model parameters, with goal-oriented optimal experimental design (GO-OED), which values EIG about specified flood predictions. We also examine how parameter EIG evolves during the event, and illustrate that parameter learning translates unevenly into reductions in predictive uncertainty across locations and lead times. Under GO-OED, point-depth targets favor nearby locations, whereas regional-average and regional maximum-depth targets can favor nonlocal locations. Weighted multi-point objectives retain similar broad spatial patterns, although their computed max-EIG locations differ. Public geospatial data further provide illustrative feasibility and contextual classifications for deployment screening. These results show how monitoring objectives shape sensor placement in optimal experimental design, and motivates an objective-first workflow that defines the intended prediction and priorities, applies field-verified restrictions, and ranks locations by EIG.

Polymeric materials are widely used in photovoltaic (PV) systems, making it essential to understand their service life to ensure reliable PV performance. The primary failure mechanism of polymeric materials in PV systems is photodegradation caused by ultraviolet (UV) radiation. Degradation modeling provides a framework for predicting service life, with a key step being the development of predictive models for degradation paths. This paper presents statistical and machine learning approaches for predicting the outdoor degradation of polymeric components in PV systems. We describe the study design and data collection process for developing predictive models based on indoor laboratory testing data, which are then extended to outdoor field conditions with time-varying environmental variables, with prediction uncertainty quantified through simulation. Deep learning (DL) methods are also explored, and results are compared across modeling approaches. The parametric statistical model demonstrates good fit and predictive performance across datasets and shows greater robustness by incorporating physical and chemical knowledge. The DL model provides flexibility in capturing complex covariate relationships and often yields accurate predictions, though it is less robust across datasets. The paper concludes with remarks on key findings and their implications for PV reliability.

Marketing Mix Models (MMMs) estimate the incremental sales effect of advertising from observational time series, yet they are rarely validated against ground truth, because ground truth is unobservable in real data. Synthetic data closes that gap in principle, but existing generators produce marketing spend exogenously - omitting the central difficulty of the estimation problem, since real budgets are planned around promotional calendars, seasons, and recent performance. This paper presents a parameterized generator, and a fixed reference instance, of a synthetic weekly retail dataset (156 weeks, three media channels) in which spend arises from four documented coordination mechanisms - quarterly budget feedback, anticipatory spending ahead of a promotional calendar, scheduled TV bursts, and algorithmic performance chasing - on a demand baseline with seasonal, quality, price, and unobserved sentiment components. Spend translates into incremental sales through two transformations, carryover and diminishing returns, instantiated here as geometric adstock and logistic saturation with known parameters; the true causal decomposition of every week's sales is recorded alongside the error-contaminated variables a practitioner would observe. Every mechanism is a parameter that can be varied or switched off, and a companion procedure simulates go-dark geo-experiments with exact treatment effects. The seeded generator and reference instance are publicly released with notebooks that reproduce every number in this paper.

Marketing Mix Models (MMMs) are widely used for marketing measurement and budget allocation, but face fundamental identification challenges: due to endogenous marketing spend decisions, MMM estimation on observational time-series data cannot recover the true causal effects of marketing. On the other hand, geo-experiments provide causal identification through randomization, but it is not clear how to use them efficiently to calibrate marketing mix models. We propose a novel structural estimation approach that recovers the complete set of MMM parameters - adstock decay ($\alpha$), saturation ($\lambda$), and effectiveness ($\beta$) - directly from geo-experimental time-series. By differencing outcomes between treatment and control regions, our method eliminates observed and unobserved confounding factors while preserving the temporal variation that identifies each parameter. We demonstrate on synthetic data that this approach recovers the true ROAS and the response curve over the range of spending covered by the experiments, together with credible estimates of the underlying parameters. Our framework enables efficient pooling across multiple experiments and provides a principled foundation for MMM calibration that fully utilizes the information experiments contain.

Members of collective institutions may leave public traces while their individual choices remain concealed. This paper separates a corpus-conditional public position from the behavioral rule linking that position to participation and secret choice. I measure the first with a dynamic ordinal state-space model and establish a likelihood-exclusion result for the second. If behavioral-link parameters enter only the distribution of an entirely unobserved outcome and are a priori independent of measurement parameters, public signals leave their posterior equal to their prior. Declared monotone mappings are therefore reported as sensitivity envelopes, not estimates or identified bounds. Cabinet-sized simulations compare dynamic measurement with a static ordinal model and transparent summaries, exposing gains from temporal pooling and failures under weak measurement or omitted dimensionality. The application reconstructs 41 pre-decision official-source records preceding a concealed Slovenian cabinet decision. Only six of 21 ministers have attributable signals. The model locates a small visible group but leaves the remainder prior-dominated. Event-dependence adjustments, coding perturbations, rolling prediction, and the subsequently disclosed record support the same conclusion. Uncertainty should track the public information environment rather than be converted into unsupported behavioral precision.

Integrating summary-level scientific knowledge into neural network models provides a practical strategy for transferring prediction models trained on adequately sampled source cohorts to underrepresented target populations, where individual-level data in the target domain are often limited or unavailable. In this study, we propose transfer prediction strategies incorporating external summary-level scientific knowledge and illustrate its application on the PRISMA Maternal and Neonatal Health Study, training a neural network model on the source data to predict adverse outcomes in the target cohorts. Besides, we also extend the existing GRU-D framework by incorporating static feature embeddings and attention weights to jointly leverage temporal and static information for improved prediction. Our approach employs soft labels derived from summary-level statistics describing the target population to fine-tune GRU-D-Static models that are initially trained on the source populations. We evaluate six maternal and neonatal outcomes, including stillbirth, preterm birth, low birth weight, small vulnerable newborn, neonatal death, and maternal near miss. Across all tested scenarios, fine-tuning using soft labels from just basic covariates substantially improved predictive performance compared with deep learning models trained on the source sample. Furthermore, the performance slightly improves more when additional covariates were incorporated into the logistic regression model or when partial input features from the target population were available for fine-tuning. These findings demonstrate that integrating existing scientific knowledge in the literature through transfer prediction of source neural network models can enhance prediction performance in underrepresented target populations, reducing reliance on large-scale data collection and supporting risk prediction in global health.

Clinical trial design relies on both statistical and clinical considerations for pre-specification of potentially practice-changing target treatment effects. As larger trials tend to be associated with high power and modest minimal detectable benefit, trial sample size is typically calibrated with reference to relevant precedents to prevent overpowering. Albeit trial sponsors and regulators are accustomed to this practice, there is scope for simplification to enhance the robustness, transparency and cross-trial consistency of this aspect of the design process. To this end, a design utility index is proposed here as a formal basis for sample size calibration, balancing the increase in power at higher sample sizes against the corresponding reduction in the magnitude of minimum detectable treatment effects, without requiring additional statistical assumptions or bespoke software. Application of utility calibration to a broad range of designs demonstrates consistency with regulatory expectations for minimum power, particularly in confirmatory settings, and effective protection against overpowering and against overly aggressive interim analyses.

A forecasting benchmark reports which method won. We show that the answer is set by the evaluator's choices before any model is fitted. We benchmark 24 forecasting methods and one textbook reference, including six 2025-era time series foundation models, on a production marketplace panel of 1,887 business customers over 67 months. We hold the data, the horizon and the period fixed, and vary only the evaluation design. Three choices each reverse or dissolve a headline conclusion. Changing the unit of analysis from the market total to the individual customer moves our production baseline from second of nineteen, beaten by nothing, to twenty-third of twenty-five. Nineteen of its twenty-four challengers beat it there. Changing how much error is pooled decides whether a Diebold-Mariano test finds anything at all. Scoring prediction intervals rather than point forecasts reorders the field almost completely, with a rank correlation of 0.02 on intermittent demand. We then measure what the deployed system gets from this. Its selection rule captures 55% of the distance between doing nothing and choosing with hindsight. The reversal is not a quirk of our data. We ran the released protocol, unchanged, on the public M5 retail panel. The same baseline shape places first at the market total and last per series, beaten by everything, and a replayed selection rule closes 64.7% of the same floor-to-ceiling distance there. Adding five zero-shot foundation models to that roster changes who wins at the total, not the shape. The bands' blind spot travels too: conformal bands under-cover most on the spikiest items. Splitting our own panel into ever smaller groups turns the contrast into a curve: the baseline's rank worsens at every level of disaggregation. We release the evaluation protocol and report an error of our own that inverted a result before we caught it.

In many applications, from government to ecology, integrating data from diverse and noisy sources is critical for downstream inference. However, a unique identifier to link records cleanly from the same entity may not exist. Entity resolution (also referred to as de-duplication or record linkage) merges such databases to identify duplicates, allowing for more complete data for inference. A multitude of methods from statistics and computer science in recent years assume a single, immutable true value for each variable used in linkage. We argue this restrictive assumption is often violated in practice, leading researchers to discard useful data and biasing downstream inference. For example, a respondent's education status may truly change between two different surveys, yet existing methods would treat at least one observation as a distorted version of the truth, even though both are correct. In this paper, we propose a novel entity resolution model that comprehensively accommodates "multiple truths" by introducing a mixture of exponential family distributions that further handles multiple data types. We provide options to fit the model with Markov chain Monte Carlo routines and variational inference for massive datasets. We demonstrate the method's value via simulation and linking a longitudinal survey of Italian household wealth.

The Hosmer-Lemeshow calibration test groups patients by predicted risk for a valid reference and loses power that more groups cannot recover. Later tests weight each record's residual alone, which makes them fragile: on an otherwise correct model, one record in a thousand with a corrupted covariate raises the false-alarm rate of Stukel's score test from 5% to 12%, and ten to 63%. Pooling on equal-size groups bounds a wrong record's influence, dividing its residual by its groupmates' variance; without pooling, the same directions are as fragile as Stukel's test after a refit. EDGE, the test proposed here, projects the grouped residuals onto a small basis of calibration shapes. Its degrees of freedom do not depend on the number of groups, so the partition becomes a setting, fine-grained for power and coarse when records may be wrong, and its reference needs no resampling. In a prespecified simulation EDGE was comparable in power to Stukel's test on clean data and 0.064 more powerful on average than the ten-group Hosmer-Lemeshow test. With one and ten exaggerated covariates in a thousand its false-alarm rates were 0.056 and 0.097, and with ten groups below 0.10 up to about ten exaggerated or five reversed. In external validation, at ten groups, its average power was that of Cox's recalibration test and the GiViTI belt, and its false-alarm rate stayed below 0.10 with ten reversed records in a thousand, where theirs did not. In a 70,000-patient cohort its verdict did not change with the partition.

Exposure to environmental contaminants is a recognized cause of neurotoxicity, contributing to the onset of a broad range of neurological conditions. In realistic settings, such exposure involves com- plex mixtures, and the combined effect of their components may differ from what their individual effects would predict. Characterizing such interactions and testing them against a principled notion of additivity is central to assessing the neurotoxicological risk. We take up these questions for two widespread and independently neurotoxic pollutants, lead (Pb) and benzo[a]pyrene (BaP), through a novel C. elegans assay in which nematodes were subjected to single and joint exposures across a range of doses. Morphological damage is quantified on an ordinal scale at the level of individual dopaminergic neurons. To analyze these data, we model the full distribution of the ordinal response as a convex mixture between an unexposed and a maximally affected profile. The weight of this mixture varies with chemical doses, modeled flexibly via monotone splines and, for the joint effect, in a radial coordinate system. Additivity is assessed via a likelihood ratio test against established null models, and is calibrated via parametric bootstrap. Applied to the C. elegans assay, our analysis reveals a localized, asymmetric synergy between Pb and BaP, concentrated where moderate BaP meets high Pb exposure.

This paper develops a Bayesian framework to estimate the parameters of the Generalized Extreme Value (GEV) distribution for the weight induced by annual maximum accumulations of snow, referred to as the snow load, using a Hamiltonian Monte Carlo (HMC) algorithm as implemented in the \enquote{extremMHMC} R package developed alongside this paper. Key to the approach is the use of strong prior distributions for the shape parameter that are appropriate in the context of snow loads, which helps to ensure robustness in the distribution parameter estimates for annual maximum snow loads despite small sample sizes. This robustness is key to ensuring that structural reliability analyses, which rely on the GEV distribution, produce physically realistic estimates of snow loads. Information on strong prior distributions is derived from existing studies on extreme rainfall and snowfall, with the novel use of hyperbolic tangent functions to transition the prior distribution parameters between low and high snow regimes. This approach enables global applicability of the strong prior approach while maintaining physical realism. Additionally, simulation studies confirm a reduction in Root Mean Square Error (RMSE) of the shape parameter estimate compared to frequentist methods, particularly for small sample sizes. Finally, a real-world application using a data from 9715 stations further demonstrates the feasibility of the proposed Bayesian framework for large scale implementation.

Environmental conditions are major drivers of malaria transmission, but epidemiological analyses are often constrained by temporal misalignment between health outcomes reported at coarse time scales and environmental exposures available at finer resolutions. Conventional approaches aggregate environmental data to match health outcomes, potentially obscuring delayed and nonlinear relationships. We propose a Bayesian spatio-temporal framework that addresses this limitation through a latent daily disease process linked to observed monthly malaria counts by temporal disaggregation. The framework integrates distributed lag nonlinear models for climatic effects, spatio-temporal random effects, and intervention covariates within a unified hierarchical model. The methodology was applied to malaria surveillance data from 161 districts in Mozambique between 2017 and 2024, integrating temperature, precipitation, relative humidity, vegetation, elevation, and malaria interventions. Compared with a conventional monthly model, the proposed framework improved predictive accuracy and uncertainty quantification while exploiting the temporal resolution of environmental data. Estimated relationships showed nonlinear associations between climatic variability and malaria incidence, including an optimal temperature range, increasing risk with positive vegetation anomalies, and nonlinear precipitation effects. By avoiding temporal aggregation of environmental exposures, the framework provides a flexible approach for investigating delayed environmental effects from routine surveillance data and can be extended to other environmentally sensitive diseases with mismatched temporal resolutions.

Timely, paddock-scale estimates of plant-available water (PAW) can support crop management in water-limited grain systems. We present the Sydney Soil Water-Energy Balance (SWEB) model, a scalable, physically consistent modelling framework that integrates satellite and climate data with soil properties to estimate crop evapotranspiration and root-zone soil moisture (RZSM) at daily, 30 m resolution. SWEB was validated nationally against different SM monitoring networks across Australia, which demonstrated robust performance across diverse grain-growing environments (correlation coefficients generally ranged from 0.70 to 0.85 for most regions). In future applications, mid-season PAW nowcasts could be combined with water-use-efficiency approaches to estimate water-limited yield potential and inform responsive management decisions, such as nitrogen fertiliser top-up recommendations.

Coccidioidomycosis is a fungal lung infection endemic to parts of the Pacific Northwest and southwestern United States, Mexico, Central America, and South America. A 2017 study by Thompson et al. reported that many analyzed Coccidioides isolates had elevated minimum inhibitory concentration values for fluconazole but comparatively low minimum inhibitory concentration values for other triazole drugs. We constructed 172 synthetic joint minimum inhibitory concentration profiles from the published drug-specific marginal frequency tables. Values sampled from each marginal distribution were paired across drugs using a Gaussian copula, and potential cross-drug patterns were examined using k-means clustering, hierarchical clustering, and Gaussian mixture models. On the selected reconstruction, a forced four-cluster partition consistently identified an upper caspofungin minimum-inhibitory-concentration tail. Across the evaluated dependence scenarios, the gap statistic favored a single cluster in 96% to 100% of nested reconstructions, providing no evidence that the reconstructed multivariate data supported a broader multicluster structure. A separate simulation study examined recovery of an upper-severity-score category defined from prespecified quantiles of a composite minimum-inhibitory-concentration score. The study included 4,000 replications across 20 prevalence and measurement-noise conditions. The multiclass adjusted Rand index ranged from moderate to high depending on the method and prevalence, whereas binary partition agreement for the upper-severity-score category approached zero at 1% prevalence even though recall remained near 1.0. Clustering conclusions therefore depended on the assumed cross-drug dependence structure and on the realized number of observations in the target category.

We conduct a retrospective analysis of COVID-19 transmission dynamics across U.S. states using a modified population-based Susceptible-Exposed-Infectious-Asymptomatic-Recovered-Deceased (SEIARD) compartmental model. The proposed framework introduces time-varying transmission, reporting, and mortality rates to capture temporal variations in public behavior and policy interventions during the pandemic. In particular, the transmission rate is modeled as a function of population mobility (derived from Google Mobility Reports), with a residual time-decay term capturing the net effect of unobserved factors such as behavioral adaptation and control measures, while reporting is linked to nationwide testing strategies. We employ a Bayesian approach to integrate multiple data sources and quantify uncertainties in model parameters. The model explicitly distinguishes between symptomatic and asymptomatic infectious individuals and links the latent epidemic states to observable quantities, including reported cases and deaths, through a dynamic reporting function. This retrospective modeling framework provides insights into state-level epidemic trajectories and supports data-driven decision-making for optimal allocation of healthcare resources and evaluation of public health interventions during future pandemics. We further apply a clustering analysis to the posterior parameter estimates to identify groups of U.S. states exhibiting similar epidemiological characteristics, revealing substantial regional heterogeneity in transmission intensity, reproduction dynamics, and mortality burden.

This paper develops the statistical foundations of a multivariate general nesting spatio-temporal (MGNST) regression framework for analyzing spatial, temporal, and cross-equation dependence among responses. Four parameter matrices represent spatial lag dependence, spatial error dependence,temporal autoregression, and contemporaneous error covariance. Their off-diagonal elements allow dependence to propagate within and across responses. Matrix restrictions yield eleven specifications encompassing multivariate spatial autoregressive models, multivariate spatial error models, vector autoregressive models with exogenous variables, and independent spatio-temporal regressions as special cases. We establish identifiability conditions using instrumental-variable rank conditions and introduce penalized likelihood estimation and information criteria based on effective degrees of freedom. Monte Carlo experiments examine three data-generating models, three sample sizes, and three levels of spatial dependence. Correct-selection rates under penalized AIC generally increase with sample size and the strength of spatial dependence, while parameter recovery improves as the sample size increases. For socioeconomic data from 198 municipalities in Japan's Kansai region, penalized AIC selects the full MGNST model, whereas penalized BIC selects a response-wise independent spatial error model. Despite selecting models of different complexity, both criteria support temporal persistence and spatial error dependence. The pAIC-selected MGNST model reduces strong spatial autocorrelation in the responses to negligible residual levels, demonstrating its usefulness for identifying and comparing multivariate spatio-temporal dependence structures.

Goal-oriented optimal experimental design (GO-OED) selects experiments according to the expected information gain (EIG) about a quantity of interest (QoI) rather than the full parameter vector. This work develops a variational GO-OED formulation for mixed discrete-continuous QoI laws arising in probabilistic mechanics when thresholding or event-based transformations map a positive-probability set of uncertain inputs to a common value while other inputs produce continuously varying responses. The motivating application is ship roll safety assessment in random waves, where the QoI is the temporal exceedance probability above a prescribed roll-angle threshold. This quantity is zero when no exceedance occurs and varies continuously over positive values otherwise. A purely continuous variational approximation does not dominate a posterior QoI law containing an atom, yielding an infinite Kullback-Leibler divergence and a trivial Barber-Agakov lower bound of $-\infty$. Scoring atom samples using continuous density values instead changes the objective and does not produce a valid lower-bound estimator. We introduce a mixed variational approximation that models the conditional atom probability and continuous component separately, with a normalizing flow used for the latter. An analytical example recovers the correct EIG landscape, while the ship roll application provides stable EIG lower-bound estimates and identifies informative wave conditions for temporal-exceedance-probability inference.

Imputation-based data fusion combines datasets by imputing missing variables in one dataset using information from the other. The validity of imputation-based data fusion relies on accurately preserving both the imputed distributions and the underlying multivariate dependence structure across data sources. However, traditional imputation-based data fusion methods often struggle to preserve the dependence structure, particularly when marginal distributions differ or when complex, nonlinear multivariate dependencies exist. Inadequate preservation of these dependencies can result in distorted joint distributions and biased inference in the fused data. To address this challenge, this study introduces a novel imputation-based data fusion framework that utilises conditional sampling from C-vine and D-vine copulas as the imputation mechanism. This approach flexibly models pairwise and higher-order dependencies while accommodating heterogeneous marginal distributions, thereby generating imputations that more accurately maintain the distribution of the multivariate data. The proposed method is validated through a simulation study using survey data and a real-data application that estimates the negative income distribution in survey data using information from administrative tax records. Results indicate that the proposed method outperforms existing approaches in preserving both the imputed distributions and the dependence structure across datasets.

Regularized Adjusted Plus-Minus (RAPM) uses stint-level lineup indicators to estimate player contributions to scoring margin. When low-minute player columns are removed, their stints remain in the data, but the design matrix no longer represents the complete lineup. Dummy RAPM restores this information using five indicators for the number of excluded players on each lineup side. Across 16 NBA seasons, chronological validation selects a 10-minute-per-appearance threshold and a dummy-to-player penalty ratio of 2.2. On held-out March-April games, Dummy RAPM reduces mean season game-margin RMSE from 12.897 to 12.856 and achieves lower RMSE in 13 of 16 seasons. The average reduction is 0.042 points, or 0.30%. Although the improvement in game-level predictive accuracy is small, it is consistent: RAPM performs better when it records how many excluded players are on each side.

We present a new data-driven learning of a Random Geometric Graph (RGG) of a multivariate dataset, where the graph is drawn in a probabilistic metric space. This graph learning works for generic datasets, irrespective of the type of the observables; their probability distributions; or size of the data. We identify a metric of the space that the graph is drawn in, as a probability distribution of a random variable that we introduce, namely, a variable that represents the disparity between the connectedness of two vertices of the graph, and the correlation between the two random variables that are attached to the respective vertex. It is the closed-form {\it{cdf}} of this disparity variable that we advance as the distance function of the host space of the learnt RGG, such that the edge exists between any two nodes, if this inter-nodal distance falls short of a chosen cutoff probability. Drawing the RGG in this probabilistic space leads to the graph being an Soft RGG, such that any edge - if it exists - exists with an identified probability. We forward a simple Rejection Sampling-based technique for learning the probability of any edge. The expected degree distribution of a vertex of this RGG is identified as local, and dependent on the inter-observable correlation matrix. If said correlation matrix is not known, it can be learnt given the data, using its closed-form posterior probability density function, that we forward. We illustrate our graph learning method by learning multiple RGGs of highly multivariate real datasets.

Despite the wide variety of existing Buy-`Til-You-Die (BTYD) models, nearly all rely upon the convenient assumption of transactions following a Poisson process. As modern customer bases grow larger and more diverse, a major gap in the marketing literature is BTYD models that can account for heterogeneity in timing patterns across millions of customers. This paper addresses that gap, introducing a family of models that assume transactions follow a Weibull renewal process and developing a highly scalable scheme for parameter estimation based on an amortized variational inference procedure. The proposed model fits to a proprietary dataset of 5 million online retail customers in 8 minutes which would take the current state-of-the-art an estimated 3-4 days. We show both theoretically and empirically that this dramatic improvement in computational performance comes with no appreciable change to either model interpretation or predictive performance. Beyond scalability, gradient-based variational inference also makes it easy to extend the model to covariates, which we illustrate on a public dataset of 4 million political donors during the 2020 US General election cycle. More generally, this paper demonstrates how to blend recent advances in approximate Bayesian inference and the tools of modern machine learning to dramatically improve the efficiency and expressivity of probabilistic models for customer base analysis.