Papers

We propose an optimal algorithm for estimating conditional average treatment effects (CATEs) when response functions lie in a reproducing kernel Hilbert space (RKHS). We study settings in which the contrast function is structurally simpler than the nuisance functions: (i) it lies in a lower-complexity RKHS with faster eigenvalue decay, (ii) it satisfies a source condition relative to the nuisance kernel, or (iii) it depends on a known low-dimensional covariate representation. We develop a unified two-stage kernel ridge regression (KRR) method that attains minimax rates governed by the complexity of the contrast function rather than the nuisance class, in terms of both sample size and overlap. We also show that a simple model-selection step over candidate contrast spaces and regularization levels yields an oracle inequality, enabling adaptation to unknown CATE regularity.

Transformer models contain substantial internal redundancy arising from coordinate-dependent representations and continuous symmetries, in model space and in head space, respectively. While recent approaches address this by explicitly breaking symmetry, we propose a complementary framework based on symmetry reduction. We reformulate representations, attention mechanisms, and optimization dynamics in terms of invariant relational quantities, eliminating redundant degrees of freedom by construction. This perspective yields architectures that operate directly on relational structures, providing a principled geometric framework for reducing parameter redundancy and analyzing optimization.

Many fairness goals are defined at a population level that misaligns with siloed data collection, which remains unsharable due to privacy regulations. Horizontal federated learning (FL) enables collaborative modeling across clients with aligned features without sharing raw data. We study federated auditing of demographic parity through score distributions, measuring disparity as a Wasserstein–Frechet variance between sensitive-group score laws, and expressing the population metric in federated form that makes explicit how silo-specific selection drives local-global mismatch. For the squared Wasserstein distance, we prove an ANOVA-style decomposition that separates (i) selection-induced mixture effects from (ii) cross-silo heterogeneity, yielding tight bounds linking local and global metrics. We then propose a one-shot, communication-efficient protocol in which each silo shares only group counts and a quantile summary of its local score distributions, enabling the server to estimate global disparity and its decomposition, with $O(1/k)$ discretization bias ($k$ quantiles) and finite-sample guarantees. Experiments on synthetic data and COMPAS show that a few dozen quantiles suffice to recover global disparity and diagnose its sources.

Despite the widespread use of boosting in structured prediction, a general theoretical understanding of aggregation beyond scalar prediction remains incomplete. We study vector-valued prediction under a target divergence and identify a geometric stability property under which aggregation amplifies weak guarantees into strong ones. We formalize this property as $(\alpha,\beta)$-stability by geometric median and show how it supports a boosting framework based on exponential reweighting and geometric-median aggregation. For vector-valued prediction, we characterize this stability property under several natural divergences: $\ell_1$ and $\ell_2$ distances for unconstrained vector-valued prediction, and TV, Hellinger, and KL for density estimation over finite probability vectors. Building on these results, we propose a generic boosting framework \geomedboost. Under a weak learner condition and $(\alpha,\beta)$-stability, we obtain exponential decay of the empirical divergence error, which then yields population guarantees through a generalization bound.

To provide a comprehensive summary of the tail distribution, the expected shortfall is defined as the average over the tail above (or below) a certain quantile of the distribution. The expected shortfall regression captures the heterogeneous covariate-response relationship and describes the covariate effects on the tail of the response distribution. Based on a critical observation that the superquantile regression from the operations research literature does not coincide with the expected shortfall regression, we propose and validate a novel optimization-based approach for the linear expected shortfall regression, without additional assumptions on the conditional quantile models. While the proposed loss function is implicitly defined, we provide a prototype implementation of the proposed approach with some initial expected shortfall estimators based on binning techniques. With practically feasible initial estimators, we establish the consistency and the asymptotic normality of the proposed estimator. The proposed approach achieves heterogeneity-adaptive weights and therefore often offers efficiency gain over existing linear expected shortfall regression approaches in the literature, as demonstrated through simulation studies.

We consider the orthogonalisation of the signature of a stochastic process as the analogue of orthogonal polynomials on path-space. Under an infinite radius of convergence assumption, we prove density of linear functions on the signature in $L^p$ functions on grouplike elements, making it possible to represent a square-integrable function on (rough) paths as an $L^2$-convergent series. By viewing the shuffle algebra as commutative polynomials on the free Lie algebra, we revisit much of the theory of classical orthogonal polynomials in several variables, such as the recurrence relation and Favard's theorem. Finally, we restrict our attention to the case of Brownian motion with and without drift, and prove that dimension-independent orthogonal signature exists with drift but not without. We end with numerical examples of how orthogonal signature polynomials of Brownian motion can be applied for the approximation of functions on paths sampled from the Wiener measure.

Latent Dirichlet Allocation (LDA) is a foundational model for discovering latent thematic structure in discrete data, but its Dirichlet prior cannot represent the rich correlations and hierarchical relationships often present among topics. We introduce the framework of Latent Dirichlet-Tree Allocation (LDTA), a generalization of LDA that replaces the Dirichlet prior with an arbitrary Dirichlet-Tree (DT) distribution. LDTA preserves LDA's generative structure but enables expressive, tree-structured priors over topic proportions. To perform inference, we develop universal mean-field variational inference and Expectation Propagation, providing tractable updates for all DT. We reveal the vectorized nature of the two inference methods through theoretical development, and perform fully vectorized, GPU-accelerated implementations. The resulting framework substantially expands the modeling capacity of LDA while maintaining scalability and computational efficiency.

Evaluating joint probabilities of potential outcomes and observed variables, and their linear combinations, is a fundamental challenge in causal inference. This paper addresses the bounding and identification of these probabilities in settings with discrete treatment and discrete ordinal outcome. We propose new families of monotonicity assumptions and formulate the bounding problem as a linear programming problem. We further introduce a new monotonicity assumption specifically to achieve identification. Finally, we present numerical experiments to validate our methods and demonstrate their application using real-world datasets.

MINFLUX (Minimal Photon Flux) is a single-molecule imaging technique capable of resolving fluorophores at a precision of <5 nm. Interpretation of the point patterns generated by this technique presents challenges due to variable emitter density, incomplete bio-labelling of target molecules and their detection, error prone measurement processes, and the presence of spurious (non-structure associated) fluorescent detections. Together, these challenges ensure structural inferences from single-molecule imaging datasets are non-trivial in the absence of strong a priori information, for all but the smallest of point patterns. In addition, current methods often require subjective parameter tuning and presuppose known structural templates, limiting reference-free discovery. We present a statistically grounded, end-to-end analysis framework. Focusing on MINFLUX derived datasets and leveraging Bayesian and spatial statistical methods, a pipeline is presented that demonstrates 1) uncertainty aware clustering of measurements into emitter groups that performs better than current gold standards, 2) rapid identification of molecular structure supergroups, and 3) reconstruction of repeating structures within the dataset without substantial prior knowledge. This pipeline is demonstrated using simulated and real MINFLUX datasets, where emitter clustering and centre detection maintain high performance (emitter subset assignment accuracy > 0.75) across all conditions evaluated, while structural inference achieves reliable discrimination (F1 approx. 0.9) at high labelling efficiency. Template-free reconstruction of Nup96 and DNA-Origami 3x3 grids are achieved.

For approximating a target distribution given only its unnormalized log-density, stochastic gradient-based variational inference (VI) algorithms are a popular approach. For example, Wasserstein VI (WVI) and black-box VI (BBVI) perform gradient descent in measure space (Bures-Wasserstein space) and parameter space, respectively. Previously, for the Gaussian variational family, convergence guarantees for WVI have shown superiority over existing results for black-box VI with the reparametrization gradient, suggesting the measure space approach might provide some unique benefits. In this work, however, we close this gap by obtaining identical state-of-the-art iteration complexity guarantees for both. In particular, we identify that WVI's superiority stems from the specific gradient estimator it uses, which BBVI can also leverage with minor modifications. The estimator in question is usually associated with Price's theorem and utilizes second-order information (Hessians) of the target log-density. We will refer to this as Price's gradient. On the flip side, WVI can be made more widely applicable by using the reparametrization gradient, which requires only gradients of the log-density. We empirically demonstrate that the use of Price's gradient is the major source of performance improvement.

In this paper, we develop a Bayesian calendar-time survival model motivated by infectious disease prevention studies occurring during an epidemic, when the risk of infection can change rapidly as the epidemic curve shifts. For studies in which a biomarker is the predictor of interest, we include the option to estimate a threshold of protection for the biomarker. If the intervention is hypothesized to have different associations with several circulating viral variants, or if the infectiousness of the dominant variant(s) changes over the course of the study, we treat infection from different variants as competing risks. We also introduce a novel method for incorporating existing epidemic curve estimates into an informative prior for the baseline hazard function, enabling estimation of the intervention's association with infection risk during periods of calendar time with minimal follow-up in one or more comparator groups. We demonstrate the strengths of this method via simulations, and we apply it to data from an observational COVID-19 vaccine study.

Efficiently routing queries to the optimal large language model (LLM) is crucial for optimizing the cost-performance trade-off in multi-model systems. However, most existing routers rely on single-model selection, making them susceptible to misrouting. In this work, we formulate LLM routing as the $\alpha$-VOR problem to minimize expected set size while controlling the misrouting risk, and propose a novel method – RACER, extending base routers to output model sets that can be subsequently aggregated for improved output. In particular, RACER constructs nested model sets via augmented scoring and utilizes finite-sample concentration bounds to calibrate a threshold that allows for both variable set sizes and abstention. We theoretically prove that RACER achieves rigorous distribution-free risk control on unseen test data in a post-hoc and model-agnostic manner. Extensive experiments verify our theoretical guarantees and demonstrate that RACER consistently enhances downstream accuracy across a wide range of benchmarks.

Despite the widespread use of ordinal measures in HCI, such as Likert-items, there is little consensus among HCI researchers on the statistical methods used for analysing such data. Both parametric and non-parametric methods have been extensively used within the discipline, with limited reflection on their assumptions and appropriateness for such analyses. In this paper, we examine recent HCI works that report statistical analyses of ordinal measures. We highlight prevalent methods used, discuss their limitations and spotlight key assumptions and oversights that diminish the insights drawn from these methods. Finally, we champion and detail the use of cumulative link (mixed) models (CLM/CLMM) for analysing ordinal data. Further, we provide practical worked examples of applying CLM/CLMMs using R to published open-sourced datasets. This work contributes towards a better understanding of the statistical methods used to analyse ordinal data in HCI and helps to consolidate practices for future work.

In meta analysis, multiple hypothesis testing and many other methods, p-values are utilized as inputs and assumed to be uniformly distributed over the unit interval under the null hypotheses. If data used to generate p-values have discrete distributions then either natural, mid- or randomized p-values are typically utilized. Natural and mid-p-values can allow for valid, albeit conservative, downstream methods since under the null hypothesis they are dominated by uniform distributions in the stochastic and convex order, respectively. Randomized p-values need not lead to conservative procedures since they permit a uniform distributions under the null hypotheses through the generation of independent auxiliary variates. However, the auxiliary variates necessarily add variation to procedures. This manuscript introduces and studies “minimally discrete” (MD) natural p-values, MD mid-p-values and “minimally randomized” (MR) p-values. It is shown that MD p-values dominate their non-MD counterparts in the stochastic and convex order, and hence lead to less conservative, yet still valid, downstream methods. Likewise, MR p-values dominate their non-MR counterparts in that they are still uniformly distributed under the null hypotheses, but the added variation attributable to the independently generated auxiliary variate is smaller. It is anticipated that results here will facilitate the construction of new meta-analysis and multiple testing methods via more efficient p-value construction, and facilitate theoretical study of existing and new methods by establishing gold standards for addressing the unavoidable detrimental “discreteness effect”.

This paper develops a hybrid likelihood (HL) method based on a compromise between parametric and nonparametric likelihoods. Consider the setting of a parametric model for the distribution of an observation $Y$ with parameter $\theta$. Suppose there is also an estimating function $m(\cdot,\mu)$ identifying another parameter $\mu$ via $E\,m(Y,\mu)=0$, at the outset defined independently of the parametric model. To borrow strength from the parametric model while obtaining a degree of robustness from the empirical likelihood method, we formulate inference about $\theta$ in terms of the hybrid likelihood function $H_n(\theta)=L_n(\theta)^{1-a}R_n(\mu(\theta))^a$. Here $a\in[0,1)$ represents the extent of the compromise, $L_n$ is the ordinary parametric likelihood for $\theta$, $R_n$ is the empirical likelihood function, and $\mu$ is considered through the lens of the parametric model. We establish asymptotic normality of the corresponding HL estimator and a version of the Wilks theorem. We also examine extensions of these results under misspecification of the parametric model, and propose methods for selecting the balance parameter $a$.

Over the last century, the adoption of novel scientific methods for conducting the U.S. census has been met with wide-ranging receptions. Some methods were quietly embraced, while others sparked decades-long controversies. What accounts for these differences? We argue that controversies emerge from $\textit{arrangements of statistical imaginaries}$, putting into tension divergent visions of the census. To analyze these dynamics, we compare reactions to two methods designed to improve data accuracy (imputation and adjustment) and two methods designed to protect confidentiality (swapping and differential privacy), offering insight into how each method reconfigures stakeholder orientations and rhetorical claims. These cases allow us to reflect on how technocratic efforts to improve accuracy and confidentiality can strengthen – or erode – trust in data. Our analysis shows how the credibility of the Census Bureau and its data stem not just from empirical evaluations of quantification, but also from how statistical imaginaries are contested and stabilized.

We present balnet, an R package for scalable pathwise estimation of covariate balancing propensity scores via logistic covariate balancing loss functions. Regularization paths are computed with Yang and Hastie (2024)'s generic elastic net solver, supporting convex losses with non-smooth penalties, as well as group penalties and feature-specific penalty factors. For lasso penalization, balnet computes a regularization path of balancing weights from the largest observed covariate imbalance to a user-specified fraction of this maximum. We illustrate the method with an application to spatial pixel-level balancing for constructing synthetic control weights for the average treatment effect on the treated, using satellite data on wildfires.

Machine-generated probability predictions are essential in modern classification tasks such as image classification. A model is well calibrated when its predicted probabilities correspond to observed event frequencies. Despite the need for multicategory recalibration methods, existing methods are limited to (i) comparing calibration between two or more models rather than directly assessing the calibration of a single model, (ii) requiring under-the-hood model access, e.g., accessing logit-scale predictions within the layers of a neural network, and (iii) providing output which is difficult for human analysts to understand. To overcome (i)-(iii), we propose Multicategory Linear Log Odds (MCLLO) recalibration, which (i) includes a likelihood ratio hypothesis test to assess calibration, (ii) does not require under-the-hood access to models and is thus applicable on a wide range of classification problems, and (iii) can be easily interpreted. We demonstrate the effectiveness of the MCLLO method through simulations and three real-world case studies involving image classification via convolutional neural network, obesity analysis via random forest, and ecology via regression modeling. We compare MCLLO to four comparator recalibration techniques utilizing both our hypothesis test and the existing calibration metric Expected Calibration Error to show that our method works well alone and in concert with other methods.

Environmental policy evaluation frequently requires thoughtful consideration of space and time in causal inference. We use novel statistical methods to analyze the causal effect of land concessions on deforestation rates in Cambodia. Standard approaches, such as difference-in-differences regression, effectively address spatiotemporally-correlated treatments under some conditions, but they are limited in their ability to account for unobserved confounders affecting both treatment and outcome. Double Spatial Regression (DSR) is an approach that uses double machine learning to address these scenarios. DSR resolves the confounding variables for both treatment and outcome, comparing the residuals to estimate treatment effectiveness. We improve upon DSR by considering time in our analysis of policy interventions with spatial effects. We conduct a large-scale simulation study using Bayesian Additive Regression Trees (BART) with spatial embeddings and find that, under certain conditions, our DSR model outperforms standard approaches for addressing unobserved spatial confounding. We then apply our method to evaluate the policy impacts of land concessions on deforestation in Cambodia.

Biological neural networks (like the hippocampus) can internally generate "replay" resembling stimulus-driven activity. Recent computational models of replay use noisy recurrent neural networks (RNNs) trained to path-integrate. Replay in these networks has been described as Langevin sampling, but new modifiers of noisy RNN replay have surpassed this description. We re-examine noisy RNN replay as sampling to understand or improve it in three ways: (1) Under simple assumptions, we prove that the gradients replay activity should follow are time-varying and difficult to estimate, but readily motivate the use of hidden state leakage in RNNs for replay. (2) We confirm that hidden state adaptation (negative feedback) encourages exploration in replay, but show that it incurs non-Markov sampling that also slows replay. (3) We propose the first model of temporally compressed replay in noisy path-integrating RNNs through hidden state momentum, connect it to underdamped Langevin sampling, and show that, together with adaptation, it counters slowness while maintaining exploration. We verify our findings via path-integration of 2D triangular and T-maze paths and of high-dimensional paths of synthetic rat place cell activity.

We propose PRISM-FCP (Partial shaRing and robust calIbration with Statistical Margins for Federated Conformal Prediction), a communication-efficient Byzantine-robust federated conformal prediction framework that uses partial model sharing to mitigate stochastic model-poisoning attacks during training and histogram-based filtering to mitigate adversarial calibration submissions. Existing robust FCP approaches primarily address adversarial behavior during calibration, leaving training-stage poisoning to separate robust-training mechanisms. PRISM-FCP instead considers the coupling between the two stages. During training, clients partially share updates by transmitting only $M$ of $D$ parameters per round. Under the stated stochastic attack model, this attenuates the expected energy of each Byzantine client's perturbation contribution to the aggregate by a factor of $M/D$ relative to full sharing. When this benefit outweighs the optimization slowdown caused by partial updates, it can reduce training error and improve interval efficiency. During calibration, the server uses client-provided characterization vectors to filter suspected Byzantine clients before estimating the conformal quantile from the retained clients. Experiments on synthetic benchmarks and the UCI Superconductivity and YearPredictionMSD datasets, including Gaussian, ALIE, and sign-flipping training-stage attacks, demonstrate near-nominal empirical coverage and favorable communication–performance tradeoffs in the studied settings.

Generalized pairwise comparisons summarize treatment effects in randomized experiments with ordinal, multivariate, or prioritized outcomes. We study optimal covariate adjustment under complete randomization, with potential outcomes and covariates fixed and treatment assignment as the sole source of randomness. We derive the exact randomization variance of the unadjusted estimator and characterize linear calibration adjustments using the between-arm difference in covariate means. The asymptotic variance is a strictly convex quadratic in the adjustment coefficient, yielding a closed-form optimum. The resulting estimator is asymptotically no less efficient than the unadjusted estimator, with an explicit precision gain that clarifies which covariates improve precision. However, natural pair-level and single-arm least squares adjustments generally use weights that differ from the allocation-dependent optimum and can substantially reduce efficiency. We also propose a within-arm variance estimator for asymptotically conservative inference, whose conservatism arises from unidentifiable joint variation in a unit's two potential scores, and show that residualizing each arm on its own projection coefficient saves the discrepancy that governs the efficiency loss of conventional adjustments. In the special case where the comparison function is the outcome difference, our proposed covariate-adjusted estimator is first-order equivalent to the fully interacted regression adjustment of Lin (Ann. Appl. Stat., 2013, 7, 295-318).

Multimorbidity in older adults is common, heterogeneous, and highly dynamic, and it is strongly associated with disability and increased healthcare utilization. However, existing approaches to studying multimorbidity trajectories are largely descriptive or rely on discrete-time models, which struggle to handle irregular observation intervals and right-censoring. We developed a continuous-time hidden multistate modeling framework to capture transitions among latent multimorbidity patterns while accounting for interval censoring and misclassification. A simulation study compared alternative model specifications under varying sample sizes and follow-up schemes, and the best-performing specification was applied to longitudinal data from the Swedish National study on Aging and Care-Kungsholmen (SNAC-K), including 2 716 multimorbid participants followed for up to 18 years. Simulation results showed that hidden multistate models substantially reduced bias in transition hazard estimates compared to non-hidden models, with fully time-inhomogeneous models outperforming piecewise approximations. Application to SNAC-K confirmed the feasibility and practical utility of this framework, enabling identification of risk factors for accelerated progression toward complex multimorbidity and revealing a gradient of mortality risk across patterns. Continuous-time hidden multistate models provide a robust alternative to traditional approaches, supporting individualized predictions and informing targeted interventions and secondary prevention strategies for multimorbidity in aging populations.

Maximum likelihood prediction (MLP) is a core task at the heart of modern large language models. Here, we study a quantum version of this task for a simplified data model consisting of independent and identically distributed samples, as a first step. The quantum maximum likelihood predictor (QMLP) is obtained by embedding of empirical probability distributions into quantum states and performing a minimization of quantum relative entropy over a given class of states. We derive non-asymptotic performance guarantees for QMLP in terms of convergence rates and concentration inequalities, both in trace norm and quantum relative entropy. Our approach provides a unified framework to handle MLP within both classical and quantum LLMs. We also consider the related problem of quantum information projection and generalize the quantum Pythagorean theorem to mixture families specified by possibly non-self-adjoint linear constraints. We further show that the Pythagorean inequality continues to hold in the infinite-dimensional setting whenever the convex information-projection problem attains a finite minimum.

Understanding how the composition of guest origin markets evolves over time is critical for destination marketing organizations, hospitality businesses, and tourism planners. We develop and apply Bayesian Dirichlet autoregressive moving average (BDARMA) models to forecast the compositional dynamics of guest origin market shares using proprietary Airbnb booking data spanning 2017–2025 across four major destination regions. Our analysis reveals substantial pandemic-induced structural breaks in origin composition, with heterogeneous recovery patterns across markets. In our analysis, the BDARMA framework achieves the lowest forecast error for EMEA and competitive performance across destination regions, outperforming standard benchmarks including naïve forecasts, exponential smoothing, and SARIMA on log-ratio transformed data in compositionally complex markets. For EMEA destinations, BDARMA achieves 27% lower forecast error than naïve methods ($p < 0.001$), with the greatest gains where multiple origin markets compete in the 5-25% share range. By modeling compositions directly on the simplex with a Dirichlet likelihood and incorporating seasonal variation in both mean and precision parameters, our approach produces coherent forecasts that respect the unit-sum constraint while capturing complex temporal dependencies. The methodology provides destination stakeholders with probabilistic forecasts of source market shares, enabling more informed strategic planning for marketing resource allocation, infrastructure investment, and crisis response.

We consider Bayesian inverse problems arising in data assimilation for dynamical systems governed by partial and stochastic partial differential equations. The space-time dependent field is inferred jointly with static parameters of the prior and likelihood densities. Particular emphasis is placed on the hyperparameter controlling the prior smoothness and regularity, which is critical in ensuring well-posedness, shaping posterior structure, and determining predictive uncertainty. Commonly it is assumed to be known and fixed a priori; however in this paper we will adopt a hierarchical Bayesian framework in which smoothness and other hyperparameters are treated as unknown and assigned hyperpriors. Posterior inference is performed using Metropolis-within-Gibbs sampling suitable to high dimensions, for which hyperparameter estimation involves little computational overhead. The methodology is demonstrated on inverse problems for the Navier-Stokes equations and the stochastic advection-diffusion equation, under sparse and dense observation regimes, using Gaussian priors with different covariance structure. Numerical results show that jointly estimating the smoothness substantially reduces the errors in uncertainty quantification and parameter estimation induced by smoothness misspecification, by achieving performance comparable to scenarios in which the true smoothness is known.

This work studies heterogeneous Multi-Objective Reinforcement Learning (MORL), where objectives can differ sharply in temporal frequency. Such heterogeneity allows dense objectives to dominate learning, while sparse long-horizon rewards receive weak credit assignment, leading to poor sample efficiency. We propose a Parallel Reward Integration with Symmetry (PRISM) algorithm that enforces reflectional symmetry as an inductive bias in aligning reward channels. PRISM introduces ReSymNet, a theory-motivated model that reconciles temporal-frequency mismatches across objectives, using residual blocks to learn a scaled opportunity value that accelerates exploration while preserving the optimal policy. We also propose SymReg, a reflectional equivariance regulariser that enforces agent mirroring and constrains policy search to a reflection-equivariant subspace. This restriction provably reduces hypothesis complexity and improves generalisation. Across MuJoCo benchmarks, PRISM consistently outperforms both a sparse-reward baseline and an oracle trained with full dense rewards, improving Pareto coverage and distributional balance: it achieves hypervolume gains exceeding 100% over the baseline and up to 32% over the oracle. The code is at \href{https://github.com/EVIEHub/PRISM}{https://github.com/EVIEHub/PRISM}.

In this paper, we present novel methodologies that incorporate auxiliary variables for multiple hypotheses testing related to the main point of interest while effectively controlling the false discovery rate. When dealing with multiple tests concerning the primary variable of interest, researchers can use auxiliary variables to set preconditions for the significance of primary variables, thereby enhancing test efficacy. Depending on the auxiliary variable's role, we propose two approaches: one terminates testing of the primary variable if it does not meet predefined conditions, and the other adjusts the evaluation criteria based on the auxiliary variable. Employing the copula method, we elucidate the dependence between the auxiliary and primary variables by deriving their joint distribution from individual marginal distributions.Our numerical studies, compared with existing methods, demonstrate that the proposed methodologies effectively control the FDR and yield greater statistical power than previous approaches solely based on the primary variable. As an illustrative example, we apply our methods to the Set4$\Delta$ mutant dataset. Our findings highlight the distinctions between our methodologies and traditional approaches, emphasising the potential advantages of our methods in introducing the auxiliary variable for selecting more genes.

Keeping pace with rapidly evolving technology is a key challenge in teaching statistics. To equip students with essential skills for the modern workplace, educators must integrate relevant technologies into the statistical curriculum where possible. University-level statistics education has experienced substantial technological change, particularly in the tools and practices that underpin teaching and learning. Statistical programming has become central to many courses, with R widely used and Python increasingly incorporated into statistics and data analytics programmes. Additionally, coding practices, database management, and machine learning now feature within some statistics curricula. Looking ahead, we anticipate a growing emphasis on artificial intelligence (AI), particularly the pedagogical implications of generative AI tools such as ChatGPT. In this article, we explore these technological developments and discuss strategies for their integration into contemporary statistics education.

We propose an online false discovery rate (FDR) controlling method based on conditional local FDR (LIS), designed for infectious disease datasets that are discrete and exhibit complex dependencies. Unlike existing online FDR methods, which often assume independence or suffer from low statistical power in dependent settings, our approach effectively controls FDR while maintaining high detection power in realistic epidemic scenarios. For disease modeling, we establish a Dynamic Bayesian Network (DBN) structure within the Susceptible-Infected-Susceptible (SIS) model, a widely used epidemiological framework for infectious diseases. Our method requires no additional tuning parameters apart from the width of the sliding window, making it practical for real-time disease monitoring. From a statistical perspective, we prove that our method ensures valid FDR control under stationary and ergodic dependencies, extending online hypothesis testing to a broader range of dependent and discrete datasets. Additionally, our method achieves higher statistical power than existing approaches by leveraging LIS, which has been shown to be more powerful than traditional $p$-value-based methods. We validate our method through extensive simulations and real-world applications, including the analysis of infectious disease incidence data. Our results demonstrate that the proposed approach outperforms existing methods by achieving higher detection power while maintaining rigorous FDR control.