Papers

Seasonality has traditionally shaped the U.S. housing market, with activity peaking in spring-summer and declining in autumn-winter. However, recent disruptions, particularly post-COVID-19, raise questions about shift in these patterns. This study analyzes housing market date (1991-2024) to examine evolving seasonality and regional heterogeneity. Using Housing Price Index (HPI), inventory and sales data from the Federal Housing Finance Agency and U.S. Census Bureau, seasonal components are extracted via the X-13-ARIMA procedure, and statistical tests assess variations across regions. The results confirm seasonal fluctuations in prices and volumes, with recent shifts toward earlier annual peak (March-April) and amplified seasonal effects. Regional variations align with differences in climate and market structure, while prices and sales volumes exhibit in-phase movement, suggesting thick-market momentum behaviour. These findings highlight key implications for policymakers, realtors and investors navigating post-pandemic market dynamics, offering insights into the timing and interpretation of housing market activities.

Counterfactual decision-making in the face of uncertainty involves selecting the optimal action from several alternatives using causal reasoning. Decision-makers often rank expected potential outcomes (or their corresponding utility and desirability) to compare the preferences of candidate actions. In this paper, we study new counterfactual decision-making rules by introducing two new metrics: the probabilities of potential outcome ranking (PoR) and the probability of achieving the best potential outcome (PoB). PoR reveals the most probable ranking of potential outcomes for an individual, and PoB indicates the action most likely to yield the top-ranked outcome for an individual. We then establish identification theorems and derive bounds for these metrics, and present estimation methods. Finally, we perform numerical experiments to illustrate the finite-sample properties of the estimators and demonstrate their application to a real-world dataset.

Accurate mortality modeling is central to actuarial science and public health, especially as mental health emerges as a significant factor in population outcomes. This paper develops and applies a Bayesian hierarchical model to analyze U.S. county-level mortality and suicide rates from 2010 to 2023. Applying a conditional autoregressive (CAR) structure to each combination of sex and age grouping, the model captures spatial and temporal trends while incorporating mental health surveillance data and socio-economic indicators. We first assess socio-economic covariates in predicting suicide. While the results vary considerably by age and sex, we find that the county-wide levels of educational attainment, housing prices, marriage rates, racial composition, household size, and poor mental health days all have significant relationships with suicide rates. We next consider the impact of various mental health indicators on all-cause and suicide-specific mortality and find that the strongest effects are observed in younger populations. The spatial and temporal correlation structures reveal substantial regional clustering and time-consistent trends in both all-cause mortality and suicide rates, supporting the use of spatio-temporal methods. Our findings highlight the value of integrating mental health surveillance data into mortality models to better identify emerging risk areas and vulnerable populations. This approach has the potential to inform public health policy, resource allocation, and targeted interventions aimed at reducing disparities in mortality and suicide across U.S. communities.

We study a sequential price competition among $N$ sellers, each influenced by the pricing decisions of their rivals. Specifically, the demand function for each seller $i$ follows the single index model $\lambda_i(\mathbf p) = \mu_i(\langle \boldsymbol \theta_{i,0}, \mathbf p \rangle)$, with known increasing link $\mu_i$ and unknown parameter $\boldsymbol \theta_{i,0}$, where the vector $\mathbf{p}$ denotes the vector of prices offered by all the sellers simultaneously at a given instant. Each seller observes only their own realized demand - unobservable to competitors - and the prices set by rivals. We propose a novel decentralized policy, PML-GLUCB, that combines penalized MLE with an upper-confidence pricing rule. Our approach (i) removes the need for coordinated front-loaded exploration phases across sellers - which is integral to previous models - making our method aligned with realistic market conditions; (ii) generalizes existing approaches that focus solely on linear demand models; (iii) accommodates both binary and real-valued demand observations. Relative to a dynamic benchmark policy, each seller achieves $\widetilde{O}(\sqrt{T})$ regret, which matches the optimal rate known in the linear setting.

A well-defined distance on the parameter space is key to evaluating estimators, ensuring consistency, and building confidence sets. While there are typically standard distances to adopt in a continuous space, this is not the case for combinatorial parameters such as graphs that represent statistical models. Defined on the graphs alone, existing proposals like the structural Hamming distance ignore the structure of the model space and can thus exhibit undesirable behaviors. We propose a model-oriented framework for defining the distance between graphs that is applicable across different graph classes. Our approach treats each graph as a statistical model and organizes the graphs in a partially ordered set based on model inclusion. This induces a neighborhood structure, from which we define the model-oriented distance as the length of a shortest path through neighbors, yielding a metric in the space of graphs. We apply this framework to probabilistic undirected graphs, causal directed acyclic graphs, {causal acyclic directed mixed graphs}, probabilistic completed partially directed acyclic graphs, and causal maximally oriented partially directed acyclic graphs. We analyze theoretical and empirical behaviors of the model-oriented distance. By exploiting the underlying poset structures, we develop algorithms for computing and bounding the proposed distance that scale to moderate-sized graphs. Finally, we showcase its utility for quantifying the robustness of adjustment sets to errors in specifying the causal graph.

The improving multi-armed bandits problem is a formal model for allocating effort under uncertainty, motivated by scenarios such as investing research effort into new technologies, performing clinical trials, and hyperparameter selection from learning curves. Each pull of an arm provides reward that increases monotonically with diminishing returns. A growing line of work has designed algorithms for improving bandits, albeit with somewhat pessimistic worst-case guarantees. Indeed, strong lower bounds of $\Omega(k)$ and $\Omega(\sqrt{k})$ multiplicative approximation factors are known for both deterministic and randomized algorithms (respectively) relative to the optimal arm, where $k$ is the number of bandit arms. In this work, we propose two new parameterized families of bandit algorithms and bound the sample complexity of learning the near-optimal algorithm from each family using offline data. We also perform empirical evaluations on standard hyperparameter tuning benchmarks. The first family we define includes the optimal randomized algorithm from prior work. We show that an appropriately chosen algorithm from this family can achieve stronger guarantees, with optimal dependence on $k$, when the arm reward curves satisfy additional properties related to the strength of concavity. Our second family contains algorithms that both guarantee best-arm identification on well-behaved instances and revert to worst-case guarantees on poorly-behaved instances.

Using 2.6 billion geolocated social-media posts (2014-2022) and a fine-tuned generative language model, we construct county-level indicators of life satisfaction and happiness for the United States. We document an apparent rural-urban paradox: rural counties express higher life satisfaction while urban counties exhibit greater happiness. We reconcile this by treating the two as distinct layers of subjective well-being, evaluative vs. hedonic, showing that each maps differently onto place, politics, and time. Republican-leaning areas appear more satisfied in evaluative terms, but partisan gaps in happiness largely flatten outside major metros, indicating context-dependent political effects. Temporal shocks dominate the hedonic layer: happiness falls sharply during 2020-2022, whereas life satisfaction moves more modestly. These patterns are robust across logistic and OLS specifications and align with well-being theory. Interpreted as associations for the population of social-media posts, the results show that large-scale, language-based indicators can resolve conflicting findings about the rural-urban divide by distinguishing the type of well-being expressed, offering a transparent, reproducible complement to traditional surveys.

Neural Differential Equations (NDEs) excel at modeling continuous-time dynamics, effectively handling challenges such as irregular observations, missing values, and noise. Despite their advantages, NDEs face a fundamental challenge in adopting dropout, a cornerstone of deep learning regularization, making them susceptible to overfitting. To address this research gap, we introduce Continuum Dropout, a universally applicable regularization technique for NDEs built upon the theory of alternating renewal processes. Continuum Dropout formulates the on-off mechanism of dropout as a stochastic process that alternates between active (evolution) and inactive (paused) states in continuous time. This provides a principled approach to prevent overfitting and enhance the generalization capabilities of NDEs. Moreover, Continuum Dropout offers a structured framework to quantify predictive uncertainty via Monte Carlo sampling at test time. Through extensive experiments, we demonstrate that Continuum Dropout outperforms existing regularization methods for NDEs, achieving superior performance on various time series and image classification tasks. It also yields better-calibrated and more trustworthy probability estimates, highlighting its effectiveness for uncertainty-aware modeling.

We study the long-standing problem of determining the number of principal components in econometric applications from a selective inference perspective. We consider i.i.d. observations from a $p$-dimensional random vector with $p<n$ and define the “true” dimensionality as the rank of the population covariance matrix. Building on the sequential testing viewpoint, we propose a data-driven procedure that estimates $\rank(\Sigma_X)$ using a statistic that depends on the eigenvalues of the sample covariance matrix. While the test statistic shares the functional form of its fixed design counterpart Choi et al. (2017), our analysis departs from the non-stochastic setting by treating the design as random and by avoiding parametric Gaussian assumptions. Under a locally defined null hypothesis, we establish asymptotically exact type~I error controls in the sequential testing procedure, with simulation results indicating empirical validity of the proposed method.

In this work we study the diffusion annealed Langevin dynamics, a score-based diffusion process recently introduced in the theory of generative models and which is an alternative to the classical overdamped Langevin diffusion. Our goal is to provide a rigorous construction and to study the theoretical efficiency of these models for general base distribution as well as target distribution. As a matter of fact these diffusion processes are a particular case of Nelson processes i.e. diffusion processes with a given flow of time marginals. Providing existence and uniqueness of the solution to the annealed Langevin diffusion leads to proving a Poincaré inequality for the conditional distribution of $X$ knowing $X+Z=y$ uniformly in $y$, as recently observed by one of us and her coauthors. Part of this work is thus devoted to the study of such Poincaré inequalities. Additionally we show that strengthening the Poincaré inequality into a logarithmic Sobolev inequality improves the efficiency of the model.

Continuous-time stochastic processes underlie many natural and engineered systems. In healthcare, autonomous driving, and industrial control, direct interaction with the environment is often unsafe or impractical, motivating offline reinforcement learning from historical data. However, there is limited statistical understanding of the approximation errors inherent in learning policies from offline datasets. We address this by linking reinforcement learning to the Hamilton-Jacobi-Bellman equation and proposing an operator-theoretic algorithm based on a simple dynamic programming recursion. Specifically, we represent our world model in terms of the infinitesimal generator of controlled diffusion processes learned in a reproducing kernel Hilbert space. By integrating statistical learning methods and operator theory, we establish global convergence of the value function under sufficiently small generator-approximation error and derive explicit finite-sample guarantees with bounds tied directly to system properties such as smoothness and stability. Our theoretical and proof-of-concept results indicate that operator-based approaches may hold promise for solving offline reinforcement learning using continuous-time optimal control.

This paper presents an experimental evaluation of parallel-in-time Kalman filters and smoothers using graphics processing units (GPUs). In particular, the paper evaluates different all-prefix-sum algorithms, that is, parallel scan algorithms for temporal parallelization of Kalman filters and smoothers in two ways: by calculating the required number of operations via simulation, and by measuring the actual run time of the algorithms on real GPU hardware. In addition, a novel parallel-in-time two-filter smoother is proposed and experimentally evaluated. Julia code for Metal and CUDA implementations of all the algorithms is made publicly available.

In this paper, we propose a new flexible family of distributions for data that consist of three angles, two angles and one linear component, or one angle and two linear components. To achieve this, we equip the recently proposed trivariate wrapped Cauchy copula with non-uniform marginals and develop a parameter estimation procedure. We compare our model to its main competitors for analyzing trivariate data and provide some evidence of its advantages. We illustrate our new model using toroidal data from protein bioinformatics of conformational angles, and cylindrical data from climate science related to buoy in the Adriatic Sea. The paper is motivated by these real trivariate datasets.

We present methods that provide all zeroes and extrema of a function that do not require differentiation. Using point process theory, we are able to describe the locations of zeroes or maxima, their number, as well as their distribution over a given window of observation. The algorithms in order to accomplish the theoretical development are also provided, and they are exemplified using many illustrative examples, for real and complex functions.

This paper introduces a heuristic framework for the Best Secretary Problem, where one item must be selected using rank information only. We develop five data-responsive rules extending classical fixed-cutoff methods: an expected-record threshold, an adaptive deviation correction, a probabilistic early-accept rule, a two-phase relaxation, and a local dynamic programming approximation. These rules adjust thresholds sequentially as information accumulates. Simulations across diverse sample sizes, distributions, and autocorrelated settings show that the heuristics match or exceed traditional optimal rules in stability and efficiency. The expected-record rule remains strong despite its simplicity, the adaptive correction performs well under asymmetry, and the adaptive and probabilistic rules reduce average stopping times. An ensemble combining multiple rules yields the most stable performance. Overall, a few intuitive parameters achieve near-optimal results, demonstrating that data-responsive heuristics can effectively extend rank-based optimal stopping to dynamic decision environments.

The challenge of distributed fusion estimation is investigated for a class of four-dimensional (4D) commutative hypercomplex signals that are $\mathbb{T}_k$-proper factorizable, within the framework of multiple-sensor networks with different fading measurement rates. The fading effects affecting each sensor's measurements are modeled as a stochastic variables with known second-order statistical properties. The estimation process is conducted exclusively based on these second-order statistics. Then, by exploiting the $\mathbb{T}_k$-properness property within a tessarine framework, the dimensionality of the problem is significantly reduced. This reduction in dimensionality enables the development of distributed fusion filtering, prediction, and smoothing algorithms that entail lower computational effort compared with real-valued approaches. The performance of the suggested algorithms is assessed through numerical experiments under various uncertainty conditions and $T_k$-proper contexts. Furthermore, simulation results confirm that $\mathbb{T}_k$-proper estimators outperform their quaternion-domain counterparts, underscoring their practical advantages. These findings highlight the potential of $\mathbb{T}_k$-proper estimation techniques for improving multi-sensor data fusion in applications where efficient signal processing is essential.

Both Hawkes processes and autoregressive processes rely on linear functionals of their past, while modeling different types of data. Since datasets arising from observations of the same phenomenon may be heterogeneous and sampled at different time scales, it is natural to study multiscale and heterogeneous processes, such as those obtained by combining Hawkes and autoregressive dynamics. In this paper, we introduce this new Hawkes autoregressive (HAR) model incorporating both continuous- and discrete-time dynamics, and establish several probabilistic results, including the existence of a stationary version, a cluster representation, as well as stability and ergodic properties.

The fractional Brownian motion (fBm) is parameterized by the Hurst exponent $H\in(0,1)$, which determines the dependence structure and regularity of sample paths. Empirical findings suggest that the Hurst exponent may be non-constant in time, giving rise to the so-called multifractional Brownian motion (mBm). The Itô-mBm is an alternative to the classical mBm, and has been shown to admit more intuitive sample path properties in case the Hurst function is rough. In this paper, we show that the Itô-mBm also allows for a simplified statistical treatment compared to the classical mBm. In particular, estimation of the local Hurst parameter $H(t)$ with Hölder exponent $\eta>0$ achieves rates of convergence which are standard in nonparametric regression, whereas similar results for the classical mBm only hold for the smoother regime $\eta>1$. Furthermore, we derive an estimator of the integrated Hurst exponent $\int_0^t H(s)\, ds$ which achieves a parametric rate of convergence, and use it to construct goodness-of-fit tests.

In change-point analysis, one aims at finding the locations of abrupt distributional changes (if any) in a sequence of multivariate observations. In this article, we propose some nonparametric methods based on averages of pairwise distances for this purpose. These distance-based methods can be conveniently used for high-dimensional data even when the dimension is much larger than the sample size (i.e., the length of the sequence). We carry out some theoretical investigations on the behaviour of these methods not only when the dimension of the data remains fixed and the sample size grows to infinity, but also in situations where the dimension diverges to infinity while the sample size may or may not grow with the dimension. Several high-dimensional datasets are analyzed to compare the empirical performance of these proposed methods against some state-of-the-art methods.

Violations of the positivity assumption can render conventional causal estimands unidentifiable, including the average treatment effect (ATE), the average treatment effect on the treated (ATT), and the average treatment effect on the controls (ATC). Shifting the inferential focus to their alternative counterparts – the weighted ATE (WATE), the weighted ATT (WATT), and the weighted ATC (WATC) – offers valuable insights into treatment effects while preserving internal validity. In this tutorial, we provide a comprehensive review of recent advances in propensity score (PS) weighting methods, along with practical guidance on how to select a primary target estimand (while other estimands serve as supplementary analyses), implement the corresponding PS-weighted estimators, and conduct post-weighting diagnostic assessments. The tutorial is accompanied by a user-friendly R package, ChiPS. We demonstrate the pertinence of various estimators through extensive simulation studies. We illustrate the flow of the tutorial on two real-world case studies: (i) Effect of smoking on blood lead level using data from the 2007-2008 National Health and Nutrition Examination Survey (NHANES); and (ii) Impact of history of sex work on HIV status among transgender women in South Africa.

The masking-one-out (MOO) procedure, masking an observed entry and comparing it versus its imputed values, is a very common procedure for comparing imputation models. We study the optimum of this procedure and generalize it to a missing data assumption and establish the corresponding semi-parametric efficiency theory. However, MOO is a measure of prediction accuracy, which is not ideal for evaluating an imputation model. To address this issue, we introduce three modified MOO criteria, based on rank transformation, energy distance, and likelihood principle, that allow us to select an imputation model that properly account for the stochastic nature of data. The likelihood approach further enables an elegant framework of learning an imputation model from the data and we derive its statistical and computational learning theories as well as consistency of BIC model selection. We also show how MOO is related to the missing-at-random assumption. Finally, we introduce the prediction-imputation diagram, a two-dimensional diagram visually comparing both the prediction and imputation utilities for various imputation models.

Continuous proportions measured on the same experimental unit often pose two challenges: interior outliers that inflate variance beyond the beta ceiling and residual dependence that invalidates independent-margin models. We introduce a Bayesian copula modeling approach that combines rectangular-beta margins, which temper interior outliers by reallocating mass from the peak to a uniform component, with a single-parameter copula to capture concordance. Gaussian, Gumbel, and Clayton copula families are fitted, and log marginal likelihoods are obtained via bridge sampling to guide model selection. Applied to a 13-year survey (2003-2016) of Azorella selago cushion plants on sub-Antarctic Marion Island, the copula models outperform independence baselines in explaining percent dead stem cover. Accounting for between-year dependence uncovers a positive west-slope effect and weakens the cushion size effect. Simulation results show negligible bias and near-nominal 95% highest posterior density coverage across a range of tail weight and dependence scenarios, confirming good frequentist properties. The method integrates readily with JAGS and provides a robust default for paired proportion data in ecology and other disciplines where bounded outcomes and occasional outliers coincide.

We address the general task of learning with a set of candidate models that is too large to have a uniform convergence of empirical estimates to true losses. While the common approach to such challenges is SRM (or regularization) based learning algorithms, we propose a novel learning paradigm that relies on stronger incorporation of empirical data and requires less algorithmic decisions to be based on prior assumptions. We analyze the generalization capabilities of our approach and demonstrate its merits in several common learning assumptions, including similarity of close points, clustering of the domain into highly label-homogeneous regions, Lipschitzness assumptions of the labeling rule, and contrastive learning assumptions. Our approach allows utilizing such assumptions without the need to know their true parameters a priori.

Wearable devices are often used in clinical and epidemiological studies to monitor physical activity behavior and its influence on health outcomes. These devices are worn over multiple days to record activity patterns, such as step counts recorded at the minute level, resulting in multi-level, longitudinal, high-dimensional, or functional data. When monitoring patterns of step counts over multiple days, devices may record excess zeros during periods of sedentary behavior or non-wear times. Additionally, it has been demonstrated that the accuracy of wearable devices in monitoring true physical activity patterns depends on the intensity of the activities and wear times. While work on adjusting for biases due to measurement errors in functional data is a growing field, relatively less work has been done to study the occurrence of excess zeros along with measurement errors and their combined influence on estimation and inference in multi-level scalar-on-function regression models. We propose semi-continuous modeling approaches to adjust for biases due to zero inflation and measurement errors in scalar-on-function regression models. We provide theoretical justifications for our proposed methods and, through extensive simulations, we demonstrated their finite sample properties. Finally, the developed methods are applied to a school-based intervention study examining the association between school day physical activity with age- and sex-adjusted body mass index among elementary school-aged children.

In Wasserstein geometry, one-dimensional location-scale models are flat both intrinsically and extrinsically-that is, they are curvature-free as well as totally geodesic in the space of probability distributions. In this study, we introduce a class of one-dimensional statistical models, termed the location-scale-shape model, which generalizes several distributions used in extreme-value theory. This model has a shape parameter that specifies the tail heaviness. We investigate the Wasserstein geometry of the location-scale-shape model and show that it is intrinsically flat but extrinsically curved.

Road traffic in developing countries, such as India, features a heterogeneous mix of vehicles operating under weak lane discipline (HWLD), encompassing both motorised and non-motorised modes with diverse sizes and manoeuvrability. These conditions lead to complex driver interactions, complicating the reliable identification of vehicle-following (VF) behaviour and leader-follower (LF) pairs. Traditional identification methods based on fixed thresholds for longitudinal and lateral proximity often misclassify non-following instances as valid LF pairs, degrading model performance. This study presents a refined and adaptive method for LF identification in HWLD traffic. It employs vehicle-type- and speed-specific desirable gap thresholds derived from the fundamental density-speed diagram to eliminate false-positive pairs. Additionally, Mexican Hat Wavelet Transform (MWT) is employed to analyse LV and SV speed profiles, verifying LV-SV interaction for LF pair identification. The three-stage filtering includes: (i) speed-gap consistency, (ii) approach/diverge detection via relative velocity sign changes and gap range, and (iii) wavelet-based speed correlation using MWT to confirm LV influence on SV. The framework effectively filters out LF pairs associated with overtaking, tailgating, and inconsistent gap dynamics, retaining only those with consistent VF behaviour and improving model accuracy. Analysis across thirteen LF combinations shows that VF dynamics depend on both SV and LV types. Symmetric pairs (e.g., CAR-CAR, AUTO-CAR) exhibit higher predictability and lower errors, while asymmetric pairs with heavy vehicles or two-wheelers show greater variability. The framework offers a robust foundation for traffic modelling and behaviour analysis.

Gradient descent dynamics on the deep matrix factorization problem is extensively studied as a simplified theoretical model for deep neural networks. Although the convergence theory for two-layer matrix factorization is well-established, no global convergence guarantee for general deep matrix factorization under random initialization has been established to date. To address this gap, we provide a polynomial-time global convergence guarantee for randomly initialized gradient descent on four-layer matrix factorization, given certain conditions on the target matrix and a standard balanced regularization term. Our analysis employs new techniques to show saddle-avoidance properties of gradient decent dynamics, and extends previous theories to characterize the change in eigenvalues of layer weights.

We show that using autonomous-flow-based generation, one can universally approximate orientation-preserving diffeomorphisms defined on the cube by Neural ODEs with rate $\mathcal{O}(P^{-1/d})$ with $P$ parameters. On the other hand, we show that by using only a single autonomous flow, the class of Neural ODEs is nowhere dense on the cube in dimension $d \ge 2$ . Under a compact-support$_\mathrm{id}$ condition on $(0,1)^d$, we show that using autonomous-flow-based generation, one can universally approximate compactly supported$_\mathrm{id}$ diffeomorphisms on $(0,1)^d$ for any dimension with rate $\mathcal{O}((\frac{P}{\log P})^{-2/d})$ with $P$ parameters and for compactly supported$_\mathrm{id}$ homeomorphisms on $(0,1)^d$ in dimension $d \geq 5$ with rate $\mathcal{O}(P^{-1/(d+1)})$ with $P$ parameters and by a composition of at most $I_d$ autonomous Neural ODEs with the same support$_\mathrm{id}$, where $I_d$ depends only on the dimension. Moreover, we show that the class of single autonomous flows compactly supported$_\mathrm{id}$ on $(0,1)^d$ is meagre in the space of compactly supported$_\mathrm{id}$ homeomorphisms on $(0,1)^d$ for $d\ge 2$. By linearly lifting the domain into one higher dimension, we obtain a universal approximation result for Lipschitz functions compactly supported on $(0,1)^d$ with rate $\mathcal{O}(P^{-1/(d+1)})$ with $P$ parameters.

Structured variational inference constitutes a core methodology in modern statistical applications. Unlike mean-field variational inference, the approximate posterior is assumed to have interdependent structure. We consider the natural setting of star-structured variational inference, where a root variable impacts all the other ones. We prove the first results for existence, uniqueness, and self-consistency of the variational approximation. In turn, we derive quantitative approximation error bounds for the variational approximation to the posterior, extending prior work from the mean-field setting to the star-structured setting. We also develop a gradient-based algorithm with provable guarantees for computing the variational approximation using ideas from optimal transport theory. We explore the implications of our results for Gaussian measures and hierarchical Bayesian models, including generalized linear models with location family priors and spike-and-slab priors with one-dimensional debiasing. As a by-product of our analysis, we develop new stability results for star-separable transport maps which might be of independent interest.

Heterogeneity in efficacy is sometimes observed across baskets in basket trials. In this study, we propose a model-free clustering framework that groups baskets based on transition probabilities derived from the trajectories of treatment response, rather than relying solely on a single efficacy endpoint such as the objective response rate. The number of clusters is not predetermined but is automatically determined in a data-driven manner based on the similarity structure among baskets. After clustering, baskets within the same cluster are analyzed using a hierarchical Bayesian model. This framework aims to improve the estimation precision of efficacy endpoints and enhance statistical power while maintaining the type~I error rate at the nominal level. The performance of the proposed method was evaluated through simulation studies. The results demonstrated that the proposed method can accurately identify cluster structures in heterogeneous settings and, even under such conditions, maintain the type~I error rate at the nominal level while improving statistical power.