Papers

This paper develops a random-effects extension of the ergodic basic affine jump-diffusion (BAJD) model for a population of independent trajectories with individual unobserved jump intensities and common structural parameters. Under continuous-time observation, each intensity is estimated by the empirical jump frequency. We establish the strong consistency and stable mixed-normal limit of this estimator as the observation horizon tends to infinity. The estimated intensities are then used as pseudo-observations to construct a Bernstein-polynomial estimator of their common density. Its bias, variance, mean integrated squared error, and pointwise asymptotic normality are derived under a sequential asymptotic framework. The theory is developed for general positive jump-size distributions with finite activity and specialized to the Gamma$(k,\lambda)$ family, including the exponential and Erlang cases. For this family, we derive a pooled estimator of the common rate parameter and establish its consistency, asymptotic normality, and first-order conditional asymptotic independence from the individual intensity estimators. We also propose a consistent plug-in estimator of the population stationary mean. The methodology is illustrated through simulations and an empirical application to financial realized-volatility data.

Rollcast is a probabilistic forecasting method for univariate time series that combines a compact set of rolling statistical anchors rather than relying on a single global model. Rolling means, medians, extrema, regression endpoints, and quantiles define candidate forecast locations and a representation of the current state. A state-dependent softmax gate learns anchor probabilities by minimizing negative log predictive density, while residual distributions retrieved from similar historical states provide local uncertainty. Recursive simulation propagates the resulting mixture over multiple forecast horizons. The method is evaluated in a Monte Carlo study covering eight data-generating processes, including autoregressive, random-walk, local-trend, threshold, regime-switching, stochastic-volatility, heavy-tailed, and variance-break dynamics. Across 2,000 independent fitted series, Rollcast is compared with the true conditional predictive distribution generated by an oracle simulator. Overall empirical coverage is 86.2% for nominal 90% intervals and 91.5% for nominal 95% intervals. Predictive intervals are on average 13.6% wider than the oracle at the 90% level, while CRPS is 14.4% higher than oracle CRPS. Performance is closest to the oracle under autoregressive, threshold, stochastic-volatility, heavy-tailed, and variance-break dynamics, while local trends and regime switching are more challenging. The results indicate that Rollcast can construct competitive probabilistic forecasts from simple, interpretable local summaries, while also identifying limitations in calibration and recursive uncertainty propagation.

Stochastic self-propagating star formation models describe galactic structure through local triggering and feedback but are commonly formulated using discrete spatial cells and time steps. We extend the spatio-temporal log-Gaussian Cox-Hawkes framework to obtain the Spatio-Temporal Self-Propagating Log-Gaussian Cox-Hawkes Process, a continuous point process model for the locations and times of star-forming events. The model combines spontaneous formation, correlated environmental effects, outwardly propagating excitation, local inhibition, differential rotation, and saturation. An observation layer transforms the conditional event intensity into idealised maps of instantaneous star-forming arm emissivity and recent young stellar surface brightness. For the selected parameter setting, the displayed realisation exhibits transient flocculent spiral-like patterns without any deterministic spiral geometry being imposed. The Monte Carlo experiment finds similar event production with and without rotation, whereas the stationary-front case produces fewer events. The displayed maps further suggest that differential rotation contributes to the winding and spatial arrangement of activity. The proposed framework provides a basis for future statistical inference from spatially and temporally resolved observations of star formation.

In this paper, we study the autocovariance matrix estimation and inference problems under heavy-tailedness, high-dimensionality, general nonlinear temporal dependence, and potentially nonstationarity of time series. We consider two types of tail-robust autocovariance matrix estimation methods: the element-wise Huber's $M$-estimator and a computationally more efficient element-wise truncated estimator. Both estimators are designed to achieve sharp error bounds in matrix max-norm. The nonasymptotic properties of these estimators are proved based on new variants of Bernstein-type inequalities under functional dependence for the potentially nonstationary processes which may be of independent interest. Moreover, we prove a high-dimensional Gaussian approximation result, as a limiting distribution, for our element-wise truncated autocovariance estimator. A Gaussian multiplier bootstrap result is also given to facilitate the practicality. Our theoretical results are nonasymptotic, which gives explicit error bounds in terms of the sample size, dimensionality, moments, and the strength of temporal dependence. Numerical evidence is provided to support our theoretical results. Finally, we illustrate the benefits of the proposed methodology for detecting change points in monthly macroeconomic data.

When a flexible prediction model is fitted on a training subsample drawn from a probability sample, the model-assisted estimator actually reported arises from one realized partition, yet existing theory quantifies uncertainty only for partition-averaged, cross-fitted, or symmetrized versions of it. We represent the training subsample as a second phase of sampling and derive, exactly and for any algorithm, a two-term variance decomposition and the variance family linking the single-partition estimator to its Rao-Blackwellized average, whose design bias it shares. For cell-mean rules under equal-probability designs the estimator is, exactly, a convex combination of a poststratified plug-in at the training-sample size and the exactly unbiased single-split estimator of Sanguiao Sande and Zhang, which locates its bias precisely. For tree-type predictors the second-phase variance is computable in closed form when the cell structure is fixed or $s$-measurable, and in that regime, for cell counts $H=o(n)$, its share of total variance grows with the number of cells – a design-based mechanism distinct from residual shrinkage and consistent with the variance underestimation documented for flexible learners. We propose an analytic and a replication variance estimator, neither altering the point estimate, and evaluate them by simulation under simple random and stratified designs, in the adaptive regime and in an $s$-measurable configuration in which the theory is exact: in the populations studied, budgeting the second phase recovers much of the coverage lost by ignoring it, and the residual shortfall is attributable to the bias of the point estimator rather than to its variance.

Quantile treatment effects (QTEs) measure the effect of a treatment on the distribution of an outcome, and their estimation at extreme quantile levels is of central interest in applications where the target quantiles lie far beyond the range of the data. For heavy-tailed potential outcomes, existing extremal QTE estimators rely on extrapolation combined with a causal extreme value index (EVI) estimator, but the resulting estimator is not invariant under a common location shift of the potential outcome distributions, even though the population QTE is. We address this issue in two steps. First, we adapt the location-invariant Fraga estimator of the EVI to the causal setting using inverse propensity score weighting. Second, we replace the original extrapolation formula with a difference-based scheme, under which the location parameter cancels when quantile differences are taken. The resulting QTE estimator is therefore location invariant. We establish the consistency and asymptotic normality of the proposed extremal QTE estimators, and provide a consistent variance estimator, leading to asymptotically valid inference. A simulation study confirms the location invariance, the stability with respect to the threshold, and the coverage of the proposed methods.

We propose a framework for estimating conditional extreme quantile treatment effects (CEQTEs) in observational studies with heavy-tailed outcomes. Our procedure first estimates intermediate conditional quantiles using inverse-probability-weighted (IPW) quantile regression and then extrapolates them to extreme levels using extreme value theory. Under a linear conditional quantile model, we show that the conditional and marginal distributions of each potential outcome share a common extreme value index (EVI), motivating two complementary Hill-type EVI estimators based on conditional and marginal information, respectively. On the theoretical front, we introduce an IPW tail quantile score process that bridges regression quantile score processes and uniform tail empirical processes while accounting for treatment assignment. We establish its functional weak convergence under mild regularity conditions, without requiring a max-domain-of-attraction condition. This result provides the probabilistic foundation for the asymptotic analysis of the proposed CEQTE estimators. Simulation studies demonstrate favorable finite-sample performance, and an application to NLSY79 data reveals substantial heterogeneity in the effect of college education on extremely high hourly wages across confounder-defined subpopulations.

Small Area Estimation (SAE) and Model-Based Geostatistics (MBG) provide complementary approaches to prevalence mapping, with their relative advantages depending on the inferential goals and characteristics of the available data. We argue that a fuller comparison should consider model interpretability, the role of epidemiologically motivated covariates, inferential objectives beyond area-level prediction, integration of data from different spatial partitions and surveys, and task-specific model validation. In particular, we question whether survey design variables should routinely be included in MBG models when their effects may instead be mediated by measurable environmental and socio-economic risk factors. We further argue that simulation-based validation, tailored to the operational objectives of prevalence mapping, can provide a more informative assessment of model performance than conventional cross-validation alone.

We develop a parametric estimation framework for self-exciting Hawkes processes whose intensity functions admit a parametric form. The estimation procedure is based on ordinary least squares. To apply the least squares estimation, we restrict to a kernel class that can be expressed as a sum of the product of a parameter and a function. We first establish the central limit theorem for the proposed estimator. We then show the consistency of the asymptotic variance estimator. Finally, we introduce a Wald test statistic and derive its asymptotic distribution. We apply the proposed methodology to trade times from high-frequency financial asset data. Our empirical results provide evidence for three heterogeneous trader types. In particular, we identify two trader types as high-frequency traders and one trader type as fundamental trader.

We introduce two new classes of nonparametric change-point tests for sequences of univariate non-negative random variables. The proposed procedures are based on the empirical modified Hankel transform and the Laplace transform, respectively, and provide new transform-based tools for detecting distributional changes. We derive the asymptotic null distributions of the corresponding test statistics and use a permutation bootstrap procedure to obtain $p$-values, since the limiting null distributions are not distribution-free. Through a finite-sample simulation study, we show that the proposed tests are well calibrated, generally more powerful than empirical-characteristic-function-based competitors, and capable of accurately estimating the change-point location, particularly in asymmetric settings. The usefulness of the proposed methodology is further demonstrated through applications to Argentina rainfall data, as well as U.S. GNP and S&P 500 absolute log-returns. The real-data examples illustrate that the proposed Hankel- and Laplace-transform-based tests are effective tools for detecting meaningful distributional changes in non-negative data.

We present a framework for the analysis and unmixing of random density mixtures in the Bayes Hilbert space. General identifiability results for mixtures in Hilbert spaces are established and applied to the Bayes Hilbert space setting. Building on these results, we propose a penalised maximum likelihood approach for the unmixing of Bayes Hilbert mixtures aimed at recovering the statistically space-efficient representation, together with a computationally efficient coordinate-wise maximisation algorithm for its implementation. The methodology is illustrated through a hyperspectral data application, where observations can be naturally embedded in the Bayes Hilbert space and analyzed in terms of distributional shape rather than amplitude. A complementary simulation study demonstrates the interpretability and practical performance of the proposed approach.

We give an exact randomization-based confidence set for the average treatment effect (ATE) in matched-pair studies with a binary outcome, requiring neither monotonicity nor any distributional assumption beyond the within-pair coin flip. At its core is an analytic solution to the worst-case allocation of attributable effects: two binomial-symmetry lemmas identify the pattern hardest to reject as a single boundary corner, so testing null hypotheses needs no integer program and no numerical search. Inverting the test via binary search yields a prediction set for the attributable effect in O(log S) Binomial tail calculations; the Bonferroni proposition of Rigdon and Hudgens (2015) produces the ATE confidence set at the same computational cost. The same corner extends without further machinery to a sensitivity analysis for matched observational studies under Rosenbaum's $\Gamma$-model. A simple formula for the design sensitivity illuminates when an observational study can hope to provide evidence for an effect.

The existing multimarker survival tests focus on time to event outcomes. However, recurrent events are common in real world clinical and biomedical studies, especially in the research of chronic and recurrent diseases. In this paper, we develop a suite of set based genetic association tests for panel count outcomes under a unified weighted V statistic framework. These tests can effectively account for genetic effect heterogeneity in panel count data. Additionally, we develop small sample corrections to the tests to enhance the accuracy of the tests under small samples. Simulation studies show that the proposed tests perform well in terms of size and power across various scenarios and have a bigger power than the existing set based tests for interval censored outcomes. A dental caries GWAS data set is analyzed to illustrate the utility of the tests.

A probability coordinate chart is a continuous strictly increasing bijection that transports observations to the open unit interval, where averaging is always well defined. For barycentric inference, the first coordinate moment is returned through the inverse chart. Higher initial coordinate moments similarly generate initial Kolmogorov moments by pullback, while centred coordinate moments remain on the probability scale. We develop inference for these quantities under fixed, intrinsic, and estimated charts. For a fixed benchmark chart, the coordinate mean is the primary inferential object. Confidence intervals are constructed on the probability scale and then transported through the inverse chart, preserving the geometry and allowing asymmetry on the observation scale. Hoeffding's inequality also provides finite-sample distribution-free intervals. The intrinsic case requires different treatment. The empirical cumulative distribution function is not an admissible chart, and its natural generalised plug-in construction collapses exactly to a middle order statistic. Feasible intrinsic inference therefore reduces to median inference. This self-induced functional is distinguished from the frozen oracle construction, which has a different asymptotic variance. For estimated location-scale charts, the asymptotic expansion contains a calibration correction for the first-order effect of learning the chart from the same sample. This yields an influence-function variance estimator and a bootstrap procedure that recalibrates the chart in every resample. Joint covariance theory is developed for vectors of initial and centred coordinate moments, with initial Kolmogorov moments obtained by inverse-chart pullback. The framework separates variation on the probability scale, amplification by the inverse chart, and uncertainty from chart calibration.

We used nested models as a framework for characterising the consequences of a selectively measured predictor in clinical prediction models. In this framework, a model containing predictors available in the target population is referred to as the restricted model, while the extended model additionally includes a selectively measured predictor. We show that non-invariance in the restricted model between selected patients and the target population decomposes into components due to omission of the additional predictor, potential residual non-invariance, and their interaction. This framework is extended to a predictor measured via multiple routes of selection resulting in collider structures. We further show that imputation based on the conditional distribution of the additional predictor in the selected population transfers the restricted-model non-invariance to the imputed extended model as imputation bias. In a forthcoming version, this framework will be illustrated using real world and simulated data. The framework provides a formal basis for understanding how selective predictor measurement affects generalisability of clinical prediction models using routinely collected health data.

Dependence can change between variable blocks without changing correlation; serial dependence complicates calibration. We propose DOMIC, a rank-based detector using density-operator mutual information (DOMI). Unit-norm random features make sample-level partial traces exact, so prefix sums compute segment statistics at per-split cost independent of segment length; a Gram form uses the exact kernel. Permutation tests are finite-sample exact under pair or block exchangeability for windows and schemes fixed in advance; validity after data-driven scheme selection is unproved. For i.i.d. pairs with feature dimension, frequency draw and bandwidth fixed, we prove a weighted chi-square limit at independence for the rank-based random-feature statistic. An additive entropy cost with Holevo gains drives segmentation. DOMI attains higher localised power than rank, copula, kernel and distance baselines in most non-Gaussian and correlation-free settings; the Gram form leads or ties in all 16 non-Gaussian settings where any statistic exceeds 0.10. On a three-break benchmark, Holevo partitioning returns exactly three breaks in 47.5% of replicates, against at most 19% for calibrated baselines. Block-calibrated primary scans reject in all three pre-specified pairs of a US financial panel, with strongest stage-two evidence for stocks and Treasury yields. One of 27 weather candidates passes the multiplicity-corrected re-test; candidate selection is unadjusted.

In the presence of interference, where the treatment assigned to one unit can affect the outcomes of others, many causal estimands depend on the treatment-assignment policy under which the experiment is conducted. This policy dependence creates a fundamental challenge for off-policy estimation, where the goal is to estimate causal quantities under a hypothetical intervention policy different from the one used to collect data. We study this problem of off-policy estimation of causal effects for heterogeneous Bernoulli policies. By representing exposure-weighted potential outcomes in the biased Fourier basis of the experimental design, we construct, for any prespecified Fourier subspace encoding the assumed interference structure, the unique minimum-$L^2$ weight that transports every function in that subspace. Global and local inverse-probability weights, linear-interference weights, and no-interference weights are special cases. The weight variance is a structured chi-square distance between the experiment and target policies. When the assumed interference structure is misspecified, the introduced bias couples the omitted outcome spectrum with the corresponding policy-shift coefficients, yielding a sharp robustness bound and a bias-variance trade-off. A Fourier-neighborhood-overlap condition gives consistency under structured interference, and we state a Doob-martingale central limit theorem for off-policy estimators. As the variance is not identified, we derive identifiable bounds and associated conservative estimators of the variance. Simulations illustrate these theoretical results for the design and analysis of experiments under network interference and design mismatch.

We establish the consistency and asymptotic normality of a two-step estimator of conditional expectiles in the context of conditional scale models. We first estimate the conditional variance parameters by quasi-maximum likelihood and then compute the unconditional expectile of the innovations using the empirical distribution of the standardized residuals. We show how replacing true innovations with standardized residuals affects the asymptotic variances of both conditional and unconditional expectile estimators. Finally, our empirical analysis reveals that conditional expectiles assess tail risk in cryptomarkets in a more robust manner than traditional quantile-based risk measures, such as value at risk and expected shortfall.

A central challenge in distributional regression is to allow the shape of the conditional distribution of the response variable to vary flexibly with covariates while retaining directly interpretable effects on its mean and standard deviation. We extend the penalized transformation model (PTM) family into a conditional-shape PTM, which assigns separate structured additive predictors to the conditional mean, standard deviation, and standardized distributional shape beyond location and scale. A covariate-dependent monotone transformation maps the standardized response to a fixed reference distribution, while affine standardization enforces mean zero and variance one for the induced standardized distribution. Thus, the first two predictors remain exactly the conditional mean and standard deviation. The shape predictor accommodates selected linear, nonlinear, group-specific, spatial, and interaction effects; regularization shrinks unsupported departures toward a reference-family location-scale model. We fit the PTM using mini-batched stochastic variational inference with model-aligned Gaussian blocks and staged optimization. In simulations, the PTM recovers smooth mean and standard-deviation effects and a covariate-dependent transition from skewness to bimodality while suppressing unnecessary shape effects. Under a deliberately misspecified design, it remains competitive with a structured additive Dirichlet-process mixture in test-set density and distribution-function accuracy, although all models show undercoverage and both flexible methods miss fine features. Applications to $13{,}425$ Norwegian water-conductivity observations and $1{,}182{,}514$ German daily-temperature observations demonstrate selective group-specific, seasonal, and spatial shape variation. Predictive performance criteria favor the conditional-shape PTM over a fixed-shape PTM and a Gaussian location-scale model.

Clinical trials often show treatment curves that diverge early and converge later, or vice versa patterns that are poorly captured by the proportional-hazards assumption. We develop a joint inferential framework for two nonparametric functionals of censored survival data: the Kaplan–Meier-based Mann–Whitney effect and a novel temporal contrast separating early and late differences. The approach provides interpretable, probability-scale effect measures and enables joint inference for global and temporal contrasts under right censoring. In simulation studies, the method outperforms the log-rank test under non-proportional hazards while maintaining nominal type-I error. A real-world application illustrates how the temporal contrast reveals clinically meaningful early treatment advantages that remain hidden in standard analyses

In this lecture note, we present a concise and self-contained derivation of the discrete-time Kalman filter equations that requires only a basic understanding of least squares estimation. The treatment is designed to minimize mathematical overhead while preserving both rigor and generality.

Computationally expensive black-box models often involve a large number of input factors with complex interactions and varying importance. Experimental design techniques can be used for quickly identifying the important factors, which can make the optimization of a complex computer model or the training of an expensive machine learning model more efficient. Existing screening designs for black-box models focus mainly on continuous factors, with the maximum one-factor-at-a-time (MOFAT) design being a recent example. In this work, we extend the design to incorporate multiple types of factors, including nominal, ordinal, and discrete-numeric. We first identify the properties leading to optimal screening, where qualitative and quantitative factors should be treated differently. Based on these properties, we propose practical algorithms to efficiently construct MOFAT designs for all types of factors. The usefulness of the design is demonstrated by both numerical experiments and an application to hyperparameter tuning in machine learning models.

Factorial designs allow simultaneous estimation of multiple main effects and interactions, but covariate imbalance can substantially reduce estimation precision. Existing rerandomization methods improve covariate balance yet do not fully exploit heterogeneous priorities across factorial effects or effect-specific covariate importance. To address these limitations, this paper proposes a data-adaptive rerandomization framework for $2^K$ factorial designs. We first develop an oracle criterion that jointly incorporates researchers' priorities over factorial effects and effect-specific covariate importance, enabling precision gains with guaranteed lower bounds. To make the oracle criterion implementable, we develop a data-adaptive procedure that learns effect-specific covariate importance from a random subset of units and applies an estimated oracle criterion to the remaining units. Unlike existing two-stage rerandomization methods for treatment-control experiments, our procedure accommodates multiple factorial effects and requires no auxiliary dataset. Under a finite-population framework, we establish design-based asymptotic theory and show that the proposed procedure preserves the oracle design's precision-prioritization property and, under suitable conditions, achieves the same asymptotic precision as the oracle design. Numerical studies demonstrate substantial efficiency gains over existing rerandomization methods.

Experiments on networks observed over time face network spillovers, temporal carryover, and dependence deliberately introduced by the design. We propose COSTA—Covariance-Optimized Spatiotemporal Treatment Allocation—a joint Bernoulli design for unit–time assignments. Under common treatment marginals and a nonnegative linear network–temporal exposure model, Horvitz–Thompson bias for the sustained all-treated versus all-control contrast is exactly the negative expected weight of an assignment cut. A covariance-level variance envelope yields an MSE bound that can be optimized directly over assignment covariance. To scale this design, we introduce a thresholded-Gaussian Kronecker parameterization that mirrors the network and temporal exposure operators while preserving valid Bernoulli marginals. We next develop inference theory for the joint effects of designed treatment dependence and interference-induced outcome dependence. Canonical correlations between latent blocks generating separated HT contributions supply the coefficients required by graph-$\psi$ central limit and network-HAC theory; a spectral-floor and far-row-mass condition gives a primitive sufficient check. The framework covers sparse, block, Kronecker, locally factored, and other structured covariance sequences satisfying these conditions. Semi-synthetic RetailRocket and MovieLens experiments show substantial default-setting RMSE reductions and well-calibrated model-assisted design-centered intervals across linear, nonlinear, and demand-substitution outcome surfaces.

Time-varying treatment effects can be estimated with linear regression paired with jackknife confidence intervals. In experiments at a large online video platform, the test based on this estimator is not much less powerful than modern tests for monotonic trends. If you suspect trends in your online experiments, consider using this estimator.

To deal with superfluous clusters and to reduce the number of clusters as often desired in practice, a modal cluster merging procedure is proposed. Based on some new properties established in this paper for Morse functions, the procedure merges clusters in a sequential manner without causing unnecessary density distortion. Each cluster is evaluated for its significance relative to the other clusters, using the Kullback-Leibler divergence or its log-likelihood approximation, by truncating the density for the cluster at an appropriate level. The least significant cluster is then merged into one of its adjacent clusters, using the novel concept of cluster adjacency defined in this paper. The resulting hierarchical clustering tree is useful for determining the number of clusters, as may be preferred by a specific user or in a general, meaningful manner. Numerical studies show that the new procedure deals well with difficult clustering problems and often produces intuitively appealing and numerically more accurate clustering results, as compared with several other popular clustering methods in the literature.

Contemporary oncology drug development increasingly requires efficient dose-optimization strategies that evaluate monotherapy (Mono) and combination therapy (Combo) while balancing activity, efficacy, and tolerability. We propose a staggered seamless phase I/II design for settings in which a novel agent is evaluated alone and in combination with an established therapy. In phase I, Mono dose finding begins first, and Combo subtrials can be opened adaptively once a prespecified combination-initiation signal based on early clinical or biological information is observed. Dose assignment uses a model-assisted rule based on toxicity and early activity, with backfilling at tolerable and potentially promising regimens. At the end of phase I, two candidate regimens are selected from the evaluated Mono and Combo regimens using an efficacy-toxicity utility based on accumulated toxicity and treatment-response data. Phase II seamlessly carries forward patients treated at the selected regimens, enrolls additional patients as needed, and applies Bayesian futility and efficacy stopping boundaries to identify a final recommended optimal biological dose (OBD). Simulation studies showed that the proposed design shortened phase I trial duration relative to the comparator designs while maintaining competitive OBD-selection performance and acceptable safety. The seamless phase II component further reduced the need for additional enrollment and supported efficient final OBD selection.

Conditional quantiles are central to asymmetric decision losses, tail-risk assessment, and interval forecasts, but Bayesian quantile regression for nonlinear time series can be difficult when the temporal feature vector is high-dimensional. We develop the quantile deep echo state network (Q-DESN), a Bayesian quantile regression model conditional on fixed features generated by a deep echo state network. Conditional on a specified reservoir construction and feature map, posterior uncertainty is assigned to the regression coefficients, likelihood, and shrinkage parameters. Single-level fits use asymmetric Laplace or quantile-fixed generalized asymmetric Laplace working likelihoods with ridge or regularized-horseshoe priors. Posterior computation uses Markov chain Monte Carlo when computationally practical and a model-specific variational Bayes approximation for analyses requiring repeated fitting. For quantile grids, we compare independent level-wise regressions followed by monotone rearrangement with a joint quantile-vector regression that shrinks adjacent quantile-specific coefficient differences. Synthetic studies, a single-origin retrospective Global Flood Awareness System (GloFAS) streamflow case, and a retrospective PriceFM comparison identify settings where this fixed-feature Bayesian regression improves finite-grid quantile scores.

Here we present a method for measuring, analyzing and reducing theoretical uncertainty. We call it metatheoretical multiverse analysis (MMA). The method is important because the reliability and replicability of scientific observation and testing of a phenomenon are a function of theoretical uncertainty thus reducing it is a method for improving theory. The MMA method requires encoding theory as propositional logic between variables (nodes) and their relationships with one another (edges). Propositional logic models are treated as causal path models in our method, which enable mathematical properties of testing and identification underlying their causal propositions. By encoding theories-as-data including their unknown components, it is possible to generate a multiverse of alternatively plausible theoretical models which can then be meta-analyzed using techniques borrowed from (statistical) multiverse analysis. We introduce metrics for studying theoretical multiverses that enable measuring and identifying the causes of metatheoretical uncertainty. This guides researchers to where it is best to invest theoretical development to reduce uncertainty. We use three simulations and a dedicated software package to demonstrate this method.

Large spatial data sets now record many correlated variables at many thousands of locations, often on domains where Euclidean distance misrepresents proximity. The central difficulty is modelling the cross-variable dependence jointly while retaining variable-level interpretation. We introduce the bigraphical Mat'ern-Whittle process, a multivariate Gaussian process that resolves this with two graphs. A spatial graph generates the Mat'ern structure of each variable through a fractional power of a graph Laplacian, so the process is valid on any topology, with per-variable range, smoothness and amplitude. A directed acyclic variable graph encodes the scientific structure: we prove that each absent edge yields an exact conditional independence between the corresponding fields. We further prove that the operator determinant does not involve the cross-dependence coefficients, which keeps matrix-free likelihood evaluation and Bayesian learning of the variable graph tractable at scale. Estimation requires only sparse matrix-vector products and scales to tens of millions of space-variable pairs. In simulations the method recovered parameters and graphs accurately, remained robust under misspecification, and halved held-out prediction error on a non-convex domain. In a spatial transcriptomics section with 19,809 cells and 1,122 genes, fitted in 75 minutes on a laptop, borrowing across the learned gene graph reduced held-out prediction error by 50 to 91 percent. Theoretical challenges, such as the achievable efficiency of estimating the variance of the nugget, are also explored.