Papers

We introduce a similarity representation framework for robust heritability estimation in Genome-Wide Association Studies (GWAS). This problem parallels the signal-to-noise ratio estimation problem in linear models with a large number of predictors. Traditional fixed- and random-effects methods for heritability estimation often impose restrictive assumptions on regression coefficients or the design (genotype) matrix. These assumptions are usually violated by the heterogeneous effects of genetic variants (regression coefficients) that depend on the genotype distribution and the correlation among genotypes due to linkage disequilibrium. This leads to the non-robust estimation of heritability in practice. To overcome these limitations, we propose a SiMILarity rEpresentation method (SMILE) which models the relationship between the outcome similarity and genetic similarity through Gram matrices. SMILE represents genetic similarity using a weighted Gram matrix of genotypes, where a data-dependent weight matrix is used to disentangle the heterogeneous variant effects from the genotype distribution. SMILE includes the classical random-effects model as a special case and improves the fixed-effects model by not requiring accurate estimation of the precision matrix or the regression coefficients. We develop a scalable implementation for efficient analysis of large biobank GWAS data. Extensive simulations and the analysis of the UK Biobank data demonstrate the robustness of the proposed method over the existing methods across a range of genetic architectures, and show that SMILE provides a versatile approach for heritability estimation.

Dynamic causal modeling (DCM) estimates directed effective connectivity, while parametric empirical Bayes (PEB) supports group inference using subject-specific posterior summaries. Standard PEB relies on Gaussian models and continuous shrinkage, making it sensitive to atypical estimates and unable to distinguish negligible from nonzero effects. We develop a robust and sparse group-DCM extension combining a Student-t likelihood with a spike-and-slab prior using a nonlocal product-moment (pMOM) slab. A normal–gamma representation of the Student-t distribution yields weights that downweight atypical subject–parameter combinations. The pMOM slab vanishes at zero, sharpening coefficient selection and yielding inclusion probabilities. We propagate first-level posterior uncertainty through block-covariance pre-whitening and estimate the model using an EM–ReML algorithm that updates weights, inclusion probabilities, effects, and variance components. A simulation study showed that Student-t weighting provided the main protection against contamination, the nonlocal prior contributed most under strong sparsity, and their combination was most beneficial when contamination and sparsity occurred together. In an openly shared mixed-gambles fMRI application, we found mild heterogeneity, concentrated inclusion support on three intrinsic self-connections, and close directional agreement with standard SPM-PEB. Our framework therefore adds interpretable element-level robustness diagnostics and sparse coefficient selection to hierarchical DCM while retaining first-level uncertainty.

Conditional dimension reduction (cDR) learns low-dimensional latent coordinates while accounting for observed covariates that represent known sources of variation in the data. Conditional Multidimensional Scaling (cMDS) is a cDR technique that works directly with dissimilarity data. Its standard squared-stress formulation, however, is sensitive to contaminated dissimilarity, since outliers can dominate the objective and distort the learned configuration. We proposed Robust Conditional Multidimensional Scaling (rcMDS) by replacing the squared-stress criterion with a Fair M-estimation objective. We developed a reweighted conditional SMACOF algorithm to optimize this objective. The proposed algorithm admits computationally tractable updates, and its stabilized objective values decrease monotonically and converge to a finite limit. Experiments on synthetic and real data show that the pro

Count-valued variables arise in many scientific and applied settings, yet explicit structural models that allow full identification of causal DAGs from observational data remain limited. The Poisson branching structural causal model (PB-SCM) provides a count-valued analogue of linear structural equation models using binomial thinning and independent Poisson exogenous variables, but its causal DAG is generally only partially identifiable. Building on this framework, we propose the Poisson thinning structural equation model (PT-SEM), which replaces binomial thinning in PB-SCM with Poisson thinning and allows node-wise exogenous distributions from diverse count-distribution families. Under node-wise regularity conditions, we establish identifiability of the causal DAG, the thinning coefficients, and the node-wise exogenous distributions. The same identification analysis extends to binomial thinning, yielding full identifiability whenever every nonsink has non-Poisson exogenous noise. We further develop a structure learning algorithm that optimizes, via dynamic programming, a BIC score based on local likelihoods evaluated at plug-in moment estimates, and establish its consistency for DAG selection. Simulations demonstrate favorable performance in DAG recovery and thinning-coefficient estimation, and a real-data application illustrates the practical utility of PT-SEM.

This paper develops fixed-smoothing (fixed-b, fixed-K) inference methods for time-series quantile regression that are robust to heteroskedasticity and autocorrelation. Our approach is uniformly valid over quantile levels and accounts for dependence both over time and across quantiles. It enables the construction of uniform confidence bands, Wald, and Sup-t tests for joint hypotheses, and tests of shape restrictions, providing a unified framework for assessing heterogeneity in quantile effects. A key challenge is that, under weak dependence, uniform inference for quantile regression processes is generally non-pivotal because the limiting distributions depend on the long-run covariance structure across quantiles. To address this issue, we develop two complementary approaches. The uniform-in-$\tau$ method estimates the covariance structure and simulates the non-pivotal limiting distribution. For certain tests involving a finite collection of quantile levels, the stack-Wald method delivers pivotal fixed-smoothing inference. We establish the asymptotic validity of both approaches. Simulation results show that the proposed methods substantially improve size control relative to existing HAC-based procedures while maintaining good power. An application to predictive quantile regressions for stock returns reveals substantial heterogeneity in predictive effects across both quantiles and forecast horizons.

Spatiotemporal areal data are commonly observed in various fields such including epidemiology, social science, economics and so on. To capture both spatial trends and temporal trends, spatiotemporal modeling is often employed, and the conditional autoregressive (CAR) model is one of the most widely used approaches for modeling areal data. This paper proposes a new framework for estimating spatiotemporal trends based on the CAR model. The proposed method provides locally adaptive temporal smoothing while yielding interpretable temporal trends by effectively utilizing information from both spatially neighboring areas and temporally adjacent time points. We also develop a Gibbs sampling algorithm and demonstrate the ability of the proposed method to adapt to to local temporal changes through numerical examples.

Principal component analysis (PCA) can rotate away from its population target when a covariance matrix is estimated from limited data. We introduce diagonal attenuation, which preserves sample cross-covariances while reducing coordinatewise sample variances. The method is revealed exactly by averaging a linear full-output reconstruction loss over random input masks; studying the correction directly extends it beyond the range attainable by masking. We isolate the part of the random coupling between retained and omitted population directions that is contributed by sample-variance errors, and show how attenuation can reduce the resulting rotation. Under balanced marginal variances, we derive an explicit expected-risk theorem, uniform over the attenuation path for all sufficiently large finite samples, and obtain the asymptotically risk-minimizing strength. For general covariances, we characterize when attenuation leaves the population PCA subspace unchanged and give a risk theorem that also accounts for changing eigengaps and the population cost when the target moves. Simulations track this tradeoff from exact preservation back to PCA. Across local image patches, speech spectra, and smartphone acceleration, both mask-derived and direct attenuation improve PCA under two fitting-sample budgets, and one of them has the largest mean gain among seven methods in every data–budget cell. The full path selects strengths beyond the mask-derived boundary on $63\%$–$95\%$ of the subsamples.

We develop a network-structured Bayesian hierarchical model for sparse association mapping between genomic alterations and quantitative treatment-response phenotypes. The framework combines a Gaussian Markov random field prior that borrows strength across pathway-connected genes, a global-local horseshoe prior inducing sparsity, and a conjugate Gibbs sampler requiring no Metropolis-Hastings steps. Though broadly applicable to high-dimensional settings with known predictor networks, we validate it using cancer cell-line drug-sensitivity data. Applied to GDSC2 ($N=951$ cell lines, $G=219$ driver genes, $D=295$ drugs), the model identifies 126 gene-drug associations (0.195% of 64{,}605 pairs), concentrated in EZH2 (45 drugs, all sensitivity-direction, mean effect $-0.911$ $\ln$IC50) and KMT2D (36 drugs, all sensitivity-direction, mean effect $-0.496$ $\ln$IC50). These markers show external support in an independent PRISM screen (1{,}518 compounds), with KMT2D achieving complete directional replication (36/36) and EZH2 partial replication (8/12). Five-fold cross-validated predictive log-likelihood confirms each prior layer's value: the full model outperforms the no-network ablation by $+3{,}109$ log-units per fold and the no-horseshoe ablation by $+14{,}039$ log-units, consistently across folds. Simulations under three scenarios show the full model achieves the highest precision and lowest false-discovery rate throughout, while the network prior improves sensitivity recovery under network-structured signal. A tissue-stratified extension identifies coherent subgroup refinements, including lung-specific EGFR-inhibitor sensitivity and skin-specific BRAF-Dabrafenib sensitivity. These results show the framework identifies sparse, interpretable, externally supported drug-sensitivity markers while enabling principled investigation of tissue-specific departures from shared effects.

We develop a Dynamic Spatial Panel Bayesian Additive Regression Trees model with Horseshoe shrinkage (DSP-BART-HS) for high-dimensional spatio-temporal panel data. We jointly evaluate the model against a comprehensive suite of structural spatial econometrics, non-parametric machine learning methods, and small-area estimators across nine data-generating scenarios spanning 200 replicates each, including irregular spatial topologies, dense policy effects, and non-linear individual-level interactions. DSP-BART-HS is the best or statistically indistinguishable from the best estimator in every scenario. Conventional region-time-aggregate comparators suffer severe performance degradation – trailing by a factor of three or more – whenever individual-level non-linearity drives outcome variance, a limitation this framework's tree-ensemble design directly addresses. The model also maintains strong predictive accuracy under a zero-training-region spatial holdout via its spatial diffusion mechanism. We demonstrate practical utility across two U.S. county-level panel applications – intergenerational economic mobility and geographic income inequality – under genuine unseen-region, random, and temporal holdouts, with results validated through repeated-split uncertainty quantification and paired significance testing against every comparator. One clear limitation emerges: under temporal extrapolation specifically, a lighter per-region autoregressive specification (MTS-CAR-X) achieves a robust, and consistent advantage. DSP-BART-HS nonetheless establishes strong predictive performance, automated variable selection, and flexible inference for hierarchical spatio-temporal panel data.

This study introduces the three-parameter Kumaraswamy–Teissier distribution (KTD) by applying the Kumaraswamy-G construction to the Teissier distribution. The proposed distribution provides additional shape flexibility for modelling positive-valued observations and accommodates different density and hazard-rate shapes. Several statistical properties, including moments, the moment generating function, quantile function, and order statistics, are derived. Maximum likelihood, maximum product spacing, and Bayesian methods are considered for parameter estimation, and simulation results demonstrate improving estimation accuracy with increasing sample size. A median-based regression model is further developed for the KTD, with its finite-sample performance examined through simulation. The practical utility of the proposed regression model is demonstrated through an application to monthly gridded rainfall data from the Northwest Himalayas, incorporating annual and semi-annual harmonic terms to capture seasonal variation. Its performance is assessed against Kumaraswamy and beta regression using residual-based goodness-of-fit diagnostics and likelihood-based model selection. The results demonstrate effectiveness of the proposed distribution for modelling heterogeneous and skewed positive-valued rainfall data.

The robustness of methods for characterizing the origin of wave-like propagating signals is often challenged by the heterogeneity of transmitting media and sensor network. Detecting an individual wave under such noisy conditions remains particularly difficult. In this paper, we propose a robust framework to estimate wave origins from spatiotemporal data collected by a sensor array. First, a Butterworth filter is designed to extract relevant signals from each sensor channel. Next, a 3D randomized Hough transform is introduced to identify candidate signals originating from the same wave front. Finally, a wave propagation model is fitted to the identified signals using the least-squares method. By combining the iterative voting mechanism of the randomized Hough transform with least-squares optimization, our proposed method demonstrates superior robustness and computational efficiency in both simulation studies and real-world data analyses.

Modern statistical models are growing increasingly complex in an effort to realistically capture system dynamics. Using standard simulation-based inference, these models may be computationally prohibitive, necessitating the use of model calibration methods. Bayesian score calibration is a computationally efficient framework for model calibration with strong theoretical guarantees. This framework learns an appropriate correction for an approximate model using a small number of simulations from the data-generating process. Currently, only a location-scale transformation has been explored, which may lack the flexibility to correct the complex error introduced by some approximate models. In this paper, we develop two flexible transformations for use in the Bayesian score calibration framework. The first is a polynomial extension, which can appropriately adjust approximate models with location-varying error. The second is a sequential application of Bayesian score calibration, which can accommodate approximate models with posteriors that have low support for the true parameter values. We also discuss an additional diagnostic for use with this framework. We demonstrate the increased flexibility these two approaches provide over Bayesian score calibration in two illustrative simulation studies.

We give a necessary and sufficient condition for the existence of power-one sequential tests in an i.i.d. composite testing problem. A level-\(\alpha\) test with power one against every alternative exists if and only if the alternatives are separated from the null by a countable family of finite-block events. We provide other equivalent conditions using randomized fixed-sample tests, bounded finite-block scores, e-processes, reduced-filtration test supermartingales, and a countable cover whose finite-block weak-$*$ closed convex hulls are positively separated in total variation. As a bonus, the constructive proof yields tests have pointwise expected sample size \(O_Q(\log(1/\alpha))\). Exactly the same conditions also characterize i.i.d.\ change detectability under optional-horizon average-run-length control: for every \(\eta>0\), they are equivalent to an alarm family \((T_\gamma)_{\gamma\ge1}\) satisfying \(\Prob_{P^\infty}(T_\gamma\le\sigma)\le \E_{P^\infty}\sigma/\gamma\) for every null law and every stopping time \(\sigma\). In fact, when these conditions hold, we can construct a single e-detector such that every null-law average run length lies between \(\gamma\) and \((1+\eta)\gamma+1\), and having robust Lorden delay \(O_Q(\log\gamma)\).

When data are grouped, hierarchical or multilevel models are commonly used to account for group-level variation with group-specific parameters. Leave-one-group-out cross-validation (LOGO-CV) is a suitable tool for evaluating predictive performance for new groups, providing an estimator of the expected log predictive density (elpd). Brute-force LOGO-CV requires one model refit per held-out group, often using computationally expensive inference algorithms such as MCMC. This is costly, particularly for large numbers of groups or complex model structures. Commonly used importance sampling approximations, intended to reduce this cost, tend to fail because the group-specific parameters of the held-out group must be integrated out. We identify two key challenges in LOGO-CV elpd estimation: approximating the LOGO posterior and computing the grouped marginal likelihood. We compare 11 strategies, including 5 newly proposed, to address them. Among others, we combine Pareto-smoothed importance sampling or adaptive importance sampling with integration techniques such as Laplace approximation, adaptive Gauss-Hermite quadrature, and bridge sampling. We evaluate these strategies in both simulation experiments and real-world case studies, which show that marginalising over the group-specific parameters substantially improves the reliability of the importance sampling approaches.

Single-cell RNA sequencing studies routinely use clustering to define putative cell types and cell states, yet the observed separation may arise from sampling variability rather than genuine biological heterogeneity. This paper studies formal significance testing of such clustering structure in high-dimensional data. Existing SigClust methods assess clustering significance through Monte Carlo simulation under a Gaussian single-cluster null, but this assumption can be unreliable for normalized gene expression data and other non-Gaussian settings. We propose SigClust-LCP, a nonparametric extension that models a single cluster by a log-concave distribution. To make this approach computationally feasible in moderate to high dimensions, we develop a score-matching estimator for log-concave projection inspired by recent generative modeling ideas. We establish theoretical guarantees for the estimator and for its use in clustering significance testing. Simulations show that SigClust-LCP controls Type-I error more reliably than existing methods across a range of unimodal and mixture distributions while retaining competitive power. In a single-cell RNA sequencing analysis of Hydra cells, the method avoids spurious subclusters within annotated cell populations and supports biologically meaningful separation across lineages and body-axis regions.

Edge selection in Gaussian graphical models is fundamentally a variable selection problem where pairwise relationships determine construct validity and variable importance in psychological networks. In psychology, network estimation relies predominantly on \(\ell_1\) regularization where uniform shrinkage systematically underestimates edge and centrality parameters. Alternative penalties overcome this bias but rely on fixed hyperparameters that do not adapt to the signal in the data. We develop a family of data-adaptive regularization penalties grounded in extreme value theory. Across 290 empirical psychological datasets, we show that absolute partial correlations are well-described by the Weibull distribution. Using this empirical regularity, we derive Weibull, Gumbel, and Exponential penalties that approximate \(\ell_0\) penalization and calibrate their hyperparameters to each dataset's noise floor. We formally prove the asymptotic properties of their static forms and conduct a large-scale simulation spanning two network topologies and various sample sizes (\(N\) = 100–10,000), demonstrating that their adaptive forms maintain high specificity while accumulating sensitivity as sample size increases with low parameter bias and high rank-order centrality congruence relative to field standards. Empirically, method choice alone determined centrality rankings at sample sizes typical in psychology. Of the three adaptive penalties, Weibull is recommended given the interpretability of its parameters.

Community detection is fundamental to understanding the modular organization in functional brain networks, yet noise in neuroimaging-derived networks and auxiliary node-level covariates pose critical challenges. Existing methods typically either assume networks are noise-free or ignore covariate information. We propose a Bayesian framework for recovering a shared latent community structure from multiple noisy network realizations and auxiliary covariates. The model combines a degree-corrected stochastic block model for the latent network, a block-structured noise model linking noisy observations to latent edges, and a covariate cluster model for node-level attributes. This specification allows anatomical or functional attributes of regions of interest to contribute information when network signals are weak or sparse. We develop an efficient Markov chain Monte Carlo algorithm for posterior sampling and select the number of communities using the widely applicable information criterion, avoiding prior specification of this quantity. Simulation studies demonstrate improved community recovery relative to existing methods across varying noise levels, covariate signal strengths, and numbers of noisy networks, with larger gains when network noise is moderate to high or only a small number of noisy networks is available. Applications to functional brain networks from the Alzheimer's Disease Neuroimaging Initiative and the Human Connectome Project identify biologically interpretable structures and capture disease-related reorganization and individual-level variation.

Estimating a low-dimensional subspace shared across noisy data matrices is a fundamental problem in multi-view matrix estimation. We study this problem under the two-view JIVE model, where each data matrix contains shared and view-specific low-rank components. We demonstrate that standard plug-in subspace intersection, including AJIVE, suffers from a second-order bias caused by direction-dependent leakage of the empirical singular vectors. We propose bias-corrected subspace intersection (BCSI), which removes this bias before estimating the shared subspace. We establish finite-sample risk bounds for BCSI that accommodate unequal view dimensions, signal strengths, and view-specific ranks and require no condition-number assumptions on the signal matrices. When the shared and view-specific ranks are comparable, these bounds match our minimax lower bounds up to universal constants. The resulting minimax rate contains a new second-order term, arising from quadratic leakage perturbations relative to the shrinking spectral gap when the view-specific subspaces are nearly aligned. This term is absent from previous JIVE minimax lower bounds. Numerical experiments demonstrate the advantage of BCSI over AJIVE when the leakage bias is pronounced. Along the way, we establish a nonasymptotic concentration result for the bias-corrected leakage Gram matrix of a rectangular spiked matrix, which may be of independent interest.

We study finite-sample valid tests of the one-sided mean hypothesis $H_0:\mu\leq 1$ against $H_1:\mu>1$ for nonnegative random variables. To do so, we develop a leave-one-out dual certificate framework, where certain pointwise inequalities imply p-value validity under the conditional mean null $\mathbb{E}[X_i\mid\mathbf{X}_{-i}]\leq 1$, and which also gives conditions that allow combining dual certificates for p-values to show that their pointwise minimum is also a valid p-value. The framework proves finite-sample validity of Wang and Zhao's nonparametric likelihood-ratio statistic $T_{\mathrm{nplr}}$, yields a new p-value $T_{\mathrm{bin}+}$ extending the Clopper–Pearson binomial test to general nonnegative random variables, and shows that the pointwise minimum $\min\{T_{\mathrm{nplr}},T_{\mathrm{bin}+}\}$ is itself a valid and more powerful p-value. We establish sharp optimality results for such testing problems in two regimes: both $T_{\mathrm{nplr}}$ and $T_{\mathrm{bin}+}$ attain a universal detectability boundary for the null $H_0$ without moment or tail assumptions, and $T_{\mathrm{bin}+}$ attains a nonparametric power lower bound under $n^{-1/2}$-local alternatives to $H_0$. Efficient algorithms and numerical experiments demonstrate substantial finite-sample power gains over existing valid methods.

We develop a finite-horizon calibration method for raw-score martingale posteriors, with von Mises–Fisher models as the main worked example. Starting from the maximum likelihood estimator, predictive paths are generated by simulating future observations from the current fitted model and updating the natural parameter by unpreconditioned score increments. The main methodological step is to separate predictive simulation from covariance calibration. Raw-score increments have Fisher-information covariance, whereas Bernstein–von Mises calibration requires inverse-information covariance. We therefore apply a terminal linear correction based on a local information estimate. For more efficient implementation, we also introduce a hybrid version that replaces the omitted tail of the infinite predictive continuation by a Gaussian approximation with matching leading-order quadratic variation. We prove fixed-n convergence, a finite-horizon approximation bound for the Gaussian tail, and a Bernstein–von Mises limit for the hybrid post-corrected sampler under local regularity and consistent terminal calibration. Simulations show that tail correction reduces truncation-induced underdispersion, and an OSCAR ocean-current example illustrates local directional uncertainty summaries.

We present an approach for determining sample sizes required to detect underestimations of the expected shortfall with a prescribed power when applying the recently proposed e-backtesting procedure. We consider scenarios in which the value-at-risk at level $p$ is always estimated correctly, while the difference between the true expected shortfall and the value-at-risk is underestimated by a given factor $r$. We show that exploiting the structure of the backtest e-statistic proposed for backtesting the expected shortfall at level $p$ enables the derivation of approximate lower bounds for the required sample sizes by considering a sequence of independent and identically distributed Bernoulli random variables. We also discuss potential limitations of this approximation and compare the resulting sample size requirements with those obtained in practical applications using Monte Carlo simulations. Furthermore, we present generalizations of the e-backtesting procedure, in particular to risk measures which constitute Bayes pairs.

Problem definition: Fixed parameters and inputs in a stochastic simulator induce a distribution over complete trajectories, not one trajectory. Calibration must assess this distribution, including variability and temporal dependence, against observations. Yet stochastic car-following models are commonly calibrated with trajectory-error objectives inherited from deterministic modelling. Methodology/results: We establish a scoring-rule theory of stochastic calibration. Strict propriety requires the data-generating distribution to uniquely minimise expected score. MRMean-I, the average run-wise error, drives separable stochastic spread to zero; MRMean-II, the error of the ensemble-mean trajectory, cannot identify a parameter that changes only spread; and MRMin, the error of the closest simulated run, has a population target that changes with ensemble size. These results are confirmed for stochastic Intelligent Driver Model extensions with additive acceleration noise and random desired headway. We recommend exact maximum likelihood when the correct transition density is available; otherwise, an unbiased simulation-based estimator of a strictly proper score. The energy score meets this requirement and gives the best held-out distributional prediction among the evaluated simulation-based objectives, although both models retain too-narrow bands and miss persistent disturbances. Implications:Strict propriety separates a valid calibration target from parameter identifiability and model adequacy. The theory applies to vector-valued outputs from stochastic transportation simulators; the car-following experiments illustrate its scope.

Multivariate correlated outcomes occur across disciplines, including ecology, social sciences, and psychometrics. This paper focuses on clustering these outcomes across observational units, specifically, finding groups of units with the same “outcome profile". Our motivation comes from bioregionalization in ecology, which aims to cluster sites into bioregions with the same species profiles, where site membership can depend on environmental or habitat covariates. To accomplish this, we propose finite mixtures of generalized estimating equations (MixGEE). Unlike existing approaches to model-based bioregionalization, MixGEE partitions sites into regions while accounting for between-species correlations through a region-specific working correlation structure. Thus, each region is characterized by a marginal species mean vector and a between-species correlation matrix. Unlike likelihood-based finite mixture models, MixGEE does not require a full joint distribution of the multivariate outcomes. Instead, we construct a pseudo-posterior probability for region membership motivated by the large-sample distribution of the estimating equation. This leads to an iterative algorithm alternating between updating these probabilities and solving weighted estimating equations. We determine the number of regions using cross-validation based on predictive performance for held-out sites and species components, and use a clustered Dirichlet random-weight bootstrap for uncertainty quantification. Simulations demonstrate reliable estimation and inference under various correlation structures and more stable selection of the number of groups than methods that ignore dependence. Applying MixGEE to presence–absence records of fish species around the Kerguelen Plateau reveals three distinct fish assemblage profiles with heterogeneous occurrence and within-site correlation patterns.

Longitudinal clinical data are increasingly available, offering opportunities to study treatment effects over time but also raising challenges related to time-varying confounding and evolving treatment decisions. We investigate how longitudinal treatment dynamics affect causal effect estimation when baseline and longitudinal methods target different causal estimands. Using a structural causal model (SCM), we simulate data with time-varying confounding, binary treatment, and an absorbing binary outcome under 9 scenarios combining functional complexity and treatment persistence. We compare three baseline and two longitudinal estimators against Monte Carlo ground-truth risk ratios (RRs) under sustained treatment regimes. Results show that treatment persistence is the main driver of estimator behaviour. High persistence reduces the discrepancy between baseline and sustained-regime estimands, making baseline estimators closer to the sustained-regime ground truth (average relative deviation of baseline IPTW/TMLE decreasing from 61% under low persistence to 10% under high persistence), whereas low persistence induces practical positivity challenges and increases the variability of longitudinal estimators, with empirical 95% interval widths increasing from 0.22 to 0.46 for longitudinal IPTW and from 0.20 to 0.36 for LTMLE when moving from high to low persistence. These findings emphasise that estimator performance should be interpreted jointly with the target intervention and the treatment process generating the observed data.

The statistical analysis of heavy-tailed data has received considerable attention because extreme observations frequently arise in many practical applications. The Pareto type-I distribution is a fundamental heavy-tailed model used in economics, finance, actuarial science, insurance, reliability, and extreme value analysis. In this paper, we propose novel goodness-of-fit tests for the Pareto distribution using a mean residual life characterization. The test statistic is constructed using U-statistic theory, and its asymptotic behaviour is established under both the null and alternative hypotheses. Its finite-sample performance is evaluated through Monte Carlo simulations using maximum-likelihood and method-of-moments estimation and compared with existing tests. The results show that the proposed test controls the nominal significance level and performs competitively in terms of power across a broad range of alternatives. Finally, the proposed methodology is illustrated using the Danish fire insurance loss and pollution datasets.

We propose a Bayesian model-based approach for clustering functional data with an unknown number of clusters and within-curve correlated observations. Cluster-specific mean functions are represented using B-spline basis expansions, while within-curve dependence is modeled through an Ornstein–Uhlenbeck covariance structure. A truncated Dirichlet process mixture is used to infer the effective number of clusters, and a variational EM algorithm is developed for efficient posterior approximation. Simulation studies show that the proposed method performs well under both correctly specified and misspecified mean-function settings and achieves higher average values of the reported clustering metrics than the competing methods in the simulation settings considered. Comparisons with MCMC indicate that the variational approximation produces clustering results and parameter estimates that are in close agreement with those obtained by MCMC, while requiring substantially lower computational cost. An application to Canadian daily temperature curves further demonstrates the practical usefulness of the method in identifying interpretable functional clusters while accounting for within-curve dependence.

This study proposes two novel bivariate distributions for jointly modeling temperature and rainfall by integrating Kumaraswamy-Teissier marginals with Clayton and Gumbel copula structures. To capture a wide range of dependence patterns, including both positive and negative associations, rotated copula variants (90{\deg}, 180{\deg}, and 270{\deg}) are incorporated along with their corresponding tail dependence characteristics. Model parameters are estimated using maximum likelihood and the inference functions for margins (IFM) approach, and their finite-sample performance is assessed through a comprehensive Monte Carlo simulation study. The proposed models are applied to monthly gridded temperature and rainfall data from the Northwest Himalayas, a region characterized by complex hydro-climatic variability. Comparative analysis demonstrates that the proposed framework outperforms several existing bivariate models and effectively captures lower-tail, upper-tail, and asymmetric dependence structures across summer and winter seasons. Based on the selected best-fitting copula models, univariate, joint, and conditional return periods are derived to quantify the risk of compound extremes. The results highlight the capability of the proposed approach to provide a more realistic representation of hydro-climatic dependence and offer a robust framework for assessing the risk of extreme temperature and rainfall events in mountainous regions.

This paper introduces a quantile-based Kumaraswamy-Teissier autoregressive moving average (KTARMA) model for positive-valued time series. Leveraging the flexibility of the extended Kumaraswamy-Teissier distribution within an observation-driven framework, the random component of the distribution is conditioned on the historical process and time-varying covariates, and is parameterized explicitly via its $\rho$-th conditional quantile, where $\rho \in (0,1)$. To capture temporal dependence, the systematic component maps an ARMA-type structure to this conditional quantile via an appropriate link function. For inference, we implement a conditional maximum likelihood framework and derive explicit analytical expressions for the resulting score vector and conditional information matrix, followed by the development of model diagnostic and forecasting procedures. The finite-sample performance of the developed estimators is evaluated through a Monte Carlo simulation study across various parameter configurations and quantile levels. Finally, the practical utility of the study is demonstrated by modeling monthly rainfall data over the Northwest Himalayas (2001-2025), where 525 grids are grouped into four homogeneous zones using a Self-Organizing Map and relevant atmospheric variables and large-scale climate indices are incorporated as predictive regressors. Out-of-sample forecasting evaluations reveal that the KTARMA model delivers highly competitive predictive performance, achieving consistently lower mean squared errors across all identified zones compared to KARMA and $\beta$ARMA models.

In medical research, it is often of interest to evaluate the predictive performance of a biomarker. Statistical approaches based on the Receiver Operating Characteristic (ROC) curve and its summary measures, such as the area under the curve (AUC) and the Youden index, are widely used to evaluate the prognostic performance of these biomarkers. In time-to-event studies, ROC analysis poses additional challenges due to change in disease status over time and the presence of censored individuals. To address these issues, time-dependent ROC curves were introduced. In this paper, we propose a non-parametric estimator of the time-dependent two-way partial AUC for right-censored data. We also discuss the partial Youden index and the associated optimal biomarker cutoff estimator for the right-censored data. We conduct an extensive simulation study to investigate the finite sample performance of the proposed estimators. The simulation study indicates that the proposed non-parametric estimators efficiently account for right censoring. Finally, we illustrate the proposed methods using two real data sets, one from the Primary Biliary Cirrhosis study and the other from the Molecular Taxonomy of Breast Cancer International Consortium trial.

In this paper we introduce the R package mccca, which implements multiple-class cluster correspondence analysis (MCCCA) proposed in M.Takagishi et al., (2022). MCCCA is a statistical method that identifies and visualizes heterogeneous tendencies specific to “classes” (e.g., gender and nationality) in a low dimensional space. In MCCCA, two kinds of variables, external and active variables, are distinguished. External variables directly define classes whereas the active variables are used to derive class-specific clusters of observations (individuals). In this paper, we show how to apply MCCCA and how to visualize its results using mccca. We illustrate the procedure for performing MCCCA by applying it to two data sets.