# Trinary Decision Trees for missing value handling

Henning Zakrisson  Henning Zakrisson Department of Mathematics Stockholm
University<EMAIL_ADDRESS>
###### Abstract.

This paper introduces the Trinary decision tree, an algorithm designed to
improve the handling of missing data in decision tree regressors and
classifiers. Unlike other approaches, the Trinary decision tree does not
assume that missing values contain any information about the response. Both
theoretical calculations on estimator bias and numerical illustrations using
real data sets are presented to compare its performance with established
algorithms in different missing data scenarios (Missing Completely at Random
(MCAR), and Informative Missingness (IM)). Notably, the Trinary tree
outperforms its peers in MCAR settings, especially when data is only missing
out-of-sample, while lacking behind in IM settings. A hybrid model, the
TrinaryMIA tree, which combines the Trinary tree and the Missing In Attributes
(MIA) approach, shows robust performance in all types of missingness. Despite
the potential drawback of slower training speed, the Trinary tree offers a
promising and more accurate method of handling missing data in decision tree
algorithms. Keywords: Missing data, Decision trees, Regularization

##  1\. Introduction

Missing values are prevalent in real data. As noted by e.g. Nijman et al.
(2022), this is often not handled or mentioned in machine learning
applications in a satisfactory way. Classification and Regression Trees
(CART), as defined by Breiman et al. (1984) provide numerous ways to handle
missing values in covariates. Since CARTs are are the foundation of many
increasingly popular machine learning algorithms such as Gradient Boosting
Machines (GBMs) (Friedman, 2001), Random Forests (Ho, 1995), and XGBoost (Chen
et al., 2015), they are still relevant today. Yet, the proposed methods of
handling missing data come with drawbacks.

The simplest way to handle missing values when training a tree is to simply
ignore them by discarding data points with any missing feature. This of course
means risking losing potentially useful information, and is not an option when
predicting on data with missing values. The perhaps second-simplest method is
using the majority rule, where data points with missing values are assigned to
the daughter node with the largest amount of data in training. Another method
is presented by Twala et al. (2008), the Missing In Attributes (MIA)
algorithm. MIA splits data points with missing covariates in the way that
minimizes the loss for the training data. This is similar, but not identical,
to assigning missing values an own category in a categorical feature. Quinlan
(1993) introduces the C4.5 algorithm for decision trees, and with that a
weighted probabilistic strategy for missing value-handling, henceforth
referred to as Fractional Case - FC. In FC, a data point with a missing value
in a split is assigned a weight of membership in both daughter nodes depending
on the distribution of the observable data in the node. For out-of-sample
data, the weights for all terminal nodes are calculated and the prediction is
given as a weighted average. Breiman et al. (1984) proposes using so-called
surrogate splits in order to find other covariates on which the data points
which lack the relevant observation can be split to form similar splits. This
requires that there are no missing values in the surrogate covariate - or that
a secondary surrogate variable is found in its place.

Cons of thee methods include losing potentially useful information (discarding
data, FC), assuming there is always information in missingness (MIA),
requiring missing values in the training data to be able to handle missing
values in out-of-sample prediction (MIA, surrogate splits) or losing
interpretability (FC).

The trinary decision tree for missing value handling (henceforth Trinary tree)
introduced in this paper has four important attributes:

  * •

It does not assume that missing data points contain any information about the
response

  * •

It can handle missing values in predictions even if it was trained on a data
set with no missing data

  * •

It maintains the interpretability of a standard decision tree

  * •

It produces locally unbiased estimators of tree node values - which the other
algorithms do not necessarily do

The first three are apparent from the algorithm, which is presented in Section
2, whereas the fourth one is proven in Section 3. In Section 4, the algorithm
is tested against its peers with real data sets.

##  2\. The Trinary tree

Consider a loss function $\mathcal{L}((y_{i})_{i\in\mathcal{I}},\delta)$,
where $\mathcal{I}$ is an index set and $\delta$ is a parameter. For
regression problems, $\delta$ is a real number, and for classification
problems, $\delta$ is a probability vector. In this paper, two loss functions
will be considered: the sum of squared errors (SSE) for regression and the
point-wise cross-entropy for classification. These are defined as

(1)
$\mathcal{L}^{\textrm{SSE}}((y_{i})_{i\in\mathcal{I}},\delta)=\sum_{i\in\mathcal{I}}(y_{i}-\delta)^{2},\qquad\mathcal{L}^{\textrm{XE}}((y_{i})_{i\in\mathcal{I}},\delta)=\sum_{i\in\mathcal{I}}\log(\delta_{y_{i}})$

respectively. Let $\delta_{\mathcal{I}}$ denote the minimizer of
$\mathcal{L}((y_{i})_{i\in\mathcal{I}},\delta)$, i.e.

(2)
$\delta_{\mathcal{I}}=\operatorname*{arg\,min}_{\delta}\mathcal{L}((y_{i})_{i\in\mathcal{I}},\delta)$

For both SSE and point-wise cross-entropy, $\delta_{\mathcal{I}}$ has a closed
form solution, namely

(3)
$\delta_{\mathcal{I}}^{\textrm{SSE}}=\frac{1}{|\mathcal{I}|}\sum_{i\in\mathcal{I}}y_{i}\qquad\left(\delta_{\mathcal{I}}^{\textrm{XE}}\right)_{k}=\frac{1}{|\mathcal{I}|}\sum_{i\in\mathcal{I}}\mathbbm{1}_{\\{y_{i}\\}=k},\quad
k=1,\dots,K$

where $K$ is the number of classes in the classification problem. A binary
decision tree is generally fitted greedily by minimizing the loss function at
each so called node separately, starting with the root node containing all
data points. For data set $(y_{i},x_{i})_{i=1}^{n}$, where
$y_{i}\in\mathcal{Y}$ is the response and $x_{ij}\in\mathcal{X}_{j}$,
$j=1,\dots,p$, is a covariate, this means finding a greedy solution to

(4) $\min_{\begin{subarray}{c}j\in\\{1,\dots,p\\},\\\
\mathcal{X}_{j}^{l},\mathcal{X}_{j}^{r}\in\mathcal{X}_{j}\end{subarray}}\left\\{\mathcal{L}((y_{i})_{i\in\mathcal{I}_{jl}},\delta_{\mathcal{I}_{jl}})+\mathcal{L}((y_{i})_{i\in\mathcal{I}_{jr}},\delta_{\mathcal{I}_{jr}})\right\\},$

where

(5)
$\mathcal{I}_{jl}=\\{i\in\mathcal{I}:x_{ij}\in\mathcal{X}_{j}^{l}\\},\qquad\mathcal{I}_{jr}=\\{i\in\mathcal{I}:x_{ij}\in\mathcal{X}_{j}^{r}\\},$

such that $\mathcal{X}_{j}^{l}\cup\mathcal{X}_{j}^{r}=\mathcal{X}_{j}$. In the
case where $\mathcal{X}\in\mathbb{R}$, $\mathcal{X}_{j}^{l}$ and
$\mathcal{X}_{j}^{r}$ are constrained to be continuous intervals. After
finding the optimal split, the procedure is repeated for the so called
daughter nodes, i.e. the nodes that contain the data points from the two sides
of the split. The procedure is generally continued until reaching some
stopping criterion, such as a maximum tree depth. The nodes that are not split
are called terminal nodes.

Data points where the chosen splitting covariate is missing can be handled in
a number of ways. The majority strategy assigns them to the daughter node with
the largest amount of data (see Algorithm 2 in Appendix A). The MIA strategy
assigns them to the daughter node that provides the largest reduction in the
loss function (see Algorithm 3 in Appendix A). The FC strategy assigns them to
both daughter nodes with weights depending on the distribution of the
observable data in the node (see Algorithm 4 in Appendix A).

In contrast, the Trinary strategy assigns them to a third daughter node, which
changes the minimization problem in (4) to

(6) $\min_{\begin{subarray}{c}j\in\\{1,\dots,p\\},\\\
\mathcal{X}_{j}^{l},\mathcal{X}_{j}^{r}\in\mathcal{X}_{j}\end{subarray}}\left\\{\mathcal{L}((y_{i})_{i\in\mathcal{I}_{jl}},\delta_{\mathcal{I}_{jl}})+\mathcal{L}((y_{i})_{i\in\mathcal{I}_{jr}},\delta_{\mathcal{I}_{jr}})+\mathcal{L}((y_{i})_{i\in\mathcal{I}_{jm}},\delta_{\mathcal{I}})\right\\},$

where $\mathcal{I}_{jm}=\\{i\in\mathcal{I}:x_{ij}\textrm{ missing}\\}$. Note
that $\delta_{\mathcal{I}}$ in the third term is the minimizer of the loss
function over the entire data set of the node. This means that the third term
evaluates the loss of the points that are not assigned to the left or right
nodes as if it retained the parameter estimate of the mother node. After
finding a split, the procedure is repeated for all three daughter nodes. For
the third node, the entire data set is used for further splitting, but
omitting the splitting covariate. Thus, the third node will grow further by
first splitting on the second-best covariate, then continue to grow. The third
node is considered to be at the same depth level as the mother node, since the
data set has not been split.

The point of the third node is to avoid making assumptions about the missing
data. By not assigning the missing data to either of the standard daughter
nodes, the Trinary tree does not contaminate the $\delta$ estimates for the
standard nodes with data points that do not belong there. Instead, the Trinary
tree uses the entire data set to estimate $\delta$ for the third node, as a
way to regularize predictions when in uncertainty about important covariates.

The Trinary tree training algorithm is summarized in Algorithm 1. A
visualization of a Trinary tree with depth $1$ is shown in Figure 1.

Let • $(y_{i},x_{i})_{i\in\mathcal{I}}$, where $y_{i}\in\mathcal{Y}$,
$x_{ij}\in\mathcal{X}_{j}$, $j=1,\ldots,p$, be the training data •
$\mathcal{L}((y_{i})_{i\in\mathcal{I}},\delta)$ be the loss given parameter
$\delta$ • $\delta_{\mathcal{I}}$ be the parameter that minimizes
$\mathcal{L}$ for responses $(y_{i})_{i\in\mathcal{I}}$ •
$\mathcal{I}_{jl}=\\{i\in\mathcal{I}:x_{ij}\in\mathcal{X}_{j}^{l}\\}$,
$\mathcal{I}_{jr}=\\{i\in\mathcal{I}:x_{ij}\in\mathcal{X}_{j}^{r}\\}$, and
$\mathcal{I}_{jm}=\\{i:x_{ij}\textrm{ missing}\\}$ • $d_{\textrm{max}}$ be the
maximum depth • $n$ be the minimum number of samples per node Define training
function $\mathcal{T}\left((y_{i},x_{i})_{i\in\mathcal{I}},d\right)$
$\rightarrow$ $h(x)$: If $d=d_{\textrm{max}}$ or $|\mathcal{I}|=n$: Output
$\displaystyle h(x)=\delta_{\mathcal{I}}$ Else: Fit Find $j$,
$\mathcal{X}_{j}^{l}$, and $\mathcal{X}_{j}^{r}$ that solve
$\min_{j,\mathcal{X}_{j}^{l},\mathcal{X}_{j}^{r}}\left\\{\mathcal{L}\left((y_{i})_{i\in\mathcal{I}_{jl}},\delta_{\mathcal{I}_{jl}}\right)\\!+\\!\mathcal{L}\left((y_{i})_{i\in\mathcal{I}_{jr}},\delta_{\mathcal{I}_{jr}}\right)\\!+\\!\mathcal{L}\left((y_{i})_{i\in\mathcal{I}_{jm}},\delta_{\mathcal{I}}\right)\right\\}$
such that $\left|\mathcal{I}_{l}\right|\geq n$,
$\left|\mathcal{I}_{r}\right|\geq n$,
$\mathcal{X}_{j}^{l}\cup\mathcal{X}_{j}^{r}=\mathcal{X}_{j}$ Grow
$\displaystyle h_{l}(x)$
$\displaystyle=\mathcal{T}\left(\left(y_{i},x_{i}\right)_{\mathcal{I}_{l}},d+1\right)$
$\displaystyle h_{r}(x)$
$\displaystyle=\mathcal{T}\left(\left(y_{i},x_{i}\right)_{\mathcal{I}_{r}},d+1\right)$
$\displaystyle h_{m}(x)$
$\displaystyle=\mathcal{T}\left(\left(y_{i},x_{i}\right)_{\mathcal{I}},d\right)$
Output $\displaystyle h(x)=\begin{cases}h_{l}(x),&x_{\cdot
j}\in\mathcal{X}_{j}^{l}\\\ h_{r}(x),&x_{\cdot j}\in\mathcal{X}_{j}^{r}\\\
h_{m}(x),&x_{\cdot j}\textrm{ missing}\end{cases}$

Algorithm 1 Trinary tree training algorithm
$\delta_{0}$$\delta_{1}$$x_{i1}\in\mathcal{X}_{1}^{l}$$\delta_{2}$$x_{i1}\in\mathcal{X}_{1}^{r}$$\delta_{0}$$\delta_{3}$$x_{i2}\in\mathcal{X}_{2}^{l}$$\delta_{4}$$x_{i2}\in\mathcal{X}_{2}^{r}$$x_{i1}\textrm{
missing}$$\delta_{0}$$x_{i2}\textrm{ missing}$ Figure 1. Visualization of a
Trinary tree with depth $1$ for a covariate with $p=2$ dimensions. Note that
since the third node is considered to be at the same depth level as the root
node, an additional split is made. Since the best performing split covariate
$j=1$ is no longer available, the second-best split covariate $j=2$ is used
for the second split. Since that covariate could also be missing, the third
node has its own third daughter node. Since there are no further covariates to
split on, that node is a terminal node.

######  Remark 1.

Most implementations of decision trees use different tricks and approximations
to reduce the computational complexity when handling high-cardinality
categorical covariates. Two common ones are one-hot encoding and target
encoding. In this paper, we do not discuss any of these methods, but the
implementation in Section 4 does use the ordering trick described in Section
9.2.4 of Hastie et al. (2009) in cases of binary classification or regression.
Note that this is not an approximation but rather a faster way to find
appropriate splits. For details, see Hastie et al. (2009).

##  3\. Tree-fitting estimate bias

In order to illustrate how the Trinary tree might be preferable to the other
methods, let us examine a simple example where the non-trinary methods
estimators are locally biased. Consider the data set
$\mathcal{D}=(X_{i},Y_{i})_{i=1}^{n}$, where $X_{i}\in\mathcal{X}$ is the
covariate and $Y_{i}\in\mathbb{R}$ is the response. Let the expected value of
$Y$ follow a tree structure, i.e. let $\mathbb{E}[Y|X=x]=h(x)$ where

$h(x)=\begin{cases}a,&x\in\mathcal{X}^{l}\\\
b,&x\in\mathcal{X}^{r}\end{cases}$

where $\mathcal{X}^{l}\cup\mathcal{X}^{r}=\mathcal{X}$, $a<b$ and
$\mathbb{P}(X\in\mathcal{X}^{r})=p$, $0<p<1$. Then, consider a censored
dataset $\widetilde{\mathcal{D}}=(\tilde{X}_{i},Y_{i})_{i=1}^{n}$ where

$\tilde{X_{i}}=\begin{cases}\texttt{nan},&\textrm{with probability }q\\\
X_{i},&\textrm{with probability }1-q\\\ \end{cases}$

where nan corresponds to a missing value, $q>0$. Now, consider fitting trees
to this for an SSE loss function, using the Majority, MIA and FC algorithms
respectively, and examine their estimates of $a$. For brevity, let
$\mathcal{I}_{l}$ and $\mathcal{I}_{r}$ be defined as in (5) and introduce
index sets

$\mathcal{I}^{o}=\left\\{i:\tilde{X}_{i}=X_{i}\right\\},\qquad\mathcal{I}^{m}=\left\\{i:\tilde{X}_{i}=\texttt{nan}\right\\},$

as well as the intersections

$\displaystyle\mathcal{I}^{o}_{l}$
$\displaystyle=\mathcal{I}^{o}\cap\mathcal{I}_{l},$
$\displaystyle\mathcal{I}^{o}_{r}$
$\displaystyle=\mathcal{I}^{o}\cap\mathcal{I}_{r},$
$\displaystyle\mathcal{I}^{m}_{l}$
$\displaystyle=\mathcal{I}^{o}\cap\mathcal{I}_{l},$
$\displaystyle\mathcal{I}^{m}_{r}$
$\displaystyle=\mathcal{I}^{m}\cap\mathcal{I}_{r}.$

For the majority rule algorithm, the estimate of $a$ is

$\hat{a}_{\textrm{Maj}}=\begin{cases}\frac{1}{\left|\mathcal{I}_{l}^{o}\right|}\sum\limits_{i\in\mathcal{I}_{l}^{o}}Y_{i},&\textrm{if
}\left|\mathcal{I}_{l}^{o}\right|>\left|\mathcal{I}_{r}^{o}\right|,\\\
\frac{1}{\left|\mathcal{I}_{l}^{o}\right|+\left|\mathcal{I}^{m}\right|}\sum\limits_{i\in\mathcal{I}_{l}^{o}\cup\mathcal{I}^{m}}Y_{i},&\textrm{else.}\\\
\end{cases}$

First note that

$\mathbb{E}[\hat{a}_{\textrm{Maj}}|\left|\mathcal{I}_{l}^{o}\right|\leq\left|\mathcal{I}_{r}^{o}\right|]=a$

together with that

$\displaystyle\mathbb{E}\left[\hat{a}_{\textrm{Maj}}|\left|\mathcal{I}_{l}^{o}\right|>\left|\mathcal{I}_{r}^{o}\right|\right]$
$\displaystyle=\mathbb{E}\left[\frac{1}{\left|\mathcal{I}_{l}^{o}\right|+\left|\mathcal{I}^{m}_{l}\right|+\left|\mathcal{I}^{m}_{r}\right|}\left(\sum\limits_{i\in\mathcal{I}_{l}^{o}}Y_{i}+\sum\limits_{i\in\mathcal{I}_{l}^{m}}Y_{i}+\sum\limits_{i\in\mathcal{I}_{r}^{m}}Y_{i}\right)\right]$
$\displaystyle=\mathbb{E}\left[\frac{1}{\left|\mathcal{I}_{l}^{o}\right|+\left|\mathcal{I}^{m}_{l}\right|+\left|\mathcal{I}^{m}_{r}\right|}\left(\sum\limits_{\mathcal{I}\in\left\\{\mathcal{I}_{l}^{o},\mathcal{I}_{l}^{m},\mathcal{I}_{r}^{m}\right\\}}\sum\limits_{i\in\mathcal{I}}\mathbb{E}\left[Y_{i}|i\in\mathcal{I}\right]\right)\right]$
$\displaystyle=\mathbb{E}\left[\frac{1}{\left|\mathcal{I}_{l}^{o}\right|+\left|\mathcal{I}^{m}_{l}\right|+\left|\mathcal{I}^{m}_{r}\right|}\left(a\left|\mathcal{I}_{l}^{o}\right|+a\left|\mathcal{I}^{m}_{l}\right|+b\left|\mathcal{I}^{m}_{r}\right|\right)\right]$
$\displaystyle\geq\frac{1}{\mathbb{E}\left[\left|\mathcal{I}_{l}^{o}\right|+\left|\mathcal{I}^{m}_{l}\right|+\left|\mathcal{I}^{m}_{r}\right|\right]}\mathbb{E}\left[a\left|\mathcal{I}_{l}^{o}\right|+a\left|\mathcal{I}^{m}_{l}\right|+b\left|\mathcal{I}^{m}_{r}\right|\right]$
$\displaystyle=a+\frac{pq}{1-p+pq}\left(b-a\right)$

where the concavity of the function

$\displaystyle g(x,y,z)=\frac{x}{x+y+z}$

is used for the inequality. Then, introduce
$\kappa=\mathbb{P}\left(\left|\mathcal{I}_{l}^{o}\right|\leq\left|\mathcal{I}_{r}^{o}\right|\right)$,
and note that

$\displaystyle\mathbb{E}\left[\hat{a}_{\textrm{Maj}}\right]$
$\displaystyle=\kappa\mathbb{E}[\hat{a}_{\textrm{Maj}}|\left|\mathcal{I}_{l}^{o}\right|\leq\left|\mathcal{I}_{r}^{o}\right|]+\left(1-\kappa\right)\mathbb{E}[\hat{a}_{\textrm{Maj}}|\left|\mathcal{I}_{l}^{o}\right|>\left|\mathcal{I}_{r}^{o}\right|]$
$\displaystyle\geq\kappa
a+\left(1-\kappa\right)\left(a+\frac{pq}{1-p+pq}\left(b-a\right)\right)$
$\displaystyle=a+\frac{\left(1-\kappa\right)pq}{1-p+pq}\left(b-a\right)$
$\displaystyle>a$

The proof for the MIA strategy is identical with the only change that $\kappa$
then means the probability that the loss is lower if the missing values are
assigned to the left node. For the Fractional Case strategy the estimate of
parameter $a$ has expected value

$\displaystyle\mathbb{E}\left[\hat{a}_{\textrm{FC}}\right]$
$\displaystyle=\mathbb{E}\left[\frac{1}{\sum\limits_{i=1}^{n}w_{i}^{l}}\sum\limits_{i=1}^{n}w_{i}^{l}Y_{i}\right]$
$\displaystyle\geq\frac{1}{\mathbb{E}\left[\sum\limits_{i=1}^{n}w_{i}^{l}\right]}\mathbb{E}\left[\sum\limits_{i=1}^{n}w_{i}^{l}Y_{i}\right]$
$\displaystyle=a+pq\left(b-a\right)$ $\displaystyle>a.$

where it can also be shown that in order for
$\mathbb{E}[\hat{Y}_{\textrm{FC}}]=\mathbb{E}[Y]$ to hold it is required that
$\mathbb{E}[\hat{b}_{\textrm{FC}}]<b$. Finally, for the Trinary tree, see that

$\displaystyle\mathbb{E}\left[\hat{a}_{\textrm{Tri}}\right]=\mathbb{E}\left[\frac{1}{\left|\mathcal{I}^{o}_{l}\right|}\sum\limits_{i\in\mathcal{I}^{o}_{l}}Y_{i}\right]=a.$

##  4\. Numerical illustration

In order to illustrate the benefits of the Trinary tree, data with increasing
missingness is created for the data sets in Table 1. All the data sets are
available online. The data sets have been chosen in order to provide a wide
array of applications, and varying characteristics of the data. The
performance of the individual data sets will not be evaluated specifically in
this paper, but rather the performance of the algorithms on the data sets as a
whole. The data sets have been minimally pre-processed, e.g. by removing any
potential missing values. This is done in order to have full control over the
missingness. The tree depth is first tuned by performing $10$-fold cross
validation on the full data set, with a maximum tree depth set to $5$. For the
classification problems, the folds are stratified so that the relative class
frequency is equal in every fold. The dataset is then censored, i.e. values
are replaced with missing values, (with missingness $q\%$ ranging from $0\%$
to $90\%$) in three different ways.

  * MCAR

$q\%$ of the data is removed from all features in the training and test set,
completely at random.

  * MCARTest

$q\%$ of the data is removed from all features in the test set, completely at
random. The training set is the full data set.

  * IM

$q\%$ of the data is removed from all features in the training and test set.
For numerical features, the largest values are removed first. For categorical
ones, they are removed on a category-by-category basis.

Thereafter, the four different tree algorithms (Majority, MIA, FC, and
Trinary) are trained and evaluated on the $10$ folds. Additionally, a fifth
tree algorithm, denoted TrinaryMIA is evaluated. This is an amalgamation of
Trinary and MIA that in every node evaluates whether a MIA-style or Trinary
Tree style split reduces the training loss the most and then picks the one
that does.

type | Name | size | features
---|---|---|---
regression | AutoMPG | 392 | 8
Black Friday | 550,068 | 6
Cement | 1,030 | 9
Life Expectancy | 138 | 17
classification | Titanic | 712 | 7
Lymphography | 142 | 19
Boston Housing | 506 | 14
Wheat seeds | 199 | 8
Table 1. Data sets for the numerical illustration

First the total loss is calculated for the entire data set with no missing
values, referred to as the full loss. Then the total loss is calculated for
the data set with increasing missingness, and the excess loss (i.e. the loss
divided by the full loss) is calculated. The average excess loss over all data
sets is presented in Figure 2. Since the MIA and majority strategies are
identical in cases where there is no missing data in the training data, MIA is
omitted from the middle figure. The same applies to Trinary and TrinaryMIA
trees.

$0$$0.2$$0.4$$0.6$$0.8$$0$$0.5$$1$$1.5$$2$MissingnessExcess
lossMCARMajorityMIAFCTrinaryTrinaryMIA$0$$0.2$$0.4$$0.6$$0.8$$0$$0.5$$1$$1.5$$2$MissingnessMCARTest$0$$0.2$$0.4$$0.6$$0.8$$0$$0.5$$1$$1.5$$2$MissingnessIM
Figure 2. Average excess loss per missingness ratio for the tree algorithms in
different kinds of missingness on all data sets

As can be seen, the TrinaryMIA tree is the best performer in the MCAR case,
followed by the Trinary tree. For higher levels of missingness all algorithms
seem to perform almost equally bad. In the MCARTest case, the Trinary tree is
the best performer, followed by the FC tree. For the IM case, the TrinaryMIA
and MIA trees perform very similarly, and are the best performers. The
performance of MIA in this setting is expected, but it seems like the
TrinaryMIA tree is able to find the appropriate splits as well.

##  5\. Concluding remarks

It is clear both from the bias calculations and the numerical illustration
that the Trinary tree has benefits over its peers in MCAR settings.
Especially, the performance of the algorithm in the case where data is only
missing in the test set is noteworthy. It is however important to remember
that assuming that missingness contains no information is an assumption in
itself - seen by the less impressive performance in the IM test. This drawback
seems to be easily overcome by the TrinaryMIA tree, which maintains
performance in all types of missingness. Surprisingly, the TrinaryMIA tree
also outperforms the Trinary tree in the MCAR case, however not by a large
margin.

The potential of using the Trinary tree algorithm as a weak learner in more
powerful machine learning algorithms, such as a GBM, is an interesting topic,
since the reguralization and missing value-handling would then be inherited by
the ensemble model.

The main drawback for the Trinary tree, if controlling for informative
missingness by using the TrinaryMIA method, is that for large data sets
(especially covariate data sets with many features and categorical features
with high cardinality) or deep trees, the training speed can suffer. For
shallow trees and data sets with a limited number of covariates, the speed is
however on par with the other methods. The implementation used in the
numerical illustration here is not optimized for speed, and since node
training can be parallelized, there is further speed potential if computation
time is of importance. It should also be noted that the TrinaryMIA test is
often slightly faster since it, if there is information in missingness in the
training data, will grow fewer nodes than a standard Trinary tree.

## Acknowledgements

During the preparation of this work the author used ChatGPT and GitHub Copilot
in order to improve language and readability. After using these tools, the
author reviewed and edited the content as needed and takes full responsibility
for the content of the publication.

## References

  * (1)
  * Breiman et al. (1984) Breiman, L., Friedman, J., Stone, C. & Olshen, R. (1984), Classification and Regression Trees, Taylor & Francis.
  * Chen et al. (2015) Chen, T., He, T., Benesty, M., Khotilovich, V., Tang, Y., Cho, H., Chen, K., Mitchell, R., Cano, I., Zhou, T. et al. (2015), ‘Xgboost: extreme gradient boosting’, R package version 0.4-2 1(4), 1–4.
  * Friedman (2001) Friedman, J. H. (2001), ‘Greedy function approximation: a gradient boosting machine’, Annals of statistics 25(5), 1189–1232.
  * Hastie et al. (2009) Hastie, T., Tibshirani, R., Friedman, J. H. & Friedman, J. H. (2009), The elements of statistical learning: data mining, inference, and prediction, Vol. 2, Springer.
  * Ho (1995) Ho, T. K. (1995), Random decision forests, in ‘Proceedings of 3rd international conference on document analysis and recognition’, Vol. 1, IEEE, pp. 278–282.
  * Nijman et al. (2022) Nijman, S., Leeuwenberg, A., Beekers, I., Verkouter, I., Jacobs, J., Bots, M., Asselbergs, F., Moons, K. & Debray, T. (2022), ‘Missing data is poorly handled and reported in prediction model studies using machine learning: a literature review’, Journal of clinical epidemiology 142, 218–229.
  * Quinlan (1993) Quinlan, J. (1993), C4.5: Programs for Machine Learning, Morgan Kaufmann series in machine learning, Elsevier Science.
  * Twala et al. (2008) Twala, B. E., Jones, M. & Hand, D. J. (2008), ‘Good methods for coping with missing data in decision trees’, Pattern Recognition Letters 29(7), 950–956.

##  Appendix A Tree algorithms

Let • $(y_{i},x_{i})_{i\in\mathcal{I}}$, where $y_{i}\in\mathcal{Y}$,
$x_{ij}\in\mathcal{X}_{j}$, $j=1,\ldots,p$, be the training data •
$\mathcal{L}((y_{i})_{i\in\mathcal{I}},\delta)$ be the loss given parameter
$\delta$ • $\delta_{\mathcal{I}}$ be the parameter that minimizes
$\mathcal{L}$ for responses $(y_{i})_{i\in\mathcal{I}}$ •
$\mathcal{I}_{jl}=\\{i\in\mathcal{I}:x_{ij}\in\mathcal{X}_{j}^{l}\\}$,
$\mathcal{I}_{jr}=\\{i\in\mathcal{I}:x_{ij}\in\mathcal{X}_{j}^{r}\\}$, and
$\mathcal{I}_{jm}=\\{i:x_{ij}\textrm{ missing}\\}$ • $d_{\textrm{max}}$ be the
maximum depth • $n$ be the minimum number of samples per node Define training
function $\mathcal{T}\left((y_{i},x_{i})_{i\in\mathcal{I}},d\right)$
$\rightarrow$ $h(x)$: If $d=d_{\textrm{max}}$ or $|\mathcal{I}|=n$: Output
$\displaystyle h(x)=\delta_{\mathcal{I}}$ Else: Fit Find $j$,
$\mathcal{X}_{j}^{l}$, and $\mathcal{X}_{j}^{r}$ that solve
$\min_{j,\mathcal{X}_{j}^{l},\mathcal{X}_{j}^{r}}\left\\{\mathcal{L}\left((y_{i})_{i\in\mathcal{I}_{jl}},\delta_{\widetilde{\mathcal{I}}_{l}}\right)+\mathcal{L}\left((y_{i})_{i\in\mathcal{I}_{jr}},\delta_{\widetilde{\mathcal{I}}_{r}}\right)\right\\}$
where
$\widetilde{\mathcal{I}}_{l}=\begin{cases}\mathcal{I}_{jl},&\left|\mathcal{I}_{jl}\right|\leq\left|\mathcal{I}_{jr}\right|\\\
\mathcal{I}_{jl}\cup\mathcal{I}_{jm},&\left|\mathcal{I}_{jl}\right|>\left|\mathcal{I}_{jr}\right|\end{cases}\qquad\widetilde{\mathcal{I}}_{r}=\begin{cases}\mathcal{I}_{jr}\cup\mathcal{I}_{jm},&\left|\mathcal{I}_{jl}\right|\leq\left|\mathcal{I}_{jr}\right|\\\
\mathcal{I}_{jr},&\left|\mathcal{I}_{jl}\right|>\left|\mathcal{I}_{jr}\right|\end{cases}$
and vice versa, such that $\left|\widetilde{\mathcal{I}}_{l}\right|\geq n$,
$\left|\widetilde{\mathcal{I}}_{r}\right|\geq n$,
$\mathcal{X}_{j}^{l}\cup\mathcal{X}_{j}^{r}=\mathcal{X}_{j}$ Grow
$\displaystyle h_{l}(x)$
$\displaystyle=\mathcal{T}\left(\left(y_{i},x_{i}\right)_{\widetilde{\mathcal{I}}_{l}},d+1\right)$
$\displaystyle h_{r}(x)$
$\displaystyle=\mathcal{T}\left(\left(y_{i},x_{i}\right)_{\widetilde{\mathcal{I}}_{r}},d+1\right)$
Output $\displaystyle h(x)=\begin{cases}h_{l}(x),&x_{\cdot
j}\in\mathcal{X}_{j}^{l}\textrm{ or }x_{\cdot j}\textrm{ missing and
}\left|\mathcal{I}_{jl}\right|\geq\left|\mathcal{I}_{jr}\right|\\\
h_{r}(x),&x_{\cdot j}\in\mathcal{X}_{j}^{r}\textrm{ or }x_{\cdot j}\textrm{
missing and
}\left|\mathcal{I}_{jl}\right|<\left|\mathcal{I}_{jr}\right|\end{cases}$

Algorithm 2 Majority Rule training algorithm

Let • $(y_{i},x_{i})_{i\in\mathcal{I}}$, where $y_{i}\in\mathcal{Y}$,
$x_{ij}\in\mathcal{X}_{j}$, $j=1,\ldots,p$, be the training data •
$\mathcal{L}((y_{i})_{i\in\mathcal{I}},\delta)$ be the loss given parameter
$\delta$ • $\delta_{\mathcal{I}}$ be the parameter that minimizes
$\mathcal{L}$ for responses $(y_{i})_{i\in\mathcal{I}}$ •
$\mathcal{I}_{jl}=\\{i\in\mathcal{I}:x_{ij}\in\mathcal{X}_{j}^{l}\\}$,
$\mathcal{I}_{jr}=\\{i\in\mathcal{I}:x_{ij}\in\mathcal{X}_{j}^{r}\\}$ •
$d_{\textrm{max}}$ be the maximum depth • $n$ be the minimum number of samples
per node Define training function
$\mathcal{T}\left((y_{i},x_{i})_{i\in\mathcal{I}},d\right)$ $\rightarrow$
$h(x)$: If $d=d_{\textrm{max}}$ or $|\mathcal{I}|=n$: Output $\displaystyle
h(x)=\delta_{\mathcal{I}}$ Else: Fit Find $j$, $\mathcal{X}_{j}^{l}$, and
$\mathcal{X}_{j}^{r}$ that solve
$\min_{j,\mathcal{X}_{j}^{l},\mathcal{X}_{j}^{r}}\left\\{\mathcal{L}\left((y_{i})_{i\in\mathcal{I}_{jl}},\delta_{\mathcal{I}_{jl}}\right)+\mathcal{L}\left((y_{i})_{i\in\mathcal{I}_{jr}},\delta_{\mathcal{I}_{jr}}\right)\right\\}$
such that $\left|\mathcal{I}_{jl}\right|\geq n$,
$\left|\mathcal{I}_{jr}\right|\geq n$,
$\mathcal{X}_{j}^{l}\cup\mathcal{X}_{j}^{r}=\mathcal{X}_{j}\cup\\{\textrm{nan}\\}$
Grow $\displaystyle h_{l}(x)$
$\displaystyle=\mathcal{T}\left(\left(y_{i},x_{i}\right)_{\mathcal{I}_{jl}},d+1\right)$
$\displaystyle h_{r}(x)$
$\displaystyle=\mathcal{T}\left(\left(y_{i},x_{i}\right)_{\mathcal{I}_{jr}},d+1\right)$
Output $\displaystyle h(x)=\begin{cases}h_{l}(x),&x_{\cdot
j}\in\mathcal{X}_{j}^{l}\\\ h_{r}(x),&x_{\cdot
j}\in\mathcal{X}_{j}^{r}\end{cases}$

Algorithm 3 Missing In Attributes training algorithm

Let • $(y_{i},x_{i},w_{i})_{i\in\mathcal{I}}$, where $y_{i}\in\mathcal{Y}$,
$x_{ij}\in\mathcal{X}_{j}$, $j=1,\ldots,p$, $w_{i}\in[0,1]$ be the training
data • $\mathcal{L}((y_{i},w_{i})_{i\in\mathcal{I}},\delta)$ be the loss given
parameter $\delta$ • $\delta_{\mathcal{I}}$ be the parameter that minimizes
$\mathcal{L}$ for responses and weights $(y_{i},w_{i})_{i\in\mathcal{I}}$ •
$\mathcal{I}_{jm}=\\{i\in\mathcal{I}:x_{ij}\textrm{ missing}\\}$,
$\mathcal{I}_{jl}=\\{i\in\mathcal{I}:x_{ij}\in\mathcal{X}_{j}^{l}\\}\cup\mathcal{I}_{jm}$,
$\mathcal{I}_{jr}=\\{i:x_{ij}\in\mathcal{X}_{j}^{r}\\}\cup\mathcal{I}_{jm}$ •
$d_{\textrm{max}}$ be the maximum depth • $n$ be the minimum total sample
weight per node Define training function
$\mathcal{T}\left((y_{i},x_{i})_{i\in\mathcal{I}},d\right)$ $\rightarrow$
$h(x)$: If $d=d_{\textrm{max}}$ or $|\mathcal{I}|\leq n$: Output
$\displaystyle h(x)=\delta_{\mathcal{I}}$ Else: Fit Find $j$,
$\mathcal{X}_{j}^{l}$, and $\mathcal{X}_{j}^{r}$ that solve
$\min_{j,\mathcal{X}_{j}^{l},\mathcal{X}_{j}^{r}}\left\\{\mathcal{L}\left((y_{i},w_{i}^{l})_{i\in\mathcal{I}_{jl}},\delta_{\mathcal{I}_{jl}}\right)+\mathcal{L}\left((y_{i},w_{i}^{r})_{i\in\mathcal{I}_{jr}},\delta_{\mathcal{I}_{jr}}\right)\right\\}$
where $w_{i}^{l}=\begin{cases}w_{i},&x_{ij}\in\mathcal{X}_{j}^{l}\\\
w_{i}\frac{\left|\mathcal{I}_{jl}\setminus\mathcal{I}_{jm}\right|}{\left|\mathcal{I}\setminus\mathcal{I}_{jm}\right|},&x_{ij}\in\\{\texttt{nan}\\}\end{cases}$
such that $\sum\limits_{i\in\mathcal{I}_{jl}}w_{i}^{l}\geq n$, and vice versa
for $r$, and $\mathcal{X}_{j}^{l}\cup\mathcal{X}_{j}^{r}=\mathcal{X}_{j}$ Grow
$\displaystyle h_{l}(x)$
$\displaystyle=\mathcal{T}\left(\left(y_{i},x_{i},w_{i}^{l}\right)_{\mathcal{I}_{jl}},d+1\right)$
$\displaystyle h_{r}(x)$
$\displaystyle=\mathcal{T}\left(\left(y_{i},x_{i},w_{i}^{r}\right)_{\mathcal{I}_{jr}},d+1\right)$
Output $\displaystyle h(x)=\begin{cases}h_{l}(x),&x_{\cdot
j}\in\mathcal{X}_{j}^{l}\\\ h_{r}(x),&x_{\cdot j}\in\mathcal{X}_{j}^{r}\\\
w_{l}h_{l}(x)+w_{r}h_{r}(x),&x_{\cdot j}\textrm{ missing}\end{cases}$

Algorithm 4 Fractional Case training algorithm