arXiv CS open access

Handling Missing Data in Probabilistic Regression Trees

2026 · arxiv_cs

arXiv CS · Papers · 2026

deep-learning machine-learning neural-networks

Conceptio Open Archive

Analyze & explore related research in Conceptio

Connect this record directly into your research workflow or AI agent. Explore citations across open repositories, query full text, and export structured proof bundles.

Explore in Search →Connect via MCP / API ↗

Handling Missing Data in Probabilistic Regression Trees Taiane Schaedler Prassa,∗ , Alisson Silva Neimaiera , Guilherme Pumia

arXiv:2608.06195v1 [stat.ML] 6 Aug 2026

Abstract Probabilistic Regression Trees (PRTrees) are a smooth and consistent alternative to classical regression trees, producing continuous predictions through probabilistic split assignments. This paper extends the PRTree framework to accommodate missing predictor values directly during tree construction, eliminating the need for prior imputation. Three strategies are proposed, each exploiting the available information differently: a uniform-probability approach, a partial-observation approach, and a dimension-reduced smoothing approach. These modifications are defined to preserve the fundamental probabilistic properties of the original methodology, including probability conservation and marginal compatibility, under arbitrary patterns of missing covariate values. The proposed methods are evaluated on several real-world datasets exhibiting different levels of missingness and are compared with classical regression trees. The results show that the effectiveness of probabilistic tree construction depends strongly on the treatment of missing observations. Across the considered datasets, the fill strategy emerged as the dominant modeling component, often exerting a larger influence on predictive performance than either the smoothing distribution or the proxy-selection criterion. In datasets where a substantial proportion of observations contained missing predictor values, the proposed methods frequently outperformed CART, while maintaining the interpretability and flexibility of tree-based models. Keywords: Probabilistic regression trees; missing data; nonparametric regression; statistical learning. MSC: 62G08; 62D10; 62H30; 62J02.

Introduction

Decision trees are a fundamental tool in statistical learning and machine learning, widely used for both classification and regression tasks. Their hierarchical structure partitions the predictor space through a sequence of recursive decisions, yielding models that are easy to interpret and implement. Among the various tree-based methods proposed in the literature, the Classification and Regression Trees (CART) algorithm of Breiman et al. (1984) remains one of the most widely used. CART recursively partitions the predictor space by selecting splits that locally optimize a prediction criterion. However, the resulting predictions are piecewise constant and may not adequately capture smooth relationships between predictors and the response variable (İrsoy et al., 2012; Linero and Yang, 2018). Several authors have proposed probabilistic alternatives that replace deterministic splits by smooth transition mechanisms. Examples include the Smooth Transition Regression Tree (STR-Tree) model of da Rosa et al. (2008) and the Soft Tree framework of İrsoy et al. (2012). In these approaches, observations are associated with tree nodes through probabilities rather than hard assignments, allowing predictions to vary continuously across the predictor space. Such probabilistic formulations reduce the discontinuities inherent to classical trees and provide greater flexibility for modeling complex regression functions. However, despite their methodological differences, neither STR-Trees nor Soft Trees provide native mechanisms for handling missing covariates. ∗ a

Corresponding author. This Version: August 7, 2026 Graduate Program in Statistics - Universidade Federal do Rio Grande do Sul. E-mails: taiane.prass@ufrgs.br (Prass); alissonneimaier@hotmail.com (Neimaier); guilherme.pumi@ufrgs.br (Pumi). ORCIDs: 0000-0003-3136-909X (Prass); 0000-0002-7524-0776 (Neimaier); 0000-0002-6256-3170 (Pumi).

More recently, Alkhoury et al. (2020) introduced Probabilistic Regression Trees (PRTrees), a framework that combines recursive partitioning with smooth probabilistic region assignments. In contrast to earlier probabilistic tree models, PRTrees are accompanied by a consistency result, establishing convergence of the estimated regression function to the true conditional expectation under suitable conditions. This theoretical property, together with the ability to produce continuous predictions, makes PRTrees an attractive alternative to classical regression trees. However, as originally proposed, PRTrees do not provide a native mechanism for handling missing predictor values, limiting their applicability in settings where incomplete covariate information is present. Despite these advantages, the original PRTree methodology is limited to fully observed predictor vectors. In its original formulation, no mechanism is provided for handling missing covariates during either model fitting or prediction. This limitation is particularly restrictive in practical applications, where incomplete predictor information is common. While missing values can always be addressed through complete-case analysis or external imputation procedures, such approaches require an additional preprocessing stage that is conceptually separate from the tree model itself. Handling missing data remains one of the most challenging problems in applied statistics. Missing values may reduce effective sample size, increase uncertainty, and complicate the interpretation of predictive models (Molenberghs et al., 2020). Recent theoretical work by Morvan et al. (2021) shows that impute-then-regress procedures can possess strong asymptotic optimality properties. However, the authors also point out that conditional imputations may introduce discontinuities in the estimated regression function. For methods such as PRTrees, whose primary motivation is to produce smooth predictions through probabilistic assignments, such discontinuities may be undesirable. For classical regression trees, several internal strategies have been proposed to mitigate the effects of missing predictors. In particular, the rpart implementation of CART (Therneau and Atkinson, 2019) employs surrogate splits (Breiman et al., 1984), which use alternative predictors when the primary splitting variable is unavailable. However, such mechanisms are closely tied to the deterministic structure of CART and do not naturally extend to probabilistic tree models, where observations are associated with regions through smooth probability functions rather than hard binary decisions. The absence of native missing-data mechanisms is therefore a practical obstacle to the broader use of probabilistic regression trees. The purpose of this paper is to address this limitation. Rather than proposing a new predictive paradigm or attempting to outperform existing missing-data methodologies, we develop an extension of the PRTree framework that can operate directly on incomplete predictor vectors. The proposed modifications are designed to preserve the probabilistic structure of PRTrees while eliminating the need for a separate preprocessing step devoted to handling missing covariates. To this end, we introduce three alternative strategies. The first adopts a uniform-probability assignment whenever relevant predictor information is unavailable. The second relies on the observed coordinates only, producing deterministic compatibility assessments based on the available information. The third projects both observations and regions onto the observed dimensions and applies the original probabilistic smoothing mechanism in the resulting lower-dimensional space. Although conceptually different, all three approaches are constructed so as to preserve key probabilistic properties of the original framework, including probability conservation across splits and compatibility with marginal probability calculations. A natural question is whether these extensions should be preferred over the simpler strategy of imputing missing values and subsequently fitting the original PRTree. The answer is far from straightforward. Perhaps the most delicate issue when employing imputation concerns the nature of the missing data mechanism — whether missing completely at random (MCAR), missing at random (MAR), or missing not at random (MNAR). In practice, determining the appropriate mechanism is 2

rarely a trivial matter, often being a point of discussion and a source of uncertainty. For MAR and MNAR settings, an additional layer of complexity emerges from the need to specify and model the missing data structure itself. Misspecification of either the missingness mechanism or its structural form can lead to biased or invalid imputations, which in turn may severely compromise the predictive performance of the downstream PRTree. One of the most attractive features of the proposed framework is that its non-parametric nature bypasses imputation entirely, avoiding all pitfalls related to missing data mechanism. Despite this compelling motivation, we emphasize that the goal of this paper is not establish superiority of the proposed framework over imputation-based alternatives, nor to provide an exhaustive comparative evaluation between the two paradigms. Such an interesting and important direction is deliberately left for future research. Rather, our aim is focused in providing a coherent extension of the PRTree framework that can operate directly on incomplete predictor vectors while preserving its probabilistic structure. From a practical point of view, the proposed methods offer an additional option for users who prefer to avoid the complexities and pitfalls of a separate missing-data treatment before fitting a probabilistic regression tree. The original PRTree algorithm (Alkhoury et al., 2020) is available only in Python. To make the methodology available in the R ecosystem, we developed the PRTree package (Prass and Neimaier, 2026), now available on CRAN. The package implements both the original algorithm and the missingdata extensions proposed in this paper, with the computationally intensive components implemented in FORTRAN and C. A companion paper (Prass et al., 2026) provides implementation details and a complete description of the package. The remainder of this paper is organized as follows. Section 2 reviews the original PRTree methodology. Section 3 introduces the proposed missing-data extensions and establishes their probabilistic properties. Section 4 discusses parameter estimation. Section 5 presents the empirical evaluation. Finally, Section 6 concludes with a discussion of the main findings and directions for future research.

Probabilistic Regression Trees: The Base Model

Let X be a p-dimensional input random vector that lies almost surely in a compact subspace X ⊂ Rp , and let Y be a response variable related to X through Y = f (X) + ε, where the error term ε satisfies E(ε) = 0 and E(ε2 ) < ∞. These moment conditions ensure that the conditional expectation is well-defined and that the squared-error loss is meaningful. Although the special case ε ∼ N (0, τ 2 ) is often assumed for analytical convenience, our results do not require the normality assumption. In this framework, the function f (X) = E(Y |X) is referred to as the regression function and represents the optimal predictor of Y given X under the mean squared error criterion. Tree-based methods provide a flexible approach to estimating f without imposing strong parametric assumptions. The CART algorithm (Breiman et al., 1984) builds regression trees via recursive binary partitioning of the feature space X . Starting from the full dataset at the root node, the algorithm iteratively selects the optimal splitting variable Xj and threshold z that divide the data into two subsets according to the rule Xj ≤ z versus Xj > z. The partitioning continues until a stopping criterion is met — such as a minimum node size, a threshold on error reduction, a maximum depth limit, or insufficient impurity decrease. The resulting tree partitions Rp into M disjoint regions

{R1 , · · · , RM }, with predictions given by the piecewise-constant function fCART (X) =

M X

cm I(X ∈ Rm ),

m=1

where I(·) denotes the indicator function. This greedy optimization procedure produces highly interpretable models that approximate f (X) through axis-aligned partitions of the input space. While computationally efficient, the resulting piecewise-constant predictions may struggle to capture smooth underlying relationships. PRTrees (Alkhoury et al., 2020) address this limitation by replacing the hard indicator functions in classical decision trees with smooth functions Ψ. The general prediction model takes the form fPR (X; Θ) =

M X

γm Ψ(X; Rm , σ),

(1)

m=1 M where Θ = ({Rm }M m=1 , γ, σ), comprises the partition regions {Rm }m=1 , the region-specific weights p M γ ∈ R , and the noise vector σ ∈ R+ controlling the smoothness of the soft assignments. For any X ∈ Rp , Rm ⊂ Rp , and σ ∈ Rp+ , the association functions Ψ are defined as

Ψ(X; Rm , σ) =

Y p k=1

−1 Z σk

vp − X p v1 − X 1 ϕ ,··· , σ1 σp Rm

dv,

(2)

where ϕ is a probability density function belonging to L2 , continuously differentiable (C 1 ), and whose Fourier transform is nonzero on Rp . The function Ψ defines smooth, probabilistic memberships between points and regions, avoiding the abrupt boundaries of conventional trees. As a result, predictions vary continuously with X while retaining interpretability. The classical regression tree is recovered as the special case Ψ(X; Rm , σ) = I(X ∈ Rm ). In practice, selecting a suitable function ϕ is a nontrivial task. Prior knowledge of the noise distribution can help narrow down viable candidates for ϕ in (2), and the choice can be further refined using cross-validation. The standard PRTree training algorithm — which assumes fully observed feature vectors — is described in detail by Alkhoury et al. (2020), who also establish its key theoretical properties, most notably consistency — the guarantee that the estimated regression function converges to the true conditional expectation as sample size increases. The consistency of PRTree is a remarkable result — not known for earlier probabilistic tree methods such as STR-Trees or Soft Trees — making PRTree particularly attractive for regression tasks. However, the original formulation provides no mechanism for handling missing covariates, restricting its use to complete datasets. In the next sections, we present an extension that preserves the probabilistic structure and interpretability of PRTrees while accommodating missing covariates directly.

Handling Missing Values in PRTrees

In what follows, we present the proposed modifications to the PRTree algorithm for missing data handling. Each proposed modification is constructed in such a way as to guarantee the production of well-defined association values, regardless of the missing data mechanism.

3.1

Notation and Setup

To fix notation, let S ⊆ {1, · · · , p} denote an arbitrary subset of indices. For any X ∈ Rp , denote by X |S the subvector containing only the coordinates indexed by S. For any rectangular region 4

Q Q Rm = pj=1 Rmj ⊂ Rp , denote by Rm|S = j∈S Rmj its projection onto these coordinates, and for any vector σ ∈ Rp+ , let σ |S be its restriction to the components in S. A notational convention used throughout the paper requires explicit mention: I(A) denotes the indicator function of event A and, by convention, is always evaluated first in any expression where it appears. This convention ensures unambiguous interpretation of formulas involving products of indicators with other mathematical objects. The key to our missing-data extension is that the PRTree smoothing function (2) naturally induces a probability measure over the set of regions {Rm }M m=1 , conditioned on the observed input X: P (Rm | X) := Ψ(X; Rm , σ),

1 ≤ m ≤ M.

Let Ψ∗ denote the adapted smoothing function in the presence of missing values. The corresponding probability measure is then P ∗ (R | X) := Ψ∗ (X; R, σ). Any valid extension to handle missing values must preserve two fundamental properties of the original construction: P1: Marginal probability compatibility for unbounded coordinates. If a region Rm does not impose any restriction on coordinate k, that is, Rmk = R for some 1 ≤ k ≤ p, then P (Rm | X) = P (Rm|S | X |S ),

with S = {1, · · · , p} \ {k}.

P2: Probability conservation under splitting. If a parent region RF is partitioned into two child regions RFL and RFR , then the law of total probability must hold: P (RFL | X) + P (RFR | X) = P (RF | X),

for all X ∈ X .

Property P1 ensures that coordinates which are unrestricted (i.e., have not been bounded by previous splits) do not affect the probability computation. Consequently, if a coordinate is unbounded, missing values in that coordinate cannot influence the probability assigned to the region. Property P2 guarantees that the total probability mass assigned to a region is exactly redistributed among its children after a split, maintaining coherence as the tree grows.

3.2

Three Strategies for Handling Missing Values

Our approach introduces three alternative strategies for handling missing values, controlled in the PRTree package by the parameter fill type (taking values 0, 1 or 2). Each strategy specifies how the association function Ψ∗ should be computed when some coordinates of X are missing: 0: Uniform probability method. When any value is missing, this strategy ignores all observed information and assigns equal weight to both child regions. It is the most conservative approach, providing robustness at the cost of discarding potentially useful data. 1: Partial observation approach. This strategy performs a “hard” assignment based on the observed coordinates: it returns 1 if the observed values are compatible with the region and 0 otherwise. Smoothing via Ψ∗ is applied only to fully observed cases, offering a balance between robustness and informativeness. 2: Smoothed projection technique. Here, both the observation and the region are projected onto the dimensions corresponding to the observed coordinates, and the adapted smoothing function Ψ∗ is applied in this reduced space. This makes full use of the available information while preserving the smooth structure of PRTrees.

This flexible framework allows practitioners to choose an approach aligned with their application’s priorities — favoring simplicity, balanced trade-offs, or maximal informativeness — without compromising the integrity of probability estimation. The adapted function Ψ∗ (and, consequently, the corresponding probability measure P ∗ ) is defined recursively as follows. Definition 3.1. Let Ψ be defined by (2). Given X ∈ Rp and any region R ⊆ Rp : Base case: if R = Rp (root node), then Ψ∗ (X; R, σ) = 1. Recursive step: if R ∈ {RFL , RFR } is a child region resulting from the split of a parent region RF , first define the set of indices S(X, R) as follows: S(X, R) := j ∈ {1, · · · , p} : RFL j ⊊ R and Xj is non-missing , i.e., the indices j such that the corresponding variable has been used in a previous split and Xj is non-missing. Next, define a proxy function H as:   Ψ(X; R, σ), if X is fully observed,       1, if fill type = 0 or S(X, R) = ∅, H(X; R, σ) = (3)   I(X ∈ R ), if fill type = 1 and S(X, R) = ̸ ∅,  |S(X,R) |S(X,R)     Ψ(X |S(X,R) ; R|S(X,R) , σ |S(X,R) ), if fill type = 2 and S(X, R) ̸= ∅. Then Ψ∗ is given recursively by Ψ∗ (X; R, σ) =

H(X; R, σ)Ψ∗ (X; RF , σ) I Ψ∗ (X; RF , σ) > 0 . H(X; RFL , σ) + H(X; RFR , σ)

From the definition of S, note that S(X, RFL ) = S(X, RFR ), and S(X, RF ) ⊆ S(X, RFL ), with equality if, and only if, either the splitting coordinate j was already in S(X, RF ) or Xj is missing. This leads to a fundamental property shared by all strategies: if a partition RF = RFL ∪RFR is created by splitting on feature j, then H(X; RFL , σ) = H(X; RFR , σ) whenever Xj is missing. When missing values occur in other coordinates but Xj is observed, each strategy behaves as previously described: fill type = 0 ignores all observed information, fill type = 1 performs a hard assignment based on observed coordinates, and fill type = 2 projects onto the observed dimensions before applying Ψ. All three methods provide well-defined probability estimates for any missing-data pattern and reduce to the standard PRTree algorithm when no values are missing. Lemma 3.1 gives a characterization of Ψ∗ , useful for practical implementation when there are missing values in the coordinates used to build the tree. Lemma 3.1 (Product representation). Fix a node Rm in the tree and X ∈ Rp . Denote by RF m the parent node of Rm and by ℜm the set of internal nodes on the path from the root (including) to Rm . Let S := S(X, Rm ) and ℜm|S ⊆ ℜm be the subset of nodes that were split using a coordinate that belongs to S. If Ψ∗ (X; RF m , σ) ̸= 0, then Ψ∗ (X; Rm , σ) = 2−|ℜm \ℜm|S |

H(X; Rc(RF ) , σ) , H(X; RFL , σ) + H(X; RFR , σ)

(4)

RF ∈ℜm|S

where, for any RF ∈ ℜm with children RFL and RFR , c(RF ) ∈ {FL , FR } is the child index on the path to Rm . 6

Proof. The result follows immediately by recursively applying the definition of Ψ∗ . An immediate consequence of Lemma 3.1 is that, for any strategy adopted (that is, any fill type), if p = 1 and X is missing, Ψ∗ (X, Rm , σ) = 2−d(m) , where d(m) is the depth of Rm . In the particular case when P (RF |X) = Ψ(X; RF , σ) can be factorized as P (RF |X) = Ψ(X; RF , σ) =

p Y

Ψ(Xj ; RF j , σj ) =

j=1

p Y

P (RF j |Xj ),

j=1

then the same holds upon replacing Ψ with Ψ∗ and the child probabilities satisfy P ∗ (RFL |X) =

P ∗ (RFL j |Xj )P ∗ (RF |X) P ∗ (RFL j |Xj ) + P ∗ (RFR j |Xj )

and

P ∗ (RFR |X) =

P ∗ (RFR j |Xj )P ∗ (RF |X) , P ∗ (RFL j |Xj ) + P ∗ (RFR j |Xj )

whenever P ∗ (RF |X) > 0. This means that, at each step, we only need to compute the marginal probabilities for the splitting feature j and then scale them by the probability of the parent node. This fact speeds up computations.

Parameter estimation

Given a training sample {(X i , Yi )}ni=1 , with X ∈ Rp , Y ∈ R, and in accordance with the empirical risk minimization principle with a quadratic loss, the estimation procedure for probabilistic regression trees aims at finding the parameters Θ solutions of ( n 2 ) M X X argmin Yi − γm Pim , with Pim := Ψ∗ (X i ; Rm , σ), (5) Θ∈Ξ

m=1

i=1

where Ξ=

o p p M {Rm }M , γ, σ : M ∈ N\{0}, R ⊆ R , γ ∈ R , σ ∈ R m m=1 + .

The n × M matrix P , with entries Pim , thus encodes the relations between each training example P X i and each region Rm . It is such that 0 ≤ Pim ≤ 1 and M m=1 Pim = 1, for all 1 ≤ i ≤ n. The estimation of the different parameters in Θ alternates in between region and weight estimates, as in standard regression trees, until a stopping criterion is met. During this process, the number of regions is increased by one at each loop and the matrix P and the weights γ are gradually updated. The vector σ can either be based on a priori knowledge or be estimated through a grid search on a validation set. Estimating γ. Given the regions {Rm }M m=1 and the vector σ, minimizing (5) with respect to γ leads to the least square estimator ( n 2 ) M n o X X γ̂ = argmin Yi − γm Pim = argmin ||Y − P γ||2 , (6) γ∈RM

i=1

γ∈RM

m=1

where Y = (Y1 , · · · , Yn )′ . If P ′ P is not singular, the solution is unique and it is given by γ̂ = (P ′ P )−1 P ′ Y . Estimating {Rm }M m=1 . Assume that M regions, referred to as current regions, have already been identified, meaning that the current tree has M leaves. As in standard regression trees, each current region Rm , 1 ≤ m ≤ M , can be decomposed into two sub-regions RFL and RFR with respect to a coordinate 1 ≤ j ≤ p and a splitting point t (threshold) that minimizes (5). Each possible split leads 7

to the update of P , that now belongs to Rn×(M +1) (the space of n by M + 1 real matrices), and γ, that now belongs to RM +1 . The corresponding γ̂ is obtained through (6). To make explicit the dependence of these quantities on the region, the splitting variable and the threshold, we shall use the notation P (m) (j, t) and γ̂ (m) (j, t). The best split for the current region Rm solves ( n 2 ) M +1 X X (m) (m) argmin Yi − γ̂k (j, t)Pik (j, t) , (m)

(j,t)∈J ×Tj

i=1

k=1 (m)

where J = {j ∈ N : 1 ≤ j ≤ p} and Tj denotes the set of splitting points for region Rm and variable j (more precisely, the set of middle points of the observations from Rm projected on the jth coordinate). At each step, the split that leads to the smallest mean square error is selected. Note that, unlike in CART, a region Rm in a PRTree need not contain any observations. This is because Ψ∗ (X; Rm , σ) defines a conditional probability (or association measure) that does not require X to actually fall within the region – it only measures how compatible X is with that region. Nevertheless, the observed data still provide a natural set of candidate thresholds for splitting, and we adopt the standard midpoint rule for computational efficiency and consistency with tree-based methods. Stopping. The process terminates when the maximum number of regions (Mmax ) is reached, the reduction in loss ∆L falls below a threshold ϵ, or other stopping criteria (e.g., depth, minimum node size) are met. Alkhoury et al. (2020) show that the PRTree learned from a training set of size n, with no missing (n) data, denoted fˆPR , is consistent in the sense that 2 (n) lim E fˆPR (X) − E(Y |X) = 0. n→∞

The consistency is a desirable theoretical development, but the proof of such result is quite involved even in the standard framework without missing data. Extending the results to account for the missing data handling mechanism in (3) is a non-trivial task and will be explored in future works.

Empirical Application

The objective of this empirical study is to evaluate the practical behavior of the proposed missingdata extensions in real regression problems containing naturally occurring missing covariates. Rather than establishing superiority over alternative missing-data methodologies, the goal is to assess whether probabilistic regression trees remain a viable modeling option when incomplete predictor information is handled internally by the algorithm. We consider a collection of benchmark regression datasets exhibiting heterogeneous sample sizes, predictor dimensions, and missing-data patterns. The proposed extensions are evaluated through cross-validation and compared with classical CART regression trees, which provide a natural baseline for assessing the practical impact of probabilistic region assignments and internal missing-data handling. Particular attention is given to the influence of the missing-data strategy, the proxy selection criterion, and the smoothing distribution used by the probabilistic model. The empirical analysis focuses on predictive accuracy and robustness across datasets with naturally occurring missing values. The objective is not to identify a universally superior specification, but rather to investigate how the proposed extensions behave under different modeling choices and whether they provide a practical 8

alternative for users who prefer not to perform a separate missing-data treatment before fitting a probabilistic regression tree.

5.1

Main empirical question

The empirical study addresses two questions. The first is whether PRTrees remain competitive with classical CART regression trees when applied to datasets containing naturally occurring missing covariate values. This comparison is based on predictive performance measured through RMSE and relative predictive error. The second question concerns the effect of the main modeling choices available in the proposed methodology. In particular, the experiments investigate how predictive performance varies with the missing-data extension (fill type), the proxy-selection criterion (proxy crit), and the smoothing density specification ϕ. The objective is to identify which combinations of these components are associated with the best predictive performance across the datasets considered in the study.

5.2

Data Acquisition and Preparation

The empirical study was conducted using real-world regression datasets obtained primarily from the UCI Machine Learning Repository. Candidate datasets were initially identified through the ucimlrepo R package (Balamuta and Truong, 2024), which provides programmatic access to dataset metadata and facilitates large-scale screening of regression problems. To ensure a diverse yet computationally manageable benchmark, datasets were filtered according to basic characteristics such as sample size, predictor dimension, and response type. An additional benchmark dataset, Ozone, available in the mlbench package (Leisch and Dimitriadou, 2026), was also included in the analysis. The final collection consists of five regression datasets: Communities and Crime (Crime), Productivity Prediction of Garment Employees (Garment), Auto MPG, Automobile, and Ozone. These datasets exhibit substantial heterogeneity in both sample size and predictor dimension, while all containing naturally occurring missing values in the predictor variables. This diversity provides a realistic setting for evaluating probabilistic regression trees under incomplete covariate information. Prior to model fitting, observations with missing response values were removed whenever necessary. Missing values in the predictor variables were preserved and handled directly by the proposed methods. Categorical predictors were excluded from the analysis so that all models were fitted using numerical covariates only. This preprocessing step avoids the introduction of additional variability associated with categorical encoding schemes and allows the empirical evaluation to focus exclusively on the effects of missing predictor information. Table 1 summarizes the main characteristics of the datasets considered in the empirical study. Table 1: Summary of the datasets used in the empirical evaluation: n and p denote the dataset dimensions; mglobal is the percentage of missing predictor cells among all n × p predictor values; nmiss is the number of observations containing at least one missing predictor value.

Dataset

Source

UCI ID

Response range

Crime Garment Auto MPG Ozone Automobile

UCI UCI UCI mlbench UCI

183 597 9 10

[0, 1] [0.23, 1.12] [9, 46.6] [1, 38] [−2, 3]

mglobal (%)

nmiss (%)

1994 1197 398 361 205

122 10 7 9 17

15.15 4.23 0.22 6.03 1.69

1675 (84.00%) 506 (42.27%) 6 (1.51%) 158 (43.77%) 46 (22.44%)

The datasets differ substantially not only in their overall proportion of missing predictor values, but also in the fraction of observations affected by missingness. For example, although the Ozone dataset contains only 6.03% missing predictor cells, approximately 44% of its observations contain at least one missing predictor value. A similar phenomenon occurs in the Crime dataset, where 15.15% of the predictor cells are missing but nearly 84% of the observations are affected. By contrast, only six observations in Auto MPG contain missing predictor values. These differences are expected to influence the potential benefits of the probabilistic mechanisms introduced in PRTree, since they determine how frequently incomplete observations must be propagated through the tree.

5.3

Experimental Configuration

Shared configuration. Predictive performance was assessed using the Root Mean Squared Error (RMSE) under 15-fold cross-validation. The same data partitions were used across all competing methods and experimental configurations in order to ensure direct comparability between predictive results. In addition to predictive accuracy, computational time was recorded for all methods and configurations. Reported runtimes correspond to the complete model-fitting procedure associated with each method. CART configuration. The rpart implementation was used with internal pruning cross-validation configured through xval = 10. The experimental configuration employed complexity parameter cp = 0.001, minimum splitting size supplied through minsplit, defined as approximately 1% of the effective training sample size, and maximum tree depth supplied through maxdepth = 30 (the maximum allowed by the package). CART control parameters were initialized from the default rpart.control() configuration and updated whenever compatible arguments were supplied through the experimental configuration. Tree pruning was subsequently performed using the complexity parameter associated with the minimum internal cross-validated prediction error in the corresponding cost-complexity table. PRTree configuration. Tree construction was controlled through cp = 0.001, meaning that a split was accepted only if the corresponding reduction in MSE exceeded 0.1%. The maximum number of terminal regions was fixed at max terminal nodes = 100, while prop min = 0.01 yielded a minimum node size of approximately 1% of the effective training sample size. Additional probabilitybased stopping restrictions were imposed through perc x = 0.1 and p min = 0.05, corresponding to the default values of the PRTree package. Candidate split search employed by node = TRUE and n candidates = max(⌊p/2⌋, 3), where p denotes the number of predictors in the corresponding dataset. Smoothing-parameter selection employed grid size = 8, corresponding to the default value of the PRTree package and generating eight positive candidate values for each component of σ. The candidate grid was augmented with the additional value tiny sigma = 0, allowing the limiting case of vanishing smoothing parameters to recover the hard-split behavior of classical CART trees. The parameter vector σ was selected using an internal validation sample approximately equal in size to a single cross-validation fold. For the smoothing density function ϕ, nine specifications were considered: Gaussian, Student-t, lognormal, and Gamma. The Student-t specification employed degrees of freedom df ∈ {3, 5}. The lognormal specification employed dispersion parameters sdlog ∈ {0.5, 1, 2}. The Gamma specification employed shape parameters shape ∈ {0.5, 2, 5}. The empirical evaluation considered the three proposed missing-data extensions associated with fill type values 0, 1, and 2, together with the three Stage 1 proxy-selection criteria "mean", "var", and "both". Combining the three missing-data strategies, three proxy criteria, and nine smoothing density specifications resulted in 81 distinct PRTree configurations for each dataset.

5.4

Relative predictive error

The datasets considered in the empirical study correspond to distinct prediction problems and differ substantially in scale and predictive characteristics. Consequently, direct comparisons of raw RMSE values across datasets are of limited interest. For this reason, each PRTree configuration was evaluated relative to the corresponding CART performance obtained on the same dataset and cross-validation fold. Let RMSEd (m, k) denote the root mean squared error obtained on dataset d by method m on fold k, and let RMSECART (k) denote the RMSE obtained by CART on the same fold. Relative predictive d error was defined as RMSEd (m, k) − RMSECART (k) d ∆d (m, k) = . CART RMSEd (k) Under this definition, ∆d (m, k) = 0 indicates identical predictive performance, negative values indicate that the corresponding PRTree configuration outperformed CART, and positive values indicate inferior performance relative to CART. The quantity ∆d (m, k) measures predictive performance relative to a common baseline within each dataset. Although the datasets remain inherently different and should not be interpreted as directly comparable prediction tasks, the relative error makes it possible to assess whether particular modeling choices, such as the missing-data handling strategy, proxy-selection criterion, or smoothing specification, tend to improve or deteriorate predictive performance across the collection of datasets considered in this study. Moreover, because all methods and configurations were evaluated using identical cross-validation partitions within each dataset, relative predictive errors can be compared directly across competing PRTree configurations.

5.5 5.5.1

Results Predictive performance relative to CART

To better understand the contribution of the probabilistic components introduced in PRTree, Figures 1 and 2 compare its predictive performance against the classical CART algorithm. Throughout this section, negative values of ∆ indicate that PRTree achieved lower prediction error than CART, whereas positive values indicate the opposite. Figure 1 provides an overall summary of the effect of the fill strategy. For each dataset, the boxplots display the distribution of ∆ values obtained across all smoothing distributions, proxy-selection criteria, and cross-validation folds associated with a given value of fill type. The horizontal red line corresponds to ∆ = 0, representing equal predictive performance between PRTree and CART. To improve readability, the vertical axis was truncated at a fixed positive threshold. Whenever observations exceeded this limit, the symbol n▲ was added above the corresponding boxplot, where n denotes the number of omitted observations. Therefore, these annotations identify configurations that produced larger prediction errors than those displayed within the plotting range. While the boxplots summarize the overall behavior of each fill strategy, they do not distinguish between smoothing distributions or proxy-selection criteria. A more detailed view is provided by Figure 2, which reports the empirical win rate of every evaluated configuration. In Figure 2, rows correspond to combinations of fill strategy and proxy-selection criterion. The labels f0, f1, and f2 denote fill type 0, 1, and 2, respectively. The letters B, M, and V refer to the combined, meanbased, and variance-based proxy criteria. Columns correspond to the smoothing distributions: N denotes the Gaussian distribution, t(3) and t(5) denote Student-t distributions, LN denotes lognormal distributions, and G denotes Gamma distributions. Each cell contains the proportion of crossvalidation folds in which PRTree achieved lower prediction error than CART. Values above 50% 11

Auto MPG 0.6

Automobile

9▲

Crime

21▲

33▲

22▲

Garment

Ozone

3▲

19▲

5▲

0.3

0.0

-0.3

Figure 1: Distribution of relative predictive errors grouped by fill type. Negative values indicate that PRTree outperformed CART.

indicate that PRTree outperformed CART more often than not, whereas values below 50% indicate the opposite. The most important result emerging from Figures 1 and 2 is the dominant role played by the fill strategy. Across all datasets, changes in fill type produce substantially larger variations in predictive performance than changes in either the smoothing distribution or the proxy-selection criterion. Configurations sharing the same fill strategy often exhibit similar win rates even when different distributions and proxy criteria are employed, whereas changing the fill strategy can completely alter the relative performance of PRTree with respect to CART. The influence of the proxy-selection criterion appears comparatively small. Within a fixed fill strategy, the rows associated with the mean-based, variance-based, and combined criteria usually display similar win rates, and no criterion consistently dominates across datasets. This suggests that the proxy criterion acts primarily as a secondary tuning parameter whose influence is considerably weaker than that of the fill mechanism. The role of the smoothing distribution is more subtle. The heatmaps indicate that distributional choices can affect predictive performance, but their influence is highly dependent on both the selected fill strategy and the characteristics of the dataset. In datasets such as Crime and Ozone, where fill type = 1 and 2 clearly outperform fill type = 0, the improvement is observed across a broad range of distributions. In these cases, changing the distribution typically modifies the magnitude of the gain rather than altering the overall ranking of the configurations. This suggests that most of the predictive improvement is attributable to the fill mechanism itself, while the smoothing distribution acts primarily as a refinement component. Nevertheless, some distributional patterns can still be identified. Student-t distributions frequently appear among the strongest configurations, particularly when combined with the most successful fill strategies. Although the experiments do not directly explain this behavior, the results suggest that heavier-tailed smoothing mechanisms may provide a more effective allocation of incomplete observations near splitting boundaries. At the same time, the superiority of Student-t distributions is not universal, and competitive results are frequently obtained using log-normal and Gamma distributions. Consequently, the choice of smoothing distribution should be viewed as a secondary tuning decision: it can improve performance once an appropriate fill strategy has been selected, but it rarely compensates for an inadequate fill mechanism. Auto MPG. The Auto MPG dataset illustrates a relatively stable scenario. The three fill strategies produce similar boxplots, all centered slightly below zero, indicating that PRTree tends to outperform 12

Auto MPG

Automobile

f2-V

73%

80%

53%

60%

53%

67%

47%

33%

f2-V

27%

33%

40%

33%

47%

53%

40%

27%

f2-M

73%

80%

53%

60%

53%

67%

47%

33%

f2-M

20%

40%

47%

67%

47%

27%

f2-B

73%

80%

53%

60%

53%

67%

47%

33%

f2-B

27%

33%

40%

33%

47%

53%

40%

27%

f1-V

67%

53%

67%

47%

80%

73%

60%

40%

f1-V

47%

27%

47%

33%

47%

60%

47%

27%

f1-M

67%

53%

67%

47%

80%

73%

60%

40%

f1-M

33%

20%

33%

53%

47%

53%

60%

33%

40%

f1-B

67%

53%

67%

47%

80%

73%

60%

40%

f1-B

47%

27%

47%

33%

47%

60%

47%

27%

f0-V

73%

60%

73%

67%

73%

60%

67%

60%

53%

f0-V

40%

47%

53%

33%

40%

47%

60%

40%

47%

f0-M

73%

60%

73%

67%

73%

60%

67%

60%

53%

f0-M

33%

40%

53%

40%

27%

47%

40%

47%

53%

f0-B

73%

60%

73%

67%

73%

60%

67%

60%

53%

f0-B

47%

53%

33%

40%

47%

60%

40%

47%

t(3)

t(5) LN(0.5)

t(3)

t(5) LN(0.5)

LN(1)

LN(2)

G(0.5)

G(2)

G(5)

Crime

LN(1)

LN(2)

G(0.5)

G(2)

G(5)

Garment

f2-V

100%

67%

80%

73%

100%

73%

80%

f2-V

40%

13%

27%

40%

67%

13%

f2-M

100%

f2-B

93%

100%

73%

87%

100%

67%

80%

93%

80%

87%

67%

f2-M

53%

60%

67%

20%

47%

40%

67%

33%

13%

80%

87%

80%

73%

f2-B

40%

33%

40%

20%

47%

60%

13%

f1-V

80%

87%

53%

87%

67%

73%

67%

f1-V

40%

47%

27%

40%

27%

f1-M f1-B

73%

87%

93%

80%

93%

80%

53%

60%

73%

60%

53%

f1-M

53%

67%

53%

27%

53%

47%

53%

47%

20%

67%

73%

67%

73%

f1-B

40%

53%

33%

47%

40%

33%

f0-V

33%

20%

27%

33%

20%

f0-M

f0-B

f0-M

47%

40%

33%

20%

33%

27%

40%

13%

f0-B

33%

20%

27%

33%

40%

20%

t(3)

t(5) LN(0.5)

t(3)

t(5) LN(0.5)

LN(1)

0.25

0.50

LN(1)

LN(2)

G(0.5)

G(2)

G(5)

LN(2)

G(0.5)

G(2)

G(5)

Ozone f2-V

93%

100%

67%

100%

67%

60%

100%

73%

f2-M

87%

93%

60%

80%

67%

60%

73%

60%

f2-B

93%

100%

67%

100%

60%

100%

73%

f1-V

73%

67%

73%

60%

73%

60%

47%

60%

73%

f1-M

87%

73%

67%

60%

53%

40%

60%

67%

73%

f1-B

87%

53%

67%

80%

60%

47%

60%

73%

f0-V

20%

27%

20%

27%

20%

27%

20%

f0-M

27%

13%

27%

20%

27%

f0-B

20%

27%

20%

t(3)

t(5) LN(0.5)

PRTree win rate

0.00

27% LN(1)

20% LN(2)

27% G(0.5)

27% G(2)

0.75

1.00

20% G(5)

Figure 2: Empirical win rate of each PRTree configuration relative to CART.

CART regardless of the specific probabilistic configuration. The heatmap reinforces this conclusion, with most configurations achieving win rates above 50% and several exceeding 70%. Although fill type=2 appears to provide the strongest overall performance, the differences among fill strategies remain moderate. One possible explanation is the extremely low prevalence of missing values in this dataset. As shown in Table 1, only six observations (1.51% of the sample) contain at least one missing predictor value, corresponding to merely 0.22% of all predictor cells. Consequently, the probabilistic components of PRTree are activated only rarely, limiting the potential impact of the fill strategy on predictive performance. This is consistent with the relatively small differences observed among the competing configurations. Automobile. The Automobile dataset presents the most balanced comparison among all considered datasets. The boxplots exhibit substantial variability and a considerable number of observations above the plotting threshold, indicating that some configurations can perform substantially worse 13

than CART. Consistently, the heatmap shows that most win rates lie between 30% and 60%, with only a few configurations exceeding these limits. No clear pattern emerges with respect to either fill strategy or smoothing distribution. This absence of a dominant structure suggests that neither method consistently outperforms the other in this setting. The characteristics of the dataset reported in Table 1 are consistent with this behavior. The overall proportion of missing predictor cells is relatively small (1.69%), and only 22.44% of the observations contain at least one missing predictor value. Consequently, the probabilistic mechanisms affect only a limited portion of the data, reducing the opportunity for PRTree to obtain substantial gains over CART. Moreover, Automobile is the smallest dataset considered in the study (n = 205), making performance estimates more sensitive to variations in the training sample. The combination of a limited amount of missing information and a relatively small sample size may therefore help explain the absence of a clear advantage for either method. Communities and Crime. The Crime dataset provides the strongest evidence of the importance of the fill mechanism. The contrast between fill type = 0 and the remaining strategies is striking. Every configuration associated with fill type = 0 loses to CART in all cross-validation folds, producing win rates equal to 0%. However, once the fill mechanism is modified, the behavior changes completely. Most configurations associated with fill type = 1 and 2 achieve win rates above 80%, and several reach 100%. Such a dramatic transition cannot be attributed to the choice of distribution or proxy criterion, since it occurs consistently across almost all configurations within each fill strategy. Instead, it strongly suggests that the way incomplete observations are propagated through the tree is the dominant factor determining predictive performance in this dataset. The characteristics of the dataset reported in Table 1 help explain this behavior. The Crime dataset exhibits both the largest proportion of missing predictor cells (15.15%) and the largest proportion of affected observations (84.00%). As a result, the treatment of incomplete observations becomes a central component of the learning procedure. Under these conditions, differences in the fill mechanism have a direct impact on a substantial fraction of the training data, which helps explain the dramatic performance gap observed between fill type = 0 and the remaining strategies. Productivity Prediction of Garment Employees. The Garment dataset exhibits the opposite behavior of Crime and Ozone. Most configurations achieve win rates below 50%, indicating that CART generally provides lower prediction errors. The heatmap reveals that this pattern is remarkably consistent across distributions and proxy-selection criteria, suggesting that the observed disadvantage is not driven by a particular probabilistic choice. Unlike Crime and Ozone, where the probabilistic treatment of missing observations produces substantial gains, the Garment dataset provides little evidence that the additional flexibility introduced by PRTree translates into systematic improvements in predictive performance. Although a few isolated configurations achieve competitive results, CART remains the most reliable choice overall. Interestingly, this result shows that a large number of incomplete observations does not automatically imply that probabilistic treatment will improve predictive performance. Although more than 42% of the observations contain missing predictor values, CART remains competitive and often superior. This suggests that the effectiveness of the probabilistic mechanisms depends not only on the amount of missing information, but also on how that missingness interacts with the underlying structure of the prediction problem. Ozone. The Ozone dataset displays a pattern similar to that observed for Crime, although less extreme. A clear monotonic improvement is visible when moving from fill type = 0 to fill type = 2. The boxplots progressively shift from predominantly positive relative errors to predominantly 14

negative ones, while the heatmap reveals a corresponding increase in win rates. Several configurations associated with fill type = 2 achieve win rates above 90%, whereas most configurations associated with fill type = 0 remain below 30%. The consistency of this pattern across smoothing distributions and proxy-selection criteria reinforces the conclusion that the fill mechanism is the primary driver of predictive performance. The dataset characteristics reported in Table 1 provide a plausible explanation for this behavior. Although only 6.03% of the predictor cells are missing, approximately 44% of the observations contain at least one missing predictor value. Consequently, a large fraction of the sample is affected by missingness, making the treatment of incomplete observations an important component of the learning process. This is consistent with the strong improvements observed when moving from fill type = 0 to the more sophisticated fill strategies. Taken together, the results reveal a clear hierarchy among the probabilistic design choices. The fill strategy is by far the most influential component, frequently determining whether PRTree outperforms or underperforms CART. The smoothing distribution occupies an intermediate position, refining the behavior of an already suitable fill strategy, whereas the proxy-selection criterion generally exerts only a minor influence on the final results. From a practical perspective, model selection should therefore proceed hierarchically: first identifying a promising fill strategy and only then refining the smoothing distribution and proxy-selection criterion. 5.5.2

Computational cost

In addition to predictive accuracy, computational cost was recorded for all methods and configurations considered in the empirical study. The reported runtimes were obtained on a desktop computer equipped with an Intel® Core™ i5-4590 CPU (3.30 GHz) and 16 GB of RAM, running Windows 10. Table 2 summarizes the observed runtimes for CART and PRTree across all datasets. Along with the computational results, the table reports the sample size n, predictor dimension p, and the number of candidate splits fully evaluated at each node, ncand = max{p/2, 3}. These quantities are useful for interpreting the observed runtimes, since the computational burden of PRTree depends both on the amount of data processed at each node and on the number of candidate splits receiving a full MSE evaluation. Table 2: Computational cost of CART and PRTree. For each dataset, the table reports the sample size n, predictor dimension p, the number of candidate splits fully evaluated at each node (ncand ), the total runtime of CART (in seconds), the cumulative runtime required to evaluate all 81 PRTree configurations (in minutes), and the minimum and maximum runtimes observed among individual PRTree configurations (in seconds).

Dataset Crime Garment Auto MPG Ozone Automobile

ncand

CART (s)

PRTree (min)

PRTree min (s)

PRTree max (s)

1994 1197 398 361 205

122 10 7 9 17

61 5 4 4 8

8 1 1 1 1

2106.05 100.60 56.23 55.45 63.10

308 22 5 10 9

7966 217 112 118 116

The computational burden of PRTree is substantially larger than that of CART. This increase is expected because each PRTree fit requires repeated evaluations of probabilistic node assignments and an internal smoothing-parameter selection procedure. For every configuration, nine candidate values of the smoothing parameter are evaluated, corresponding to the eight values generated by grid size = 8 together with the additional value tiny sigma = 0. Since the tree structure depends on the

smoothing parameter, each candidate value requires a complete reconstruction of the tree. An additional source of computational cost arises from the split-selection procedure. At each node, all possible splits are first ranked using the selected proxy criterion. The best ncand = max{p/2, 3} candidate splits are then fully evaluated using the MSE criterion. Consequently, the computational effort associated with split selection depends both on the number of observations contributing to the node statistics and on the number of candidate splits receiving a full MSE evaluation. The runtimes reported in Table 2 suggest that both sample size and ncand play an important role in determining computational cost. The Crime dataset required more than 35 hours to evaluate all configurations, whereas the remaining datasets required between approximately one and two hours. This behavior is consistent with the fact that Crime combines the largest sample size (n = 1994) with the largest value of ncand (61). By contrast, although Automobile evaluates more candidate splits than Garment (8 versus 5), its substantially smaller sample size (205 versus 1197 observations) results in a lower overall runtime. Taken together, these results suggest that the computational burden of PRTree is driven by the interaction between the amount of data processed and the number of candidate splits receiving a full MSE evaluation. Missingness may also affect computational cost because incomplete observations require additional processing during the calculation of the proxy measures used to rank candidate splits. As the proportion of missing values increases, these calculations become more frequent throughout tree construction. Nevertheless, the runtimes reported in Table 2 do not suggest that missingness alone is the dominant driver of computational cost. The largest differences in runtime are more closely associated with the combined effects of sample size and the number of candidate splits receiving a full MSE evaluation. To investigate how the probabilistic design choices affect computational cost, Figure 3 reports the average runtime of each evaluated configuration. The notation follows the same convention adopted in Figure 2: rows correspond to combinations of fill strategy and proxy-selection criterion, while columns correspond to the smoothing distributions and their associated parameters. Several patterns emerge from Figure 3. First, the choice of smoothing distribution has a substantial impact on runtime. Across all datasets, the Student-t distributions are consistently the most computationally demanding alternatives. The configurations associated with t(3) and t(5) systematically require longer execution times than the corresponding Gaussian, log-normal, and Gamma configurations. In contrast, the log-normal distributions are typically among the fastest options, while the Gamma distributions occupy an intermediate position. Second, the fill strategy also affects computational cost. Within a given dataset, runtimes generally increase when moving from fill type = 0 to fill type = 2. This pattern is particularly visible for the Crime dataset, where some fill type = 2 configurations require several times more computation than their fill type = 0 counterparts. Similar, although less pronounced, patterns can also be observed in the remaining datasets. Unlike the large differences in runtime observed across datasets, which are largely associated with predictor dimension and the corresponding value of ncand , these within-dataset variations are likely related to the additional probabilistic allocations required by more sophisticated fill mechanisms. Comparing Figures 2 and 3 reveals that the most computationally expensive configurations are not necessarily those with the strongest predictive performance. In the Crime and Ozone datasets, many of the best-performing configurations are indeed associated with fill type = 2, which is also the most expensive fill strategy. However, within a fixed fill strategy, the distributions that require the longest runtimes do not systematically achieve the highest win rates. Student-t distributions are often among the slowest configurations, yet Gamma and log-normal distributions frequently attain comparable predictive performance at substantially lower computational cost. Taken together, the results suggest that different factors govern runtime at different levels of the 16

Auto f2-V

0.6

Automobile 1.6

1.7

0.2

0.8

0.3

0.1

f2-V

0.8

1.6

0.3

0.4

0.9

0.4

0.1

f2-M

0.6

1.6

1.7

0.2

0.8

0.4

0.1

f2-M

0.8

1.6

0.3

0.4

0.9

0.4

0.2

f2-B

0.6

1.6

1.7

0.2

0.8

0.4

0.1

f2-B

0.9

1.6

0.3

0.4

0.8

0.4

0.1

f1-V

0.7

1.6

1.7

0.2

0.9

0.4

0.1

f1-V

0.8

1.6

1.7

0.3

0.4

0.5

0.8

0.4

0.2

f1-M

0.8

1.7

1.8

0.3

0.9

0.3

0.1

f1-M

0.9

1.5

1.6

0.3

0.4

0.8

0.4

0.2

f1-B

0.7

1.7

1.8

0.2

0.3

0.2

0.9

0.4

0.1

f1-B

0.8

1.6

1.5

0.3

0.4

0.8

0.4

0.2

f0-V

0.7

1.6

1.8

0.2

0.9

0.4

0.1

f0-V

0.9

1.7

1.8

0.5

0.6

0.4

1.1

0.5

0.2

f0-M

0.8

1.8

1.9

0.3

0.2

0.9

0.4

0.1

f0-M

1.1

1.8

0.5

1.1

0.5

0.2

f0-B

0.7

1.6

1.9

0.3

0.2

1.0

0.4

0.1

f0-B

1.0

1.8

1.9

0.5

1.0

0.5

0.3

t(3)

t(5)

LN(.5)

LN(1)

LN(2)

G(.5)

G(2)

G(5)

t(3)

t(5)

LN(.5)

LN(1)

LN(2)

G(.5)

G(2)

G(5)

Crime

Garment

f2-V

37.6

122.5

129.0

19.8

19.5

17.2

67.9

31.0

10.1

f2-V

0.4

2.0

0.7

0.6

2.2

1.3

1.0

f2-M

36.3

117.7

123.2

18.4

17.7

13.9

57.9

27.9

9.0

f2-M

1.0

3.4

3.6

0.7

0.5

1.8

1.2

1.0

f2-B

40.5

132.8

129.5

20.1

18.3

15.1

60.6

29.9

9.3

f2-B

0.5

2.0

2.1

0.8

0.7

0.6

2.2

1.4

1.0

f1-V

12.0

30.5

31.1

9.3

9.5

9.6

18.5

14.3

8.4

f1-V

0.5

1.9

0.6

1.6

1.1

f1-M

11.5

29.4

30.1

8.4

7.8

7.5

15.9

11.8

7.7

f1-M

1.1

3.4

3.5

0.6

0.7

0.5

1.3

1.2

f1-B

12.0

30.5

31.2

8.9

9.6

9.5

18.2

14.1

8.4

f1-B

0.5

2.0

0.6

1.6

1.1

f0-V

10.1

29.3

30.7

6.4

6.3

5.5

13.9

10.5

7.1

f0-V

0.4

1.6

1.5

0.7

0.6

1.7

1.5

1.3

f0-M

10.0

28.2

29.9

6.0

5.7

5.1

13.8

10.2

6.4

f0-M

0.7

2.7

2.6

0.6

0.4

1.2

f0-B

10.2

29.2

30.9

6.7

6.4

5.6

14.3

11.4

7.6

f0-B

0.4

1.5

1.6

0.7

0.6

1.5

1.3

1.2

t(3)

t(5)

LN(.5)

LN(1)

LN(2)

G(.5)

G(2)

G(5)

t(3)

t(5)

LN(.5)

LN(1)

LN(2)

G(.5)

G(2)

G(5)

Ozone f2-V

0.6

1.8

1.9

0.3

1.0

0.5

0.2

f2-M

0.7

1.9

2.0

0.3

1.0

0.5

0.2

f2-B

0.6

1.8

1.9

0.3

1.0

0.5

0.2

f1-V

0.7

1.5

1.8

0.4

0.3

0.8

0.6

0.4

f1-M

0.6

1.6

1.8

0.4

1.0

0.6

0.4

f1-B

0.6

1.5

1.6

0.4

0.3

0.9

0.6

0.5

f0-V

0.4

1.0

1.1

0.2

0.5

0.4

0.3

f0-M

0.4

1.1

1.2

0.2

0.5

0.4

0.3

f0-B

0.3

1.0

1.2

0.2

0.6

0.4

0.3

t(3)

t(5)

LN(.5)

LN(1)

LN(2)

G(.5)

G(2)

G(5)

Time (min)

Figure 3: Average runtime (minutes) of each PRTree configuration.

analysis. Across datasets, computational cost is driven primarily by predictor dimension through its effect on the number of candidate splits that receive a full MSE evaluation. Within a given dataset, runtime is further influenced by the fill strategy and, to a lesser extent, by the choice of smoothing distribution. Although some of the most accurate configurations are also among the most computationally demanding, the relationship between runtime and predictive performance is far from monotonic. Consequently, computational cost should be viewed as a consequence of specific algorithmic design choices rather than as a direct indicator of predictive quality.

Conclusion

This paper introduced three different strategies to adapt probabilistic regression trees allowing it to handle missing predictor values directly during tree construction. Rather than relying on external imputation procedures or discarding incomplete observations, the proposed approach incorporates prob17

abilistic mechanisms into the split-selection process, allowing incomplete observations to contribute to model estimation throughout the recursive partitioning procedure. The empirical results demonstrate that PRTree can provide substantial predictive improvements over CART when missing predictor values affect a considerable fraction of the data. The magnitude of these improvements depends strongly on the adopted fill strategy. Across all datasets considered in this study, the fill mechanism emerged as the dominant modeling component, often producing larger performance differences than those associated with either the proxy-selection criterion or the choice of smoothing distribution. In contrast, the latter components primarily acted as refinement mechanisms, producing more localized improvements once an appropriate fill strategy had been selected. The experiments also revealed that the effectiveness of the probabilistic framework depends on the amount of missing values present in the data. Datasets with a large proportion of observations containing missing predictor values, such as Crime and Ozone, benefited substantially from the proposed methodology. Conversely, when missing values affected only a small fraction of the observations, as in Auto MPG, the differences between PRTree and CART became considerably smaller. These findings suggest that the potential gains obtained from probabilistic tree construction are closely related to the extent to which incomplete observations participate in the learning process. From a computational perspective, PRTree is substantially more demanding than CART. The additional cost arises from repeated tree reconstructions during smoothing-parameter selection and from the evaluation of multiple candidate splits using the full MSE criterion. The experiments indicate that computational cost is driven primarily by the interaction between sample size and the number of candidate splits selected for full evaluation. Nevertheless, the most computationally expensive configurations were not systematically associated with the best predictive performance, indicating that increased computational effort alone does not guarantee improved accuracy. Several directions for future research remain open. The current implementation considers a fixed family of proxy measures and smoothing distributions, but alternative probabilistic specifications may further improve predictive performance or computational efficiency. Extensions to ensemble methods, such as random forests, constitute a particularly promising direction, as the probabilistic framework proposed here is naturally compatible with bootstrap aggregation and randomized split selection. Additional work is also needed to investigate alternative missing-data mechanisms and to develop more computationally efficient procedures for candidate-split evaluation. Overall, the results indicate that probabilistic tree construction provides a viable and effective alternative for regression problems involving incomplete predictor information. By integrating missingvalue handling directly into the tree-building process, PRTree offers a flexible framework that preserves the interpretability of regression trees while improving robustness in the presence of missing covariate values.

Acknowledgments G. Pumi and T.S. Prass gratefully acknowledge the financial support received by the Conselho Nacional de Desenvolvimento Cientı́fico e Tecnológico – CNPq Brasil – Bolsa de Produtividade em Pesquisa Proc. 303281/2025-1 (Pumi) and 305886/2025-8 (Prass). A.S. Neimaier gratefully acknowledge the financial support received by the Coordenação de Aperfeiçoamento de Pessoal de Nı́vel Superior – Brasil (CAPES) – Finance Code 001.

References Alkhoury, S., Devijver, E., Clausel, M., Tami, M., Gaussier, E., Oppenheim, G., 2020. Smooth and consistent probabilistic regression trees, in: Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., Lin, H. (Eds.), Advances in Neural Information Processing Systems, Curran Associates, Inc.. pp. 11345–11355. URL: https://dl.acm.org/doi/pdf/10.5555/3495724.3496676. Balamuta, J.J., Truong, P., 2024. Explore UCI ML Repository Datasets. https://CRAN.R-project.org/package=ucimlrepo. r package version 0.0.2.

URL:

Breiman, L., Friedman, J., Stone, C.J., Olshen, R., 1984. Classification and Regression Trees. Chapman and Hall/CRC. da Rosa, J.C., Veiga, A., Medeiros, M.C., 2008. Tree-structured smooth transition regression models. Computational Statistics & Data Analysis 52, 2469–2488. Leisch, F., Dimitriadou, E., 2026. mlbench: Machine Learning Benchmark Problems. https://CRAN.R-project.org/package=mlbench. r package version 2.1-8.

URL:

Linero, A.R., Yang, Y., 2018. Bayesian regression tree ensembles that adapt to smoothness and sparsity. Journal of the Royal Statistical Society Series B: Statistical Methodology 80, 1087–1110. Molenberghs, G., Fitzmaurice, G., Kenward, M., Tsiatis, A., Verbeke, G., 2020. Handbook of Missing Data Methodology. Chapman & Hall/CRC Handbooks of Modern Statistical Methods, Taylor & Francis Group. Morvan, M.L., Josse, J., Scornet, E., Varoquaux, G., 2021. What’s a good imputation to predict with missing values?, in: Proceedings of the 35th International Conference on Neural Information Processing Systems, Curran Associates Inc.. pp. 11530–11540. Prass, T.S., Neimaier, A.S., 2026. Probabilistic Regression https://CRAN.R-project.org/package=PRTree. R package version 1.1.0.

Trees.

URL:

Prass, T.S., Neimaier, A.S., Pumi, G., 2026. Prtree: An R package for probabilistic regression trees with built-in missing data handling. URL: https://doi.org/10.48550/arXiv.2510.03634. Therneau, T., Atkinson, B., 2019. rpart: Recursive Partitioning and Regression Trees. https://CRAN.R-project.org/package=rpart. R package version 4.1-15.

URL:

İrsoy, O., Yıldız, O.T., Alpaydın, E., 2012. Soft decision trees, in: Proceedings of the 21st International Conference on Pattern Recognition (ICPR2012), pp. 1819–1822.

Handling Missing Data in Probabilistic Regression Trees

Analyze & explore related research in Conceptio

Related documents