Papers

The polynomial maximization method estimates signal parameters in non-Gaussian noise from a finite description of the noise by its moments and cumulants, without a density. Known accuracy results assume that these noise characteristics are known, or use a design in which estimating them jointly costs no accuracy. In practice they are estimated from the residuals of a preliminary fit, and the effect of this substitution on accuracy is unknown. We treat the feasible scheme as a two-step estimator and, for a general smooth signal model in a fixed design, find the increment of its asymptotic covariance in closed form: $\Sigma_{\rm fe}-\Sigma_{\rm or}=c_2(1-g_{\rm or}) Q^{-1}\bar{d}\bar{d}^{\top}Q^{-1}$, a positive semidefinite matrix of rank at most one. Here $c_2$ is the noise variance, $g_{\rm or}\le 1$ the oracle variance reduction factor, and $Q$ and $\bar{d}$ the second-moment matrix and mean of the signal gradient over the design. Estimating the fractional moments that form the weights costs no first-order accuracy, because that block is orthogonal; the whole increment comes from estimating the sign-preserving centering means. The increment lies along the single direction $Q^{-1}\bar{d}$, so linear combinations of the parameters orthogonal to it keep oracle accuracy. In the trace metric the variance reduction factor becomes $g_{\rm fe}=g_{\rm or}+(1-g_{\rm or})R_D$, where $R_D=\|Q^{-1}\bar{d}\|^2/{\rm tr} Q^{-1}$ depends on the design alone. We give four sufficient conditions for first-order oracle equivalence; noise symmetry is not among them. Monte Carlo experiments with a harmonic signal confirm the direction and size of the increment. A covariance estimate that ignores the increment undercovers where non-Gaussianity gives the largest gain, symmetric noise included.

Clinical trials with survival endpoints lose information when participants are censored before death is observed. We develop target-preserving estimators that use posttreatment disease history, such as recurrence or progression, to recover information lost to censoring for marginal survival and restricted mean survival effects. The difficulty is that the intermediate event is downstream of treatment: naive adjustment can change the causal estimand, and the useful information enters only through the observed coarsening. We derive observed-data influence functions with and without recurrence history and obtain an exact gain identity. The identity shows that efficiency improvement is driven by the censoring hazard, the split of the alive risk set into recurrence states, and the residual-survival separation between those states. In the no-covariate illness-death model, the Aalen-Johansen estimator realizes the recurrence-augmented efficient score after standardization to the marginal target. With covariates, correctly specified Cox-Breslow transition hazards provide a root-N plug-in benchmark, while a hazard-induced one-step estimator gives rate robustness and, under primitive transition and censoring learner rates, canonical inference with a second-order product remainder. A semi-synthetic metastatic breast cancer study calibrated from digitized progression-free survival and overall survival curves illustrates the gain identity. The framework applies broadly to censored time-to-event studies with informative intermediate histories.

Dose optimisation of two-agent combinations in early-phase oncology trials must jointly balance toxicity, efficacy, and pharmacologically meaningful drug exposure, yet existing model-assisted combination designs use only binary toxicity–efficacy outcomes and leave pharmacokinetic (PK) data unused for dosing decisions. We propose PKComb-BOIN12, a utility-based Bayesian optimal interval design embedding a continuous PK exposure metric into two-agent dose optimisation via a PK-guided escalation rule and a PK-sufficiency constraint on the final admissible set. Two contributions distinguish this design from a direct combination extension of existing single-agent PK-informed designs. First, PK-based elimination uses a coordinatewise order-ideal (down-set) removal rule, so that eliminating a combination for insufficient exposure is licensed exactly by the assumed monotonicity of exposure in each agent's dose, rather than an unjustified scalar summary of combined dose intensity. Second, the pre-specified target exposure is treated as transportable-but-verifiable: a hierarchical calibration mechanism updates the assumed monotherapy-derived target using accruing combination-arm PK data, with a formal bound on the resulting selection bias under transportability violations. A time-to-event extension, TITE-PKComb-BOIN12, accommodates late-onset outcomes. A factorial simulation study spanning realistic scenarios, including one calibrated to a published combination trial, decomposes the source of PKComb-BOIN12's gains over its non-PK predecessor and evaluates robustness to non-monotone exposure, PK-variance misspecification, and target-transportability violations. An R Shiny application supports design evaluation and real-time conduct.

Spatially balanced sampling designs, such as the local pivotal methods (LPM), improve representation of auxiliary-variable space, but design-based variance estimation is challenging because second-order inclusion probabilities are generally intractable and may be zero or near zero. We propose two approaches for estimating the variance of the Horvitz-Thompson mean and the finite-population variance. The first estimates joint inclusion probabilities by Monte Carlo simulation and uses pairwise independence tests to modify the resulting variance estimators. The second, a model-assisted bootstrap (MABS), uses an interpolating Gaussian process to construct a synthetic finite population and estimates variance by resampling under the original LPM design. In simulations, the modified empirical approach performs well when responses are strongly associated with auxiliary variables, but its performance for the Horvitz-Thompson mean deteriorates as this association weakens. MABS exhibits more stable performance across noise levels and sampling settings, including under LPM2 and non-Gaussian noise. For finite-population variance estimation, however, no method shows a uniform advantage. These results illustrate the complementary strengths of the proposed approaches for inference under spatially balanced sampling.

Fréchet regression provides a framework for predicting random object-valued responses by minimizing the conditional expected squared distances. We propose semi-supervised kernel ridge Fréchet regression (SS-KRFR) to exploit unlabeled covariates when labeled samples are scarce. The method learns kernel spectral features from both labeled and unlabeled covariates, then utilizes labeled responses to estimate the conditional Fréchet function by kernel ridge regression. We derive prediction error bounds that distinguish the roles of labeled and unlabeled samples, and establish polynomial convergence rates under suitable conditions. Under additional kernel assumptions, a sharper analysis relaxes the sample size requirements for learning the features. We also develop hybrids that combine features, predictions, or both from models fitted with and without unlabeled covariates, establishing uniform stability over mixing weights in Hadamard spaces at common tuning and a separate rate guarantee for prediction mixing. Simulations and a real-data application show semi-supervised gains and improvements over competing Fréchet regression methods in several settings.

Multi-step forecasting systems, including modern artificial intelligence (AI) models, often issue overlapping predictions for the same future time from different forecast origins. These predictions can disagree even when produced by the same model, a phenomenon we call cross-origin incoherence. We propose a Bayesian reconciliation framework that treats overlapping forecasts as noisy and potentially biased measurements of a common latent trajectory. The method operates as a model agnostic post-processing layer, requiring no access to the forecasting models architecture, parameters, or training data. A Bayesian state-space model jointly reconciles forecasts targeting the same future time while accounting for horizon-specific bias and uncertainty. Controlled simulations demonstrate substantial reductions in cross-origin incoherence and accurate latent trajectory recovery under the assumed data generating process. An empirical application to daily S and P 500 volatility forecasts reduces 90% prediction interval coverage from 94.8% to 88.3%, closer to nominal coverage, but increases point forecast error. An additional LSTM experiment demonstrates applicability to modern neural forecasting systems.

We present a unified, transformation-free linear regression framework tailored for ordinal compositional data, such as distributions across educational levels or aggregate Likert-scale survey responses. Traditional log-ratio approaches often obscure interpretability, ignore the inherent ordering of categories, and struggle with boundary zero values. To overcome these limitations, we propose a deterministic model constrained to the space of column-stochastic transformation matrices, naturally preserving the fundamental geometry of the simplex. By adopting the weighted 1-Wasserstein distance as the loss function, our method directly embeds the ordinal nature of the data into the estimation process. We provide a computationally efficient and globally optimal solution explicitly formulated as a Linear Programming (LP) problem. The framework systematically addresses four distinct regression scenarios, introducing a novel ordinal tensor product to rigorously handle interactions across multiple compositional predictors. Furthermore, we equip the model with a suitable regularization strategy and novel diagnostic tools, including a Wasserstein-based coefficient of determination ($R^2_W$) and an Order Preservation Index (OPI). Extensive simulation studies and two empirical applications demonstrate the robust predictive performance and interpretability of the proposed methodology.

Patient world models used for clinical trial simulation can agree on transition kernels and arm-specific risks, yet disagree on the fraction of patients harmed by switching treatment—the counterfactual quantity that matters for intervention-aware reasoning. We audit this reliability gap in a two-stage shared-response SCM: a categorical intermediate health state is followed by common terminal care. Under independent stages, the sharp harm interval has closed-form endpoints for at most three intermediate states, with an exactness boundary at four states. Declared dependence and response-mismatch budgets yield calibrated outer bounds when stage independence or complete mediation is relaxed; in a symmetric three-state model the entire sensitivity frontier is sharp, $[0,\min\{1/2,1/3+(\rho+\delta)/2\}]$, and shows exactly how budgets erase the gain over endpoint-only bounds. Two eight-variable response LPs propagate interventional uncertainty for finite-sample audits. Exact witnesses verify attainability. On public clinical simulators (EpiCare; sepsis), native configurations show little resolved stage dependence and no additional joint-compatibility gain over pairwise transport—honest negative results for reliability claims. All experiments are locally reproducible; guarantees remain conditional on the stated causal model. The results provide a concrete protocol for deciding when a patient world model is safe to trust for counterfactual harm.

Two-sample Mendelian randomization with summary statistics has become a popular MR design in Epidemiology. Its ratio estimate has an average causal effect interpretation only under strong untestable assumptions including linearity. We derive what the per-variant ratio, and the pooled and robust estimators built from it, estimate when a continuous exposure has a nonlinear effect, under two conditions on how genetic variation and unmeasured factors shape exposure: no unmeasured common effect modifier, and an exposure control function. Under either condition, each variant's ratio is a weighted average of the slope of the dose-response curve, with weights that depend on how that variant moves the distribution of exposure and that are nonnegative when the exposure distributions at different genotypes do not cross. Variants that shift everyone's exposure equally target nearly the same quantity; variants acting at particular doses or on the spread of exposure target different ones, some outside the range of any individual effect. Pooled estimates therefore have no fixed interpretation, heterogeneity tests and MR-Egger misread nonlinearity as pleiotropy, and summary statistics cannot reveal the problem, even when the dose-response curve is identifiable from individual-level data. For continuous exposures, two-sample estimates rest on a linearity assumption that should be stated and defended; individual-level exposure data resolve the problem and should be prioritized.

Finite discrete count data are ubiquitous in fields such as biomedicine and social sciences and these data frequently exhibit over-dispersion, rendering the binomial distribution inadequate for modeling. Although the beta-binomial distribution mitigates this by assigning a subjective beta prior to the success probability, its reliance on a pre-specified parametric shape may result in biased representations of prior information, failing to capture the complex and heterogeneous characteristics of real-world data. To address these limitations, this paper proposes a novel Bernstein-binomial model. By leveraging the uniform approximation properties of Bernstein distribution, we construct a highly flexible and general prior framework for the success probability in binomial distribition, which dynamically adapts to complex data structures, avoiding the subjectivity and restrictions associated with conventional unimodal priors. A comprehensive theoretical framework for both frequentist and Bayesian inferences is established. Furthermore, the model is extended to a regression setting to account for covariate effects. Extensive simulation studies and real-world data applications confirm that the proposed Bernstein-binomial model significantly outperforms existing models. It offers superior fitting accuracy, robustness, and out-of-sample predictive performance when modeling over-dispersed counts with complex latent prior structures.

Noisy matrix completion is a fundamental problem in statistical learning and has attracted a substantial amount of interest over the last two decades. In a variety of applications, the missingness pattern is highly informative, yet it has received relatively less attention. Most existing methods either overlook this source of information or oversimplify its generating process, with few attempting to model the association between the data matrix and its informative missingness. In this study, we propose a general modeling framework that jointly models the data matrix and its missingness pattern using a pair of low-rank models coupled by a flexible shared linear structure. This joint low-rank approach incorporates the informative missingness to improve matrix completion and can flexibly capture various associations between the two modes of data. We develop an efficient joint estimation procedure for this framework based on projected gradient descent, and establish local convergence guarantees that unveil its computational and statistical properties. We further demonstrate, through extensive simulation studies and real-world data analysis, that our proposed approach outperforms competing methods, achieving substantial improvements in prediction accuracy on missing entries.

We present a Bayesian framework for assessing the calibration of nested uncertainty sets from simulation studies. Calibration is typically evaluated by repeatedly simulating data from a known model and comparing empirical coverage probabilities with their nominal values, an approach that often requires hundreds or thousands of independent model evaluations to obtain reliable estimates. Such computational demands can be prohibitive when the forward models are expensive to evaluate, for example in data assimilation using a large Earth-system model, or may not be applicable, for exmaple when the output of the simulation is a single or multiple uncertainty sets rather than a full distribution. We present a method that leverages the nested structure of uncertainty sets across nominal coverage levels to assess the calibration of these models. The method, which is based on a Dirichlet-Multinomial model, allows for testing both exact and region-based notions of model calibration in an efficient manner using Bayes factors. Unlike methods that require a full predictive distribution, this method works for procedures that directly produce prediction or confidence intervals. We demonstrate the method on a data assimilation problem containing model form error and show that it is capable of detecting miscalibration even when there are only a limited number of simulation studies available.

Sampling from a probability distribution with access only to its log-density is a foundational problem in statistics, machine learning, and the physical sciences. This problem is studied in the data-free score-based diffusion literature, which samples the target density via estimating scores of densities produced by convolution with Gaussian random variables and implementing a reverse SDE. In this framework, score estimators are constructed via Markov chain Monte Carlo, importance sampling, or rejection sampling. The estimators often operate far beyond the regime where mixing time and error control can be established and have large bias and variance when estimating the underlying score. In our method we propose post-processing: splitting a fixed computation budget across multiple evaluations, querying the estimator at different points, and fitting a single pooled estimate via linear combination to minimize bias and variance at a target location. The final value is an affine map of a vector of multiple score estimates. We show this reduces score estimation error when compared to the computation-matched base estimator queried at a single point. We demonstrate improved sample quality over other score estimation methods matched for computation cost on difficult target densities such as log-concave densities with large condition number, multimodal densities, and high-dimensional densities.

The log predictive density score (LPDS) and related estimators of the expected log predictive density (ELPD) have become the de facto standard for Bayesian model comparison. The LPDS is a natural out-of-sample measure of predictive density performance that accounts for parameter uncertainty, is less dependent on the prior than the marginal likelihood, and is easily computed sequentially using standard Monte Carlo from posterior draws. Our article demonstrates, however, that the standard Monte Carlo estimate of the LPDS can be extremely variable and severely biased when the test data has outliers, leading to misleading conclusions and suboptimal model choices. Estimating the LPDS as the difference between the log marginal likelihood of the full dataset and the training data is shown to be a dramatically more efficient approach. This estimator is also shown to be a much more effective fallback for observations that have been flagged as problematic by the widely used PSIS-LOO model comparison framework in Stan.

Health technology assessment often aims to balance evidence across multiple treatments and outcomes and assess treatment benefit-harm profiles. Multivariate network meta-analysis (mvNMA) is well suited to this task, as it synthesizes evidence across treatments and outcomes. However, the potential of mvNMA has been overlooked, leading decision-makers to rely on outcome-specific NMAs that ignore outcome correlations. Based on independent NMA models, previous work proposed methods such as multivariate P-scores or spie charts to assess treatment benefit-harm profiles by yielding across-outcomes treatment hierarchies. In this article, we propose an mvNMA-based framework for assessing treatment benefit-harm profiles. Our framework uses a hybrid-correlation mvNMA model and incorporates DuMouchel priors for treatment effects, enabling estimation of all treatment effects even when outcomes are missing for some treatments. We then extend common NMA ranking metrics, including SUCRA, the probability of being best, and median ranks, to the mvNMA setting. Based on these metrics, we adapt the VIKOR algorithm for mvNMA. This deterministic multicriteria decision method quantifies overall and worst-case treatment performance and, by balancing these aspects, produces an across-outcomes ranking while identifying treatments offering the optimal compromise between these two aspects. To facilitate use of our framework, we developed the R package mvnma. We applied our framework to two clinical datasets: one comparing eleven physical therapy treatments across two outcomes and another comparing nine antidepressants across five outcomes. Our framework provides an alternative strategy for assessing treatment benefit-harm profiles across multiple outcomes and encourages multivariate approaches in applied clinical and policy settings.

Penalized function-on-function regression, as implemented in the R package refund, treats all grid values of all curves as independent. Residual curves are smooth, so values within a curve are dependent. We study what this does to point estimates and pointwise confidence intervals, in a simulation with Gaussian, count and binary responses and in a plasmode study built on nine real datasets. Under within-curve dependence, model-based 95% intervals cover 0.63–0.82. A curve-clustered sandwich covariance with a leverage adjustment for the penalized hat matrix (CL2), computed after the usual REML fit, brings coverage close to nominal in the simulation; in the plasmode study it is a few points short on eight datasets and 0.84–0.91 for the coefficient surface on the ninth. Satterthwaite degrees of freedom recover most of the remaining shortfall except on that dataset, also with only 20 curves, at the cost of wider intervals. REML, however, undersmooths the coefficient surface: its error is two to four times that of a fit whose smoothing parameters are chosen by curve-blocked neighbourhood cross-validation (NCV), and its CL2 intervals for the surface are wide. Intervals centred at the smoother NCV estimate are narrow but undercover, the familiar tension between estimation and pointwise inference in nonparametric regression. We recommend the NCV estimate to describe the shape of the surface and REML + CL2 with Satterthwaite degrees of freedom for inference. For a surface that is zero everywhere, under dependent errors, these intervals exclude zero at about 5% of the grid points, and the NCV estimate is five to seven times closer to zero than REML's. For fitted means, intercepts and scalar effects, REML + CL2 intervals are close to nominal, except for binary intercepts and scalar effects and for misregistered counts. We state which defaults of refund the evidence supports.

Personalized federated learning can improve estimation for high-dimensional time series while preserving client-specific dynamics. However, parameter heterogeneity creates a fundamental ambiguity between shared and personalized structure. We introduce the principle of personalization diversity, under which genuinely personalized effects recur in only a limited fraction of clients. Based on this principle, we propose PerFeCT-VAR, Personalized Federated Vector Autoregression via Frequency-Capped Thresholding, which decomposes each client-specific transition matrix into shared low-rank dynamics, shared sparse links, and personalized sparse departures. A cross-client frequency cap yields a sharp shared–personalized separation threshold and motivates frequency-capped thresholding for federated estimation. Our theory further accounts for heterogeneous client distributions through a distributional compatibility condition and establishes joint linear convergence of the shared and personalized components. Under sufficient personalization diversity, the shared dynamics retain federated total-sample-size gains while the personalized components achieve client-level accuracy. Simulations and an application to multi-store retail revenue data demonstrate the predictive and interpretive benefits of the proposed framework.

High-dimensional mean testing is challenging under strong dependence and heavy tails, since classical Hotelling statistics are ill posed and regularized versions based on sample moments remain sensitive to outliers. This paper develops a robust regularized Hotelling framework for one-sample mean inference under elliptical distributions. The proposed HRST statistic combines a Hettmansperger-Randles spatial location estimator with the inverse of a trace-normalized shrinkage Tyler scatter matrix, and calibrates the resulting quadratic form through ridge-resolvent random-matrix theory rather than sparse covariance or precision-matrix assumptions. In the regime $p/n\to y\in(0,\infty)$, allowing a vanishing lower bulk edge and finite-rank diverging spikes, we establish feasible asymptotic normality for each shrinkage level, derive local-alternative distributions with explicit noncentrality parameters, and prove joint Gaussian limits over finite shrinkage grids. These results justify an adaptive Cauchy combination test across shrinkage levels. Simulations and a paired tumor-normal gene-expression study show that HRST maintains size and improves robustness under heavy-tailed distributions while retaining power under strong correlation.

Binary classification often requires controlling both type I and type II error rates at prespecified levels. When the class distributions overlap substantially, the two error limits may be incompatible if every observation must be classified. We study Neyman-Pearson classification with indecisions (NPI), which permits uncertain observations to remain unclassified and assesses selective type I and type II errors among definitive classifications within each class. For an independently trained score, we develop a model-free calibration procedure that combines all candidate threshold pairs passing calibration into a rule whose decision region contains that of every passing candidate. We couple the entire threshold search to a common martingale within each class and compute its finite-horizon crossing probability exactly, avoiding a multiplicity adjustment for the number of candidate pairs. The resulting classifier satisfies both population selective error limits simultaneously with user-specified high probability in finite samples, without assumptions on the accuracy of the fitted score. We bound the additional indecision due to calibration and evaluate the procedure in simulations and applications to recidivism and credit-default prediction.

We propose confidence intervals with near-nominal coverage of variance-covariance parameters in linear mixed models, even when interest or nuisance parameters are near the boundary, that is, when some random effects have small variances or are strongly correlated. In such settings, our simulations show standard Wald, likelihood ratio, and score intervals with nominal level 95% can have coverage as low as 45% and as high as 99%, while the proposed intervals are near nominal. The proposed intervals invert a modified profile score statistic based on the restricted likelihood, with nuisance parameters estimated on an extended parameter set. In a longitudinal study of socialization in children with autism, where the estimated correlation between the random intercept and slope is near minus one, our intervals show that children who start higher improve more slowly. In a genetic study of white blood cell count in mice, they distinguish the contributions of variant classes even though several intervals include zero. Compared to universal inference intervals, which are valid in finite samples, ours are 1.7 to 3.1 times narrower, and finite where one universal interval is not, at a forty-seventh of the computing time. We establish asymptotically correct coverage uniformly in the parameter under a condition known to hold with independent clusters and crossed random effects.

Spatial clustering increasingly involves complex objects, including probability distributions, networks, and structured matrices, for which conventional parametric likelihoods may be difficult to specify. We construct a generalized likelihood by exponentiating a cluster-specific Bregman loss on fixed representations of the responses, avoiding the need to specify a parametric sampling model for each response class. Each cluster is associated with an unknown representative, and the generalized likelihood measures within-cluster homogeneity through the Bregman divergence between the represented observations and that representative. A prior family matched to the Bregman geometry yields conjugate generalized posterior updates and permits exact marginalization of the cluster-specific representatives. The resulting collapsed scores combine within-cluster Bregman dispersion with uncertainty in the representatives, while a spanning-tree partition prior restricts posterior partition support to spatially contiguous clusters. This collapsed representation also leads to a tractable MCMC algorithm for posterior computation. The framework accommodates multiple complex-object response classes within a common inferential construction. We establish posterior consistency for the spatial partition and cluster-specific representatives under infill-domain asymptotics. Simulation studies with different complex-object responses demonstrate accurate recovery of spatial partitions and cluster-specific representatives. We further analyze racial-composition distributions across spatial units in Houston, illustrating geographically coherent clusters with distinct demographic compositions.

Imperfect measurements allow different underlying distributions to generate the same observations. Restrictions linking unobserved components can make it difficult to characterize all compatible distributions and compute the resulting range of parameter values. We ask how to construct computable parameter bounds and tighten them without additional data or stronger assumptions. We develop a general framework for partial identification under imperfect measurement based on distributionally robust optimization. We use measurement restrictions to construct a neighborhood around a benchmark distribution, ensuring that every compatible distribution lies within it. Optimizing over this neighborhood gives parameter bounds through a common procedure that can be adapted to different measurement problems. The choice of benchmark affects how informative these bounds are. We establish conditions under which changing it can tighten the bounds while preserving every parameter value consistent with the same data and assumptions. Applications to networks, missing outcomes, and regression illustrate how the framework tightens bounds and can recover the full range of compatible values.

We study pointwise estimation of a smooth binomial mixing density and characterize its minimax rates. Assuming that the mixing density is $s$-Hölder smooth, we first consider the homogeneous setting with a common number of binomial trials $t$. We derive matching lower and upper bounds that reveal three regimes depending on the relative sizes of the sample size $n$ and the number of trials $t$. When $t$ is small, the finite number of trials creates an identification barrier that persists regardless of sample size; at intermediate $t$, more data reduce sampling uncertainty, while the number of trials limits how finely the mixing density can be recovered; and when $t$ is sufficiently large, the usual nonparametric density estimation rate is recovered. Our lower bounds exploit the polynomial structure of the binomial mixture model, while attainability is achieved by a suitably regularized orthogonal series estimator. We then extend the analysis to heterogeneous trial parameters, where the minimax rate depends on the full trial profile through the effective sample size available at each polynomial degree. These results characterize the fundamental statistical limits of learning smooth populations of binomial probabilities under both homogeneous and heterogeneous trials.

Matching is widely used to mimic randomized experiments by forming matched sets in which treated and control units differ only randomly with respect to observed covariates. However, when the study population consists of interconnected units from a single network or a small number of networks, matching solely on observed covariates may produce matched units that are more closely connected in the network than would occur by chance. This increased network proximity within matched sets can induce spurious associations between treatment and outcome when both variables exhibit similar autocorrelation patterns on the network. To reduce spurious associations while preserving the validity of causal comparisons, we propose a new matching method that matches units with similar covariates subject to additional network-proximity constraints. For post-matching inference, we propose a randomization-based procedure for testing the sharp null hypothesis of no causal effect. The inference uses the asymptotic normal approximation and accommodates statistical dependence among test statistics obtained from each matched set without requiring explicit specifications of their correlation structures. We demonstrate the validity and utility of the proposed methods through simulation studies and apply them to the Framingham Heart Study. The matching method and subsequent inference procedure are implemented in the R package netmatchRI.

We consider the estimation of the invariant distribution $\pi$ of a Markov kernel on a discrete state space (or of the mean value $\pi(f)$ of a functional $f$ under $\pi$) when solving the invariance equation is not feasible. Look-ahead estimators exploit the knowledge of the kernel by propagating the empirical occupation measure through $k$ steps of the chain. When the kernel is known, as in Markov chain Monte Carlo, we compare their mean squared error with that of the empirical estimator, for independent or Markovian data, stationary or not. Our main result establishes that, for Markovian data, the one-step estimator always improves on the empirical one asymptotically, and that the improvement grows with $k$ for reversible kernels. A non-reversible counterexample shows that reversibility cannot simply be dropped. In the statistical setting where the kernel belongs to a parametric family with unknown parameter, we introduce plug-in look-ahead estimators, which perform the look-ahead step with an estimated kernel. When the parameter is estimated at the parametric rate, we give a necessary and sufficient condition under which the plug-in look-ahead estimator outperforms the empirical one. The theoretical results are illustrated by extensive simulation studies, which also shed further light on the behaviour of the estimators at hand. Both settings are finally applied to real data arising from geographical mobility problems in the Lyon metropolitan area, France.

The LSJM R package provides a series of functions to estimate models based on location-scale linear mixed model and location-scale joint model theory. For longitudinal data, it includes the estimation of location-scale mixed models with subject-specific residual variance that can be time-dependent or distinguish within-visit and between-visits variability. When survival data is considered, implemented joint models handle competing risks as well as semi-competing risks possibly with interval-censored data through the use of a illness-death sub-model. They incorporate flexible dependence structure between the longitudinal marker and the times-to-events through the current value, the slope and/or the residual variability of the marker. Maximum likelihood estimators are obtained using the Marquardt-Levenberg algorithm with stringent convergence criteria. The package also provides several post-fit functions including goodness-of-fit analyses and individual dynamic predictions. Some examples on real data allows to illustrates the use of this package.

Motivated by a real-life problem, namely, inferring genomic regions affected by selection, we studied some properties of a two-step method that performs this task based on the FST. This approach exhibits computational simplicity as an appealing feature. Step one involves fitting a conjugate Bayesian model to find point estimates of FST and selecting the markers to be considered in step two, where they are grouped using such estimates. Subjective priors are frequently used in this problem and large hyperparameters may be encountered; hence, there is interest in addressing their impact. Preliminary analyses suggested that, the number of markers selected in step one tends to zero when the hyperparameters have the same value and tend to infinity, or are different but tend to infinity and their difference tends to zero. Thus, we studied how hyperparameters affect the posterior mean, variance, and predictive density of allele frequencies, and the posterior distribution of FST, formally. Moreover, an empirical approach based on simulated and real data was performed. Our main findings are summarized in two theorems that set the asymptotic behavior of the method, providing a foundation to understand the aforementioned phenomena. As to the empirical studies, the groups created in the second step, the estimated posterior means and the number of selected markers showed high sensitivity to hyperparameters. Albeit the mathematical intractability of the posterior distribution of FST, our results revealed asymptotic properties explaining two features of the method that lacked a formal elucidation.

Aggregated relational data replaces a fully observed network with a compressed version, where survey respondents are asked questions such as “How many of your friends belong to group $G$?”. We consequently only observe a matrix $R = OA$ of the adjacency matrix $A$ with a fixed design $O$. In this paper we understand the effects of this compression on community detection under degree-corrected stochastic block models. First, we give guarantees for when it is possible to recover the full network parameters, and then only just the communities under conditions on a diagonal corrected version of $\mathbb{E}[R]$, which can be used to inform experimental design. Secondly, we show that spectral clustering with single linkage exactly recovers communities when $n\rho_n\gtrsim K^2+K\log n$ when there are $K$ balanced communities under “uniformly spread” or near-orthogonal designs, and characterize more generally how the design influences rates of convergence. We also introduce second-order “friends-of-friends” queries and develop the corresponding theory.

Identifying patient subgroups with heterogeneous treatment effects is central to precision medicine, yet existing approaches face a tension. Black-box methods estimate individualized effects well but yield no interpretable subgroups, while tree-based methods produce explicit partitions but are unstable and less accurate in moderate-sample trials. We propose CausalSynthTree, a co-supervised method in which black-box causal estimators guide construction of an interpretable tree. It partitions the covariate space into cells, fits cell-wise conditional average treatment effect models from both observed data and synthetic data labeled by black-box causal teachers, and synthesizes them into a tree through a treatment-effect disparity criterion. Each leaf carries a sparse linear model, so the tree path defines a subgroup and the leaf model shows which covariates modify the effect. Honest inference for leaf-wise effects remains valid despite co-supervised augmentation and data-driven subgroup selection. In simulations it closes much of the gap between the two families. From moderate sample sizes onward it is more accurate than the black-box learners whenever the effect is linear or has a single threshold, and trails them only on a finer partition than the tree it grows. It controls spurious splits and attains nominal coverage throughout. In the ACTG175 HIV trial, where standard interaction tests detect no effect modification, the method reports no subgroups while the competing tree methods report several. In an adjuvant colon cancer trial it reports a borderline covariate in part of the resamples, together with a within-region age gradient that is consistent with an independent test.

Partially observed multivariate time series pose two coupled challenges: recovering missing measurements and forecasting future trajectories when observation patterns are informative. Existing methods often assume ignorable missingness or use the mask only as an auxiliary input. We propose LIFTS (Latent-driven Imputation and Forecasting for partly observed Time Series), a probabilistic framework in which a latent process jointly governs the complete measurement trajectory and observation mask, allowing observation patterns to inform latent-state inference for both forecasting and imputation. LIFTS combines flexible neural-network parameterizations with a forward-backward architecture: the forward pass performs filtering and autoregressive forecasting, while the backward pass performs smoothing and conditional multiple imputation. The model and inference components are trained jointly using a masked energy score and self-masking. We establish forward-backward distributional representations for the general framework and, for an explicit structural subclass, provide full-law identification and consistency guarantees. Simulations and an application to PhysioNet 2019 show that LIFTS provides accurate point predictions, improved distributional accuracy, and well-calibrated predictive intervals; it outperforms competing methods particularly when training data are limited and has substantially faster inference than the diffusion-based baseline.