We have {$\mathbf{X}' = \mathbf{x}$} = {$\mathbf{X}' = \mathbf{x}$, $M_1 = 0$} = {$\mathbf{X} = \mathbf{x}$}, and consequently,

$$ \mathbb{E}[Y | \mathbf{X}' = \mathbf{x}] = \mathbb{E}[Y | \mathbf{X} = \mathbf{x}]. \quad (14) $$

Second Subcase: assume $\mathbb{P}[M_1 = 1 | X_2 = x_2, \dots, X_d = x_d] > 0$.

We have

$$ \begin{aligned} \mathbb{P}[\mathbf{X}' \in B(\mathbf{x}, \rho)] &= \mathbb{E}[\mathbf{1}_{\mathbf{X}' \in B(\mathbf{x}, \rho)} \mathbf{1}_{M_1=0}] + \mathbb{E}[\mathbf{1}_{\mathbf{X}' \in B(\mathbf{x}, \rho)} \mathbf{1}_{M_1=1}] \\ &= \mathbb{E}[\mathbf{1}_{\mathbf{X} \in B(\mathbf{x}, \rho)} \mathbf{1}_{M_1=0}] + \mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbf{1}_{M_1=1}], \end{aligned} $$

and

$$ \mathbb{E}[f(\mathbf{X})\mathbf{1}_{\mathbf{X}' \in B(\mathbf{x}, \rho)}] = E[f(\mathbf{X})\mathbf{1}_{\mathbf{X} \in B(\mathbf{x}, \rho)}\mathbf{1}_{M_1=0}] \\ + E[f(\mathbf{X})\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)}\mathbf{1}_{M_1=1}]. $$

Therefore,

$$ \begin{aligned} \mathbb{E}[Y | \mathbf{X}' \in B(\mathbf{x}, \rho)] &= \frac{\mathbb{E}[f(\mathbf{X})\mathbf{1}_{\mathbf{X}' \in B(\mathbf{x}, \rho)}]}{\mathbb{P}[\mathbf{X}' \in B(\mathbf{x}, \rho)]} \\ &= \frac{\mathbb{E}[f(\mathbf{X})\mathbf{1}_{\mathbf{X} \in B(\mathbf{x}, \rho)}\mathbf{1}_{M_1=0}] + \mathbb{E}[f(\mathbf{X})\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)}\mathbf{1}_{M_1=1}]}{\mathbb{E}[\mathbf{1}_{\mathbf{X} \in B(\mathbf{x}, \rho)}\mathbf{1}_{M_1=0}] + \mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)}\mathbf{1}_{M_1=1}]} \end{aligned} \quad (15) $$

The terms in (15) involving $M_1 = 0$ satisfy

$$ \mathbb{E}[\mathbf{1}_{\mathbf{X} \in B(\mathbf{x}, \rho)} \mathbf{1}_{M_1=0}] \le \mu(B(\mathbf{x}, \rho)) \le \frac{\pi^{d/2}}{\Gamma(\frac{d}{2} + 1)} \|g\|_{\infty} \rho^d, \quad (16) $$

and

$$ \begin{aligned} |\mathbb{E}[f(\mathbf{X})\mathbf{1}_{\mathbf{X} \in B(\mathbf{x}, \rho)}\mathbf{1}_{M_1=0}|| & \le \mathbb{E}[|f(\mathbf{X})|\mathbf{1}_{\mathbf{X} \in B(\mathbf{x}, \rho)}|] \\ & \le \frac{\pi^{d/2}}{\Gamma(\frac{d}{2} + 1)} \|g\|_{\infty} \|f\|_{\infty} \rho^d. \end{aligned} \quad (17) $$

The second term of the denominator in (15) can be bounded from below,

$$ \begin{aligned} & \mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbf{1}_{M_1=1}] \\ &= \mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbb{P}[M_1 = 1 | X_2, \dots, X_d]] \\ & \ge \frac{\pi^{(d-1)/2}}{\Gamma(\frac{d-1}{2} + 1)} (\inf_{[0,1]^d} g) \rho^{d-1} \eta. \end{aligned} \quad (18) $$
We have {$\mathbf{X}' = \mathbf{x}$} = {$\mathbf{X}' = \mathbf{x}$, $M_1 = 0$} = {$\mathbf{X} = \mathbf{x}$}, and consequently,

$$ \mathbb{E}[Y | \mathbf{X}' = \mathbf{x}] = \mathbb{E}[Y | \mathbf{X} = \mathbf{x}]. \quad (14) $$

Second Subcase: assume $\mathbb{P}[M_1 = 1 | X_2 = x_2, \dots, X_d = x_d] > 0$.

We have

$$ \begin{aligned} \mathbb{P}[\mathbf{X}' \in B(\mathbf{x}, \rho)] &= \mathbb{E}[\mathbf{1}_{\mathbf{X}' \in B(\mathbf{x}, \rho)} \mathbf{1}_{M_1=0}] + \mathbb{E}[\mathbf{1}_{\mathbf{X}' \in B(\mathbf{x}, \rho)} \mathbf{1}_{M_1=1}] \\ &= \mathbb{E}[\mathbf{1}_{\mathbf{X} \in B(\mathbf{x}, \rho)} \mathbf{1}_{M_1=0}] + \mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbf{1}_{M_1=1}], \end{aligned} $$

and

$$ \mathbb{E}[f(\mathbf{X})\mathbf{1}_{\mathbf{X}' \in B(\mathbf{x}, \rho)}] = E[f(\mathbf{X})\mathbf{1}_{\mathbf{X} \in B(\mathbf{x}, \rho)}\mathbf{1}_{M_1=0}] \\ + E[f(\mathbf{X})\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)}\mathbf{1}_{M_1=1}]. $$

Therefore,

$$ \begin{aligned} \mathbb{E}[Y | \mathbf{X}' \in B(\mathbf{x}, \rho)] &= \frac{\mathbb{E}[f(\mathbf{X})\mathbf{1}_{\mathbf{X}' \in B(\mathbf{x}, \rho)}]}{\mathbb{P}[\mathbf{X}' \in B(\mathbf{x}, \rho)]} \\ &= \frac{\mathbb{E}[f(\mathbf{X})\mathbf{1}_{\mathbf{X} \in B(\mathbf{x}, \rho)}\mathbf{1}_{M_1=0}] + \mathbb{E}[f(\mathbf{X})\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)}\mathbf{1}_{M_1=1}]}{\mathbb{E}[\mathbf{1}_{\mathbf{X} \in B(\mathbf{x}, \rho)}\mathbf{1}_{M_1=0}] + \mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)}\mathbf{1}_{M_1=1}]} \end{aligned} \quad (15) $$

The terms in (15) involving $M_1 = 0$ satisfy

$$ \mathbb{E}[\mathbf{1}_{\mathbf{X} \in B(\mathbf{x}, \rho)} \mathbf{1}_{M_1=0}] \le \mu(B(\mathbf{x}, \rho)) \le \frac{\pi^{d/2}}{\Gamma(\frac{d}{2} + 1)} \|g\|_{\infty} \rho^d, \quad (16) $$

and

$$ \begin{aligned} |\mathbb{E}[f(\mathbf{X})\mathbf{1}_{\mathbf{X} \in B(\mathbf{x}, \rho)}\mathbf{1}_{M_1=0}]| &\le \mathbb{E}[|f(\mathbf{X})|\mathbf{1}_{\mathbf{X} \in B(\mathbf{x}, \rho)}]| \\ &\le \frac{\pi^{d/2}}{\Gamma(\frac{d}{2} + 1)} \|g\|_{\infty} \|f\|_{\infty} \rho^d. \end{aligned} \quad (17) $$

The second term of the denominator in (15) can be bounded from below,

$$ \begin{aligned} &\mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbf{1}_{M_1=1}] \\ &= \mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbb{P}[M_1 = 1 | X_2, \dots, X_d]] \\ &\ge \frac{\pi^{(d-1)/2}}{\Gamma(\frac{d-1}{2} + 1)} (\inf_{[0,1]^d} g) \rho^{d-1} \eta. \end{aligned} \quad (18) $$
The second term of the numerator in (15) verifies

$$
\begin{align*}
& \mathbb{E}[f(\mathbf{X})\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbf{1}_{M_1=1}] \\
&= \mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbb{E}[f(\mathbf{X})\mathbf{1}_{M_1=1} | X_2, \dots, X_d]] \\
&= \mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbb{E}[f(\mathbf{X}) | X_2, \dots, X_d] \mathbb{E}[\mathbf{1}_{M_1=1} | X_2, \dots, X_d]].
\end{align*}
$$

If $\mathbb{E}[f(\mathbf{X})|X_2 = x_2, \dots, X_d = x_d] > 0$, by uniform continuity of $f$ and $g$,

$$
\begin{align*}
& \mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbb{E}[f(\mathbf{X})|X_2, \dots, X_d] \mathbb{E}[\mathbf{1}_{M_1=1}|X_2, \dots, X_d]] \\
& \geq \mathbb{E}[f(\mathbf{X})|X_2 = x_2, \dots, X_d = x_d] \frac{\pi^{(d-1)/2}}{\Gamma(\frac{d-1}{2} + 1)} (\inf_{[0,1]^d} g) \rho^{d-1} \eta.
\end{align*}
$$

Similarly, if $\mathbb{E}[f(\mathbf{X})|X_2 = x_2, \dots, X_d = x_d] < 0$, we have

$$
\begin{align*}
& \mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbb{E}[f(\mathbf{X})|X_2, \dots, X_d] \mathbb{E}[\mathbf{1}_{M_1=1}|X_2, \dots, X_d]] \\
& \leq \mathbb{E}[f(\mathbf{X})|X_2 = x_2, \dots, X_d = x_d] \frac{\pi^{(d-1)/2}}{\Gamma(\frac{d-1}{2} + 1)} (\inf_{[0,1]^d} g) \rho^{d-1} \eta \\
& \leq 0.
\end{align*}
$$

Hence, if $\mathbb{E}[f(\mathbf{X})|X_2 = x_2, \dots, X_d = x_d] \neq 0$

$$
\begin{align}
& |\mathbb{E}[f(\mathbf{X})\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbf{1}_{M_1=1}]| \nonumber \\
& \geq |\mathbb{E}[f(\mathbf{X})|X_2 = x_2, \dots, X_d = x_d]| \frac{\pi^{(d-1)/2}}{\Gamma(\frac{d-1}{2} + 1)} (\inf_{[0,1]^d} g) \rho^{d-1} \eta. \tag{19}
\end{align}
$$

Gathering inequalities (16)-(19) and using equation (15), we have, if $\mathbb{E}[f(\mathbf{X})|X_2 = x_2, \dots, X_d = x_d] \neq 0$

$$
\begin{align*}
\lim_{\rho \to 0} \mathbb{E}[Y | \mathbf{X}' \in B(\mathbf{x}, \rho)] &= \lim_{\rho \to 0} \frac{\mathbb{E}[f(\mathbf{X}) \mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbf{1}_{M_1=1}]}{\mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbf{1}_{M_1=1}]} \\
&= \mathbb{E}[f(\mathbf{X}) | X_2 = x_2, \dots, X_d = x_d, M_1=1].
\end{align*}
$$

Finally, if $\mathbb{E}[f(\mathbf{X})|X_2 = x_2, \dots, X_d = x_d] = 0$ then by uniform continuity of $f$, there exists $\varepsilon_\rho$
such that $\varepsilon_\rho \to 0$ as $\rho \to 0$ satisfying,

$$
|\mathbb{E}[f(\mathbf{X})\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbf{1}_{M_1=1}]| \le \varepsilon_\rho \rho^{d-1} \frac{\pi^{d/2}}{\Gamma(\frac{d}{2}+1)} \|g\|_{\infty},
$$
The second term of the numerator in (15) verifies

$$
\begin{align*}
& \mathbb{E}[f(\mathbf{X})\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbf{1}_{M_1=1}] \\
&= \mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbb{E}[f(\mathbf{X})\mathbf{1}_{M_1=1}|X_2, \dots, X_d]] \\
&= \mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbb{E}[f(\mathbf{X})|X_2, \dots, X_d] \mathbb{E}[\mathbf{1}_{M_1=1}|X_2, \dots, X_d]].
\end{align*}
$$

If $\mathbb{E}[f(\mathbf{X})|X_2 = x_2, \dots, X_d = x_d] > 0$, by uniform continuity of $f$ and $g$,

$$
\begin{align*}
& \mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbb{E}[f(\mathbf{X})|X_2, \dots, X_d] \mathbb{E}[\mathbf{1}_{M_1=1}|X_2, \dots, X_d]] \\
& \geq \mathbb{E}[f(\mathbf{X})|X_2 = x_2, \dots, X_d = x_d] \frac{\pi^{(d-1)/2}}{\Gamma(\frac{d-1}{2} + 1)} (\inf_{[0,1]^d} g) \rho^{d-1} \eta.
\end{align*}
$$

Similarly, if $\mathbb{E}[f(\mathbf{X})|X_2 = x_2, \dots, X_d = x_d] < 0$, we have

$$
\begin{align*}
& \mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbb{E}[f(\mathbf{X})|X_2, \dots, X_d] \mathbb{E}[\mathbf{1}_{M_1=1}|X_2, \dots, X_d]] \\
& \leq \mathbb{E}[f(\mathbf{X})|X_2 = x_2, \dots, X_d = x_d] \frac{\pi^{(d-1)/2}}{\Gamma(\frac{d-1}{2} + 1)} (\inf_{[0,1]^d} g) \rho^{d-1} \eta \\
& \leq 0.
\end{align*}
$$

Hence, if $\mathbb{E}[f(\mathbf{X})|X_2 = x_2, \dots, X_d = x_d] \neq 0$

$$
\begin{align}
& |\mathbb{E}[f(\mathbf{X})\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbf{1}_{M_1=1}]| \nonumber \\
& \geq |\mathbb{E}[f(\mathbf{X})|X_2 = x_2, \dots, X_d = x_d]| \frac{\pi^{(d-1)/2}}{\Gamma(\frac{d-1}{2} + 1)} (\inf_{[0,1]^d} g) \rho^{d-1} \eta. \tag{19}
\end{align}
$$

Gathering inequalities (16)-(19) and using equation (15), we have, if $\mathbb{E}[f(\mathbf{X})|X_2 = x_2, \dots, X_d = x_d] \neq 0$

$$
\begin{align*}
\lim_{\rho \to 0} \mathbb{E}[Y | \mathbf{X}' \in B(\mathbf{x}, \rho)] &= \lim_{\rho \to 0} \frac{\mathbb{E}[f(\mathbf{X}) \mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbf{1}_{M_1=1}]}{\mathbb{E}[\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)} \mathbf{1}_{M_1=1}]} \\
&= \mathbb{E}[f(\mathbf{X}) | X_2 = x_2, \dots, X_d = x_d, M_1=1].
\end{align*}
$$

Finally, if $\mathbb{E}[f(\mathbf{X})|X_2 = x_2, \dots, X_d = x_d] = 0$ then by uniform continuity of $f$, there exists $\epsilon_\rho$
such that $\epsilon_\rho \to 0$ as $\rho \to 0$ satisfying,

$$
|\mathbb{E}[f(\mathbf{X})\mathbf{1}_{(X_2, \dots, X_d) \in B((x_2, \dots, x_d), \rho)}\mathbf{1}_{M_1=1}|] \leq \epsilon_\rho \rho^{d-1} \frac{\pi^{d/2}}{\Gamma(\frac{d}{2} + 1)} \|g\|_{\infty},
$$
hence

$$
\begin{align*}
\lim_{\rho \to 0} \mathbb{E}[Y | \mathbf{X}' \in B(\mathbf{x}, \rho)] &= 0 \\
&= \mathbb{E}[f(\mathbf{X}) | X_2 = x_2, \dots, X_d = x_d] \\
&= \mathbb{E}[f(\mathbf{X}) | X_2 = x_2, \dots, X_d = x_d, M_1 = 1],
\end{align*}
$$

since $M_1 \perp X_1|(X_2, \dots, X_d)$. Consequently, for all $\mathbf{x} \in [0, 1]^d$ such that $x_1 = \alpha$,

$$
\lim_{\rho \to 0} \mathbb{E}[Y | \mathbf{X}' \in B(\mathbf{x}, \rho)] = \mathbb{E}[f(\mathbf{X}) | X_2 = x_2, \dots, X_d = x_d, M_1 = 1]. \quad (20)
$$

Combining equations (13), (14) and (20), the prediction given by the mean imputation followed
by learning is, for all $\mathbf{x}' \in \mathbb{R}^d$,

$$
\begin{align*}
f_{\mathbf{X}'}^{*}(\mathbf{x}') ={}& \mathbb{E}[Y | X_2 = x_2, \dots, X_d = x_d, M_1 = 1] \mathbf{1}_{x_1' = \alpha} \mathbf{1}_{\mathbb{P}[M_1=1|X_2=x_2, \dots, X_d=x_d]>0} \\
& + \mathbb{E}[Y | \mathbf{X} = \mathbf{x}'] \mathbf{1}_{x_1' = \alpha} \mathbf{1}_{\mathbb{P}[M_1=1|X_2=x_2, \dots, X_d=x_d]=0} \\
& + \mathbb{E}[Y | X_2 = x_2, \dots, X_d = x_d, M_1 = 0] \mathbf{1}_{x_1' \neq \alpha},
\end{align*}
$$

which concludes the proof.

B Decision trees: an example of empirical risk minimization with missing
data

Aside from using procedures that require imputation prior to learning in order to handle missing values (as theoretically discussed in the previous section), one can use procedures intrinsically able to handle missing values. Among them, decision trees offer a natural way for empirical risk minimization with missing values [Saar-Tsechansky and Provost, 2007, Twala et al., 2008]. This is due to their ability to handle the half-discrete nature of X̃, as they rely on greedy decisions rather than smooth optimization.

We first present the different approaches available to handle missing values in tree-based meth-
ods in Sections B.2 and B.3. We then compare them theoretically in Section B.4, highlighting
the interest of using the “missing incorporated in attribute” approach in particular when the
missing values are informative. In Section 5, we will compare numerically the two different
strategies: imputation prior to learning and learning directly with missing data.
hence

$$
\begin{align*}
\lim_{\rho \to 0} \mathbb{E}[Y | \mathbf{X}' \in B(\mathbf{x}, \rho)] &= 0 \\
&= \mathbb{E}[f(\mathbf{X}) | X_2 = x_2, \dots, X_d = x_d] \\
&= \mathbb{E}[f(\mathbf{X}) | X_2 = x_2, \dots, X_d = x_d, M_1 = 1],
\end{align*}
$$

since $M_1 \perp X_1|(X_2, \dots, X_d)$. Consequently, for all $\mathbf{x} \in [0, 1]^d$ such that $x_1 = \alpha$,

$$
\lim_{\rho \to 0} \mathbb{E}[Y | \mathbf{X}' \in B(\mathbf{x}, \rho)] = \mathbb{E}[f(\mathbf{X}) | X_2 = x_2, \dots, X_d = x_d, M_1 = 1]. \quad (20)
$$

Combining equations (13), (14) and (20), the prediction given by the mean imputation followed
by learning is, for all $\mathbf{x}' \in \mathbb{R}^d$,

$$
\begin{align*}
f_{\mathbf{X}'}^{*}(\mathbf{x}') ={}& \mathbb{E}[Y | X_2 = x_2, \dots, X_d = x_d, M_1 = 1] \mathbf{1}_{x_1' = \alpha} \mathbf{1}_{\mathbb{P}[M_1=1|X_2=x_2, \dots, X_d=x_d]>0} \\
& + \mathbb{E}[Y | \mathbf{X} = \mathbf{x}'] \mathbf{1}_{x_1' = \alpha} \mathbf{1}_{\mathbb{P}[M_1=1|X_2=x_2, \dots, X_d=x_d]=0} \\
& + \mathbb{E}[Y | X_2 = x_2, \dots, X_d = x_d, M_1 = 0] \mathbf{1}_{x_1' \neq \alpha},
\end{align*}
$$

which concludes the proof.

B Decision trees: an example of empirical risk minimization with missing
data

Aside from using procedures that require imputation prior to learning in order to handle missing values (as theoretically discussed in the previous section), one can use procedures intrinsically able to handle missing values. Among them, decision trees offer a natural way for empirical risk minimization with missing values [Saar-Tsechansky and Provost, 2007, Twala et al., 2008]. This is due to their ability to handle the half-discrete nature of X̃, as they rely on greedy decisions rather than smooth optimization.

We first present the different approaches available to handle missing values in tree-based meth-
ods in Sections B.2 and B.3. We then compare them theoretically in Section B.4, highlighting
the interest of using the “missing incorporated in attribute” approach in particular when the
missing values are informative. In Section 5, we will compare numerically the two different
strategies: imputation prior to learning and learning directly with missing data.
## B.1 Tree construction with CART

CART (Classification And Regression Trees, Breiman et al. 1984) is one of the most popular tree algorithm, originally designed for complete data sets. It recursively builds a partition of the input space $\mathcal{X} = \mathbb{R}^d$, and predicts by aggregating the observation labels (average or majority vote) inside each terminal node, also called leaf. For each node $A = \prod_{j=1}^d [a_{j,L}, a_{j,R}] \subset \mathbb{R}^d$, CART algorithm finds the best split $(j^*, z^*)$ among the set of all eligible splits $S = \{(j, z), j \in [1, d]\}$, where $z \in \mathbb{R}$, $z_j \in [a_{j,L}, a_{j,R}]$, where a split is defined by the variable $j$ along which the split is performed and the position $z$ of the split. More precisely, the best split $(j^*, z^*)$ in any node $A$ is the solution of the following optimization problem

$$ (j^*, z^*) \in \underset{(j,z) \in S}{\operatorname{argmin}} \mathbb{E} \left[ \left(Y - \mathbb{E}[Y|X_j \le z, \mathbf{X} \in A]\right)^2 \cdot \mathbb{1}_{X_j \le z, \mathbf{X} \in A} + \left(Y - \mathbb{E}[Y|X_j > z, \mathbf{X} \in A]\right)^2 \cdot \mathbb{1}_{X_j > z, \mathbf{X} \in A} \right]. \quad (21) $$

For each cell $A$, problem (21) can be rewritten as

$$ f^* \in \underset{f \in \mathcal{P}_c}{\operatorname{argmin}} \mathbb{E} \left[ (Y - f(\mathbf{X}))^2 \mathbb{1}_{\mathbf{X} \in A} \right], \quad (22) $$

where $\mathcal{P}_c$ is the set of piecewise-constant functions on $A \cap \{x_j \le s\}$ and $A \cap \{x_j > s\}$ for $(j, s) \in S$. Therefore the optimization problem (22) amounts to solving a least square problem on the subclass of functions $\mathcal{P}_c$. Thus, by minimizing the mean squared error, the CART procedure targets the quantity $\mathbb{E}[Y|\mathbf{X}]$. In the presence of missing values, this criterion must be adapted and several ways to do so have been proposed. Section B.2 and Section B.3 detail the existing criteria that can be used when dealing with missing values.

## B.2 Splitting criterion discarding missing values

A simple option to extend CART methodology in presence of missing values is to select the best split only on the available cases for each variable. More precisely, for any node $A$, the best split in presence of missing values is a solution of the new optimization problem

$$ (j^*, z^*) \in \underset{(j,z) \in S}{\operatorname{argmin}} \mathbb{E} \left[ \left(Y - \mathbb{E}[Y|X_j \le z, M_j = 0, \mathbf{X} \in A]\right)^2 \cdot \mathbb{1}_{X_j \le z, M_j = 0, \mathbf{X} \in A} + \left(Y - \mathbb{E}[Y|X_j > z, M_j = 0, \mathbf{X} \in A]\right)^2 \cdot \mathbb{1}_{X_j > z, M_j = 0, \mathbf{X} \in A} \right], \quad (23) $$

which is nothing but problem (21) computed, for each $j \in \{1, \dots, d\}$, on observed values "$M_j = 0$" only. Note that, by convention, the problem (23) is optimised empirically along
## B.1 Tree construction with CART

CART (Classification And Regression Trees, Breiman et al. 1984) is one of the most popular tree algorithm, originally designed for complete data sets. It recursively builds a partition of the input space $\mathcal{X} = \mathbb{R}^d$, and predicts by aggregating the observation labels (average or majority vote) inside each terminal node, also called leaf. For each node $A = \prod_{j=1}^d [a_{j,L}, a_{j,R}] \subset \mathbb{R}^d$, CART algorithm finds the best split $(j^*, z^*)$ among the set of all eligible splits $S = \{(j, z), j \in [1, d]\}$, where $z \in \mathbb{R}$, $z_j \in [a_{j,L}, a_{j,R}]$, where a split is defined by the variable $j$ along which the split is performed and the position $z$ of the split. More precisely, the best split $(j^*, z^*)$ in any node $A$ is the solution of the following optimization problem

$$ (j^*, z^*) \in \underset{(j,z) \in S}{\operatorname{argmin}} \mathbb{E} \left[ \left(Y - \mathbb{E}[Y|X_j \le z, \mathbf{X} \in A]\right)^2 \cdot \mathbb{1}_{X_j \le z, \mathbf{X} \in A} + \left(Y - \mathbb{E}[Y|X_j > z, \mathbf{X} \in A]\right)^2 \cdot \mathbb{1}_{X_j > z, \mathbf{X} \in A} \right]. \quad (21) $$

For each cell $A$, problem (21) can be rewritten as

$$ f^* \in \underset{f \in \mathcal{P}_c}{\operatorname{argmin}} \mathbb{E} \left[ (Y - f(\mathbf{X}))^2 \mathbb{1}_{\mathbf{X} \in A} \right], \quad (22) $$

where $\mathcal{P}_c$ is the set of piecewise-constant functions on $A \cap \{x_j \le s\}$ and $A \cap \{x_j > s\}$ for $(j, s) \in S$. Therefore the optimization problem (22) amounts to solving a least square problem on the subclass of functions $\mathcal{P}_c$. Thus, by minimizing the mean squared error, the CART procedure targets the quantity $\mathbb{E}[Y|\mathbf{X}]$. In the presence of missing values, this criterion must be adapted and several ways to do so have been proposed. Section B.2 and Section B.3 detail the existing criteria that can be used when dealing with missing values.

## B.2 Splitting criterion discarding missing values

A simple option to extend CART methodology in presence of missing values is to select the best split only on the available cases for each variable. More precisely, for any node $A$, the best split in presence of missing values is a solution of the new optimization problem

$$ (j^*, z^*) \in \underset{(j,z) \in S}{\operatorname{argmin}} \mathbb{E} \left[ \left(Y - \mathbb{E}[Y|X_j \le z, M_j = 0, \mathbf{X} \in A]\right)^2 \cdot \mathbb{1}_{X_j \le z, M_j = 0, \mathbf{X} \in A} + \left(Y - \mathbb{E}[Y|X_j > z, M_j = 0, \mathbf{X} \in A]\right)^2 \cdot \mathbb{1}_{X_j > z, M_j = 0, \mathbf{X} \in A} \right], \quad (23) $$

which is nothing but problem (21) computed, for each $j \in \{1, \dots, d\}$, on observed values "$M_j = 0$" only. Note that, by convention, the problem (23) is optimised empirically along
the variables *j* that have at least two observed entries (we do not allow splits to be performed along variables that have less than two observed entries). This splitting strategy is described in Algorithm 1. As the missing values were not used to calculate the criterion, it is still necessary to specify to which cell they are sent. The solution consisting in discarding missing data at each step would lead to a drastic reduction of the data set and is therefore not viable. The different methods to propagate missing data down the tree are explained below.

<caption>**Algorithm 1 Splitting strategy**</caption>

1: **Input:** a node A, a sample $D_n$ of observations falling into A.
2: For each variable $j \in \{1, \dots, d\}$ and position $z$, compute the CART splitting criterion on observed values "$M_j = 0$" only.
3: Choose a split $(j^*, z^*)$ that minimizes the previous criterion (see optimization problem 23).
4: Split the cell A accordingly. Two new cells $A_L$ and $A_R$ are created.
5: **Output:** Split $(j^*, z^*)$, $A_L$, $A_R$.

**Surrogate splits** Once the best split is chosen via Algorithm 1, surrogate splits search for a split on another variable that induces a data partition close to the original one. More precisely, for a selected split $(j_0^*, z_0^*)$, to observations send down the tree with no $j_0^*$th variable, a new stump, i.e., a tree with one cut, is fitted to the response $1_{X_{j_0^*} \le z_0^*}$, using variables $(X_j)_{j \ne j_0^*}$. The split $(j_1^*, z_1^*)$ which minimizes the misclassification error is selected, and observations are split accordingly. Those that lack both variables $j_0^*$ and $j_1^*$ are split with the second best, $j_2^*$, and so on until the proposed split has a worse misclassification error than the blind rule of sending all remaining missing values to the same daughter, the most populated one. To predict, the training surrogates are kept. This construction is described in Algorithm 2 and is the default method in rpart [Therneau et al., 1997]. Surrogate method is expected to be appropriate when there are relationships between covariates.

**Probabilistic splits** Another option is to propagate missing observations according to a Bernoulli distribution $\mathcal{B}(\frac{\#L}{\#L+\#R})$, where $\#L$ (resp. $\#R$) is the number of points already on the left (resp. right) (see Algorithm 3 in Appendix B.5). This is the default method in C4.5 [Quinlan, 2014].

**Block propagation** A third option is to choose the split on the observed values, and then send all incomplete observations as a block, to a side chosen by minimizing the error (see Algorithm 4 in Appendix B.5). This is the method used in LightGBM [Ke et al., 2017].

Note that Hothorn et al. [2006] proposed conditional trees, a variant of CART which also uses surrogate splits, but adapts the criterion (23) to missing values. Indeed, this criterion
the variables *j* that have at least two observed entries (we do not allow splits to be performed along variables that have less than two observed entries). This splitting strategy is described in Algorithm 1. As the missing values were not used to calculate the criterion, it is still necessary to specify to which cell they are sent. The solution consisting in discarding missing data at each step would lead to a drastic reduction of the data set and is therefore not viable. The different methods to propagate missing data down the tree are explained below.

<caption>**Algorithm 1 Splitting strategy**</caption>

1: **Input:** a node A, a sample $\mathcal{D}_n$ of observations falling into A.
2: For each variable $j \in \{1, \dots, d\}$ and position $z$, compute the CART splitting criterion on observed values "$M_j = 0$" only.
3: Choose a split $(j^*, z^*)$ that minimizes the previous criterion (see optimization problem 23).
4: Split the cell A accordingly. Two new cells $A_L$ and $A_R$ are created.
5: **Output:** Split $(j^*, z^*)$, $A_L$, $A_R$.

*Surrogate splits* Once the best split is chosen via Algorithm 1, surrogate splits search for a split on another variable that induces a data partition close to the original one. More precisely, for a selected split $(j_0^*, z_0^*)$, to observations send down the tree with no $j_0^*$th variable, a new stump, i.e., a tree with one cut, is fitted to the response $1_{X_{j_0^*} \le z_0^*}$, using variables $(X_j)_{j \ne j_0^*}$. The split $(j_1^*, z_1^*)$ which minimizes the misclassification error is selected, and observations are split accordingly. Those that lack both variables $j_0^*$ and $j_1^*$ are split with the second best, $j_2^*$, and so on until the proposed split has a worse misclassification error than the blind rule of sending all remaining missing values to the same daughter, the most populated one. To predict, the training surrogates are kept. This construction is described in Algorithm 2 and is the default method in *rpart* [Therneau et al., 1997]. Surrogate method is expected to be appropriate when there are relationships between covariates.

*Probabilistic splits* Another option is to propagate missing observations according to a Bernoulli distribution $\mathcal{B}(\frac{\#L}{\#L+\#R})$, where $\#L$ (resp. $\#R$) is the number of points already on the left (resp. right) (see Algorithm 3 in Appendix B.5). This is the default method in C4.5 [Quinlan, 2014].

*Block propagation* A third option is to choose the split on the observed values, and then send all incomplete observations as a block, to a side chosen by minimizing the error (see Algorithm 4 in Appendix B.5). This is the method used in LightGBM [Ke et al., 2017].

Note that Hothorn et al. [2006] proposed conditional trees, a variant of CART which also uses surrogate splits, but adapts the criterion (23) to missing values. Indeed, this criterion
implies a selection bias: it leads to underselecting the variables with many missing values due to multiple comparison effects [Strobl et al., 2007]. As a result, it favors variables where many splits are available, and therefore those with fewer missing values. Conditional trees are based on the calculation of a linear statistic of association between Y and each feature $X_j$, $T = (X_j, Y)$. Then, its distribution under the null hypothesis of independence between Y and $X_j$ is estimated by permutation, and the variable with the smallest p-value is selected. As illustrated in Appendix C.1, conditional trees are meant to improve the selection of the splitting variables but do not ensure an improvement in prediction performance.

## B.3 Splitting criterion with missing values: MIA

A second class of methods uses missing values to compute the splitting criterion itself. Consequently, the splitting location depends on the missing values, contrary to all methods presented in Section B.2. Its most common instance is “missing incorporated in attribute” (MIA, Twala et al. 2008, ). More specifically, MIA considers the following splits, for all splits $(j, z)$:

- $\{\tilde{X}_j \le z \text{ or } \tilde{X}_j = \text{NA}\} \text{ vs } \{\tilde{X}_j > z\}$,

- $\{\tilde{X}_j \le z\} \text{ vs } \{\tilde{X}_j > z \text{ or } \tilde{X}_j = \text{NA}\}$,

- $\{\tilde{X}_j \ne \text{NA}\} \text{ vs } \{\tilde{X}_j = \text{NA}\}$.

In a nutshell, for each possible split, MIA tries to send all missing values to the left or to the right, and compute for each choice the corresponding error (right-hand side in 21, as well as

**Algorithm 2 Surrogate strategy**

1: **Input:** a node A, a sample $\mathcal{D}_n$ of observations falling into A, a split $(j_0^*, z_0^*)$ produced by Algorithm 1

2: Create the variable $1_{X_{j_0^*} \le z_0^*}$

3: Let $\mathcal{J} = \{1, \dots, d\} - \{j_0\}$

4: **while** Missing data have not been sent down the tree **do**

5: Compute the misclassification error $\varepsilon^*$ corresponding to sending all remaining missing observation on the most populated side.

6: **for** all $j \in \mathcal{J}$ **do**

7: Fit a tree with one split along j, on data in $\mathcal{D}_n$ with observed values $X_j$ and $X_{j_0^*}$, in order to predict $1_{X_{j_0^*} \le z_0^*}$.

8: For $j \in \mathcal{J}$, let $\varepsilon_j$ be the misclassification error of the previous trees

9: Let $j_{min} \in \text{argmin}_{j \in \mathcal{J}} \varepsilon_j$

10: **end for**

11: **if** $\varepsilon_{j_{min}} < \varepsilon^*$ **then**

12: Use the tree built on $j_{min}$ to send data with missing values on $j_0$ into $A_L$ or $A_R$, depending on the tree prediction

13: $\mathcal{J} \leftarrow \mathcal{J} \cup \{j_{min}\}$

14: **else**

15: Send all remaining missing values to the most populated cell ($A_L$ or $A_R$)

16: **end if**

17: **end while**
implies a selection bias: it leads to underselecting the variables with many missing values due to multiple comparison effects [Strobl et al., 2007]. As a result, it favors variables where many splits are available, and therefore those with fewer missing values. Conditional trees are based on the calculation of a linear statistic of association between Y and each feature $X_j$, $T = (X_j, Y)$. Then, its distribution under the null hypothesis of independence between Y and $X_j$ is estimated by permutation, and the variable with the smallest p-value is selected. As illustrated in Appendix C.1, conditional trees are meant to improve the selection of the splitting variables but do not ensure an improvement in prediction performance.

## B.3 Splitting criterion with missing values: MIA

A second class of methods uses missing values to compute the splitting criterion itself. Consequently, the splitting location depends on the missing values, contrary to all methods presented in Section B.2. Its most common instance is “missing incorporated in attribute” (MIA, Twala et al. 2008, ). More specifically, MIA considers the following splits, for all splits $(j, z)$:

- $\{\tilde{X}_j \le z \text{ or } \tilde{X}_j = \text{NA}\} \text{ vs } \{\tilde{X}_j > z\}$,

- $\{\tilde{X}_j \le z\} \text{ vs } \{\tilde{X}_j > z \text{ or } \tilde{X}_j = \text{NA}\}$,

- $\{\tilde{X}_j \ne \text{NA}\} \text{ vs } \{\tilde{X}_j = \text{NA}\}$.

In a nutshell, for each possible split, MIA tries to send all missing values to the left or to the right, and compute for each choice the corresponding error (right-hand side in 21, as well as

<caption>**Algorithm 2** Surrogate strategy</caption>

1: **Input:** a node A, a sample $\mathcal{D}_n$ of observations falling into A, a split $(j_0^*, z_0^*)$ produced by Algorithm 1
2: Create the variable $\mathbb{1}_{X_{j_0}^* \le z_0^*}$
3: Let $\mathcal{J} = \{1, \dots, d\} - \{j_0\}$
4: **while** Missing data have not been sent down the tree **do**
5: Compute the misclassification error $\varepsilon^*$ corresponding to sending all remaining missing observation on the most populated side.
6: **for** all $j \in \mathcal{J}$ **do**
7: Fit a tree with one split along j, on data in $\mathcal{D}_n$ with observed values $X_j$ and $X_{j_0}^*$, in order to predict $\mathbb{1}_{X_{j_0}^* \le z_0^*}$.
8: For $j \in \mathcal{J}$, let $\varepsilon_j$ be the misclassification error of the previous trees
9: Let $j_{min} \in \text{argmin}_{j \in \mathcal{J}} \varepsilon_j$
10: **end for**
11: **if** $\varepsilon_{j_{min}} < \varepsilon^*$ **then**
12: Use the tree built on $j_{min}$ to send data with missing values on $j_0$ into $A_L$ or $A_R$, depending on the tree prediction
13: $\mathcal{J} \leftarrow \mathcal{J} \cup \{j_{min}\}$
14: **else**
15: Send all remaining missing values to the most populated cell ($A_L$ or $A_R$)
16: **end if**
17: **end while**
the error associated to separating the observed values from the missing ones. Finally, it chooses
the split among the previous ones with the lowest error (see Algorithm 5 in Appendix B.5).
Note that block propagation can be seen as a greedy way of successively optimizing the choices
in two first options. However, as we show in Prop. 1, these successive choices are sub-optimal.
Missing values are treated as a category by MIA, which is thus nothing but a greedy algorithm
minimizing the square loss between Y and a function of $\tilde{X}$ and consequently targets the
quantity (4) which separate $E[Y|\tilde{X}]$ into $2^d$ terms. However, it is not exhaustive: at each step,
the tree can cut for each variable according to missing or non missing and selects this cut when
it is relevant, i.e. when it minimizes the prediction error. The final leaves can correspond to a
cluster of missing values patterns (observations with missing values on the two first variables
for instance and any missing patterns for the other variables).

MIA is thought to be a good method to apply when missing pattern is informative, as this
procedure allows to cut with respect to missing/ non missing and uses missing data to compute
the best splits. Note this latter property implies that the MIA approach does not require a
different method to propagate missing data down the tree. Notably, MIA is implemented in the
R packages partykit [Hothorn and Zeileis, 2015] and grf [Tibshirani et al., 2020], as well as in
XGBoost [Chen and Guestrin, 2016] and for the HistGradientBoosting models in scikit-learn
[Pedregosa et al., 2011].

**Remark 2 Implementation:** A simple way to implement MIA consists in duplicating the incomplete columns, and replacing the missing entries once by $+\infty$ and once by $-\infty$ (or an extreme out-of-range value). This creates two dummy variables for each original one containing missing values. Splitting along a variable and sending all missing data to the left (for example) is the same as splitting along the corresponding dummy variable where missing entries have been completed by $-\infty$. Alternatively, MIA can be with two scans on a feature's values in ascending and descending orders [Chen and Guestrin, 2016, Alg 3].

**Remark 3 Implicit imputation:** Whether it is in the case where the missing values are propagated in the available case method (Section B.2), or incorporated in the split choice in MIA, missing values are assigned either to the left or the right interval. Consequently, handling missing values in a tree can be seen as implicit imputation by an interval value.

## B.4 Theoretical comparison of CART versus MIA

We now compare theoretically the positions of the splitting point at the root and the prediction
errors on simple examples with MCAR values. Proposition 1 computes the splitting position
of MIA and CART, and highlights that the splitting position of MIA varies even for MCAR
the error associated to separating the observed values from the missing ones. Finally, it chooses
the split among the previous ones with the lowest error (see Algorithm 5 in Appendix B.5).
Note that block propagation can be seen as a greedy way of successively optimizing the choices
in two first options. However, as we show in Prop. 1, these successive choices are sub-optimal.
Missing values are treated as a category by MIA, which is thus nothing but a greedy algorithm
minimizing the square loss between Y and a function of $\tilde{X}$ and consequently targets the
quantity (4) which separate $E[Y|\tilde{X}]$ into $2^d$ terms. However, it is not exhaustive: at each step,
the tree can cut for each variable according to missing or non missing and selects this cut when
it is relevant, i.e. when it minimizes the prediction error. The final leaves can correspond to a
cluster of missing values patterns (observations with missing values on the two first variables
for instance and any missing patterns for the other variables).

MIA is thought to be a good method to apply when missing pattern is informative, as this
procedure allows to cut with respect to missing/ non missing and uses missing data to compute
the best splits. Note this latter property implies that the MIA approach does not require a
different method to propagate missing data down the tree. Notably, MIA is implemented in the
R packages partykit [Hothorn and Zeileis, 2015] and grf [Tibshirani et al., 2020], as well as in
XGBoost [Chen and Guestrin, 2016] and for the HistGradientBoosting models in scikit-learn
[Pedregosa et al., 2011].

**Remark 2 Implementation:** A simple way to implement MIA consists in duplicating the incomplete columns, and replacing the missing entries once by $+\infty$ and once by $-\infty$ (or an extreme out-of-range value). This creates two dummy variables for each original one containing missing values. Splitting along a variable and sending all missing data to the left (for example) is the same as splitting along the corresponding dummy variable where missing entries have been completed by $-\infty$. Alternatively, MIA can be with two scans on a feature's values in ascending and descending orders [Chen and Guestrin, 2016, Alg 3].

**Remark 3 Implicit imputation:** Whether it is in the case where the missing values are propagated in the available case method (Section B.2), or incorporated in the split choice in MIA, missing values are assigned either to the left or the right interval. Consequently, handling missing values in a tree can be seen as implicit imputation by an interval value.

## B.4 Theoretical comparison of CART versus MIA

We now compare theoretically the positions of the splitting point at the root and the prediction
errors on simple examples with MCAR values. Proposition 1 computes the splitting position
of MIA and CART, and highlights that the splitting position of MIA varies even for MCAR
Fig. 5: Split position chosen by MIA and CART criterion, depending on the fraction $p$ of missing values on $X_1$.

<figure><img src="image_6.png" /></figure>

missing data. Proposition 2 then compares the risk of the different splitting strategies: probabilistic split, block propagation, surrogate split, and MIA. We prove that MIA and surrogate splits are the two best strategies, one of which may be better than the other depending on the dependence structure of covariates.

**Proposition 1** Let $p \in [0, 1]$. Consider the regression model

$$ \begin{cases} Y = X_1 \\ X_1 \sim U([0, 1]) \end{cases}, \quad \begin{cases} \mathbb{P}[M_1 = 0] = 1 - p \\ \mathbb{P}[M_1 = 1] = p \end{cases}, $$

where $M_1 \perp (X_1, Y)$ is the missingness pattern on $X_1$. Let $C_{MIA}(1, s, q, p)$ be the value of the splitting MIA criterion computed on $X_1$ at threshold $s$ such that $(1, s) \in S$, and $q \in \{L, R\}$, where $q$ stands for the side where missing values are sent. Therefore,

1. The best split $s^*$ given by the CART criterion (23) is $s^* = 1/2$.

2. The best splits $s^*_{MIA,L}(p)$ and $s^*_{MIA,R}(p)$ given by the MIA procedure (described in Section B.3), assuming that all missing values are sent to the left node (resp. to the right node), satisfy

$$ s^*_{MIA,L}(p) = \underset{s \in [0,1]}{\operatorname{argmin}} C_{MIA}(1, s, L, p), \quad (24) $$

where

$$ C_{MIA}(1, s, L, p) = \frac{1}{3} - \frac{1}{p + (1-p)s} \left( \frac{p}{2} + \frac{(1-p)s^2}{2} \right)^2 - (1-p)(1-s) \left( \frac{1+s}{2} \right)^2, $$

and $s^*_{MIA,R}(p) = 1 - s^*_{MIA,L}(p)$.

The proof is given in Appendix A. Proposition 1 shows that the split given by optimizing the CART criterion does not depend on the percentage $p$ of missing values since the pattern is independent of $(X, Y)$. A numerical solution to equation (24) is displayed in Figure 5. When there are no missing values ($p = 0$), the split occur at $s = 1/2$ as expected. When $p$ increases, the threshold does not correspond anymore to the one calculated using observed values only
Fig. 5: Split position chosen by MIA and CART criterion, depending on the fraction $p$ of missing values on $X_1$.

<figure><img src="image_7.png" /></figure>

missing data. Proposition 2 then compares the risk of the different splitting strategies: probabilistic split, block propagation, surrogate split, and MIA. We prove that MIA and surrogate splits are the two best strategies, one of which may be better than the other depending on the dependence structure of covariates.

**Proposition 1** Let $p \in [0, 1]$. Consider the regression model

$$ \begin{cases} Y = X_1 \\ X_1 \sim U([0, 1]) \end{cases}, \quad \begin{cases} \mathbb{P}[M_1 = 0] = 1-p, \\ \mathbb{P}[M_1 = 1] = p \end{cases}, $$

where $M_1 \perp (X_1, Y)$ is the missingness pattern on $X_1$. Let $C_{MIA}(1, s, q, p)$ be the value of the splitting MIA criterion computed on $X_1$ at threshold $s$ such that $(1, s) \in S$, and $q \in \{L, R\}$, where $q$ stands for the side where missing values are sent. Therefore,

1. The best split $s^*$ given by the CART criterion (23) is $s^* = 1/2$.

2. The best splits $s^*_{MIA,L}(p)$ and $s^*_{MIA,R}(p)$ given by the MIA procedure (described in Section B.3), assuming that all missing values are sent to the left node (resp. to the right node), satisfy

$$ s^*_{MIA,L}(p) = \underset{s \in [0,1]}{\operatorname{argmin}} C_{MIA}(1, s, L, p), \quad (24) $$

where

$$ C_{MIA}(1, s, L, p) = \frac{1}{3} - \frac{1}{p + (1-p)s} \left( \frac{p}{2} + \frac{(1-p)s^2}{2} \right)^2 - (1-p)(1-s) \left( \frac{1+s}{2} \right)^2, $$

and $s^*_{MIA,R}(p) = 1 - s^*_{MIA,L}(p)$.

The proof is given in Appendix A. Proposition 1 shows that the split given by optimizing the CART criterion does not depend on the percentage $p$ of missing values since the pattern is independent of $(X, Y)$. A numerical solution to equation (24) is displayed in Figure 5. When there are no missing values ($p = 0$), the split occur at $s = 1/2$ as expected. When $p$ increases, the threshold does not correspond anymore to the one calculated using observed values only
as it is influenced by the missing entries even in the MCAR setting. It may be surprising that the splitting position change in a MCAR setting but the missing values correspond to data from the whole interval [0, 1] and thus introduce noise in the side they are sent to. This implies that the cell that receives missing values must be bigger than usual so that the noise of the missing data is of the same order than the noise of original observations belonging to the cell. Recall that, since the threshold in MIA is chosen by taking into account missing values, it is straightforward to propagate a new element with missing values down the tree.

Recall that the quadratic risk $R$ of a function $f^*$ is defined as $R(f^*) = E[(Y - f^*(X))^2]$. Proposition 2 enables us to compare the risk of a tree with a single split computed with the different strategies. It highlights that even in the simple case of MCAR, MIA gives more accurate predictions than block propagation or probabilistic split.

**Proposition 2** Consider the regression model

$$ \left\{ \begin{array}{l} Y = X_1 \\ X_1 \sim U([0, 1]), \\ X_2 = X_1 \mathbb{1}_{W=1} \end{array} \right. , \quad \left\{ \begin{array}{l} \mathbb{P}[W=0] = \eta \\ \mathbb{P}[W=1] = 1-\eta \end{array} \right. , \quad \left\{ \begin{array}{l} \mathbb{P}[M_1=0] = 1-p, \\ \mathbb{P}[M_1=1] = p, \end{array} \right. $$

where $(M_1, W) \perp (X_1, Y)$. The random variable $M_1$ is the pattern of missingness for $X_1$ and $W$ stands for the link between $X_1$ and $X_2$. Let $f_{MIA}^*$, $f_{block}^*$, $f_{prob}^*$, $f_{surr}^*$ be respectively, the theoretical prediction resulting from one split according to MIA, CART with block propagation and CART with probabilistic splitting strategy, and a single split, where missing data are handled via surrogate split (in the infinite sample setting). We have

$$ R(f_{MIA}^*) = \min_{s \in [0,1]} C_{MIA}(1, s, L, p) \mathbb{1}_{p \le \eta} + \min_{s \in [0,1]} C_{MIA}(1, s, L, \eta) \mathbb{1}_{p > \eta}, $$

$$ R(f_{block}^*) = C_{MIA}(1, 1/2, L, p) = C_{MIA}(1, 1/2, R, p) $$

$$ R(f_{prob}^*) = -\frac{p^2}{16} + \frac{p}{8} + \frac{1}{48}, $$

$$ R(f_{surr}^*) = \frac{1}{48} + \frac{6}{48} \eta p. $$

where $C_{MIA}(1, s, L, p)$ is defined in Proposition 1. In particular,

$$ R(f_{MIA}^*) \le R(f_{block}^*) \quad \text{and} \quad R(f_{MIA}^*) \le R(f_{prob}^*). $$

Proof is given in Appendix A. Figure 6 depicts the risk of each estimate, in the context of proposition 2, resulting from a split computed via one of the four methods described above. Only surrogate and MIA risks depend on the value $\eta$ which measures the independence between $X_1$ and $X_2$. As proved, the risk of probabilistic split and block propagation is larger than that
as it is influenced by the missing entries even in the MCAR setting. It may be surprising that the splitting position change in a MCAR setting but the missing values correspond to data from the whole interval [0, 1] and thus introduce noise in the side they are sent to. This implies that the cell that receives missing values must be bigger than usual so that the noise of the missing data is of the same order than the noise of original observations belonging to the cell. Recall that, since the threshold in MIA is chosen by taking into account missing values, it is straightforward to propagate a new element with missing values down the tree.

Recall that the quadratic risk R of a function $f^*$ is defined as $R(f^*) = \mathbb{E}[(Y - f^*(X))^2]$. Proposition 2 enables us to compare the risk of a tree with a single split computed with the different strategies. It highlights that even in the simple case of MCAR, MIA gives more accurate predictions than block propagation or probabilistic split.

**Proposition 2** Consider the regression model

$$ \left\{ \begin{array}{l} Y = X_1 \\ X_1 \sim U([0, 1]), \\ X_2 = X_1 \mathbb{1}_{W=1} \end{array} \right. , \left\{ \begin{array}{l} \mathbb{P}[W=0] = \eta \\ \mathbb{P}[W=1] = 1-\eta \end{array} \right. , \left\{ \begin{array}{l} \mathbb{P}[M_1=0] = 1-p, \\ \mathbb{P}[M_1=1] = p, \end{array} \right. $$

where $(M_1, W) \perp (X_1, Y)$. The random variable $M_1$ is the pattern of missingness for $X_1$ and $W$ stands for the link between $X_1$ and $X_2$. Let $f_{MIA}^*$, $f_{block}^*$, $f_{prob}^*$, $f_{surr}^*$ be respectively, the theoretical prediction resulting from one split according to MIA, CART with block propagation and CART with probabilistic splitting strategy, and a single split, where missing data are handled via surrogate split (in the infinite sample setting). We have

$$ R(f_{MIA}^*) = \min_{s \in [0,1]} C_{MIA}(1, s, L, p) \mathbb{1}_{p \le \eta} + \min_{s \in [0,1]} C_{MIA}(1, s, L, \eta) \mathbb{1}_{p > \eta}, $$

$$ R(f_{block}^*) = C_{MIA}(1, 1/2, L, p) = C_{MIA}(1, 1/2, R, p) $$

$$ R(f_{prob}^*) = -\frac{p^2}{16} + \frac{p}{8} + \frac{1}{48}, $$

$$ R(f_{surr}^*) = \frac{1}{48} + \frac{6}{48} \eta p. $$

where $C_{MIA}(1, s, L, p)$ is defined in Proposition 1. In particular,

$$ R(f_{MIA}^*) \le R(f_{block}^*) \quad \text{and} \quad R(f_{MIA}^*) \le R(f_{prob}^*). $$

Proof is given in Appendix A. Figure 6 depicts the risk of each estimate, in the context of proposition 2, resulting from a split computed via one of the four methods described above. Only surrogate and MIA risks depend on the value $\eta$ which measures the independence between $X_1$ and $X_2$. As proved, the risk of probabilistic split and block propagation is larger than that
<figure><img src="image_7.png" /><figcaption>Fig. 6: Theoretical risk of the splitting methods as a function of $p$, for three values of $\eta$ parameter that controls the amount of coupling between $X_1$ and $X_2$ in the model of Proposition 2.</figcaption></figure>

of MIA. Besides, surrogate split is better than MIA if the link between $X_1$ and $X_2$ is strong (small values of $\eta$) and worse if this link is weak (high values of $\eta$).

## B.5 Splitting algorithms

**Algorithm 3 Probabilistic split**

1: **Input:** a node $A$, a sample $D_n$ of observations falling into $A$, a split $(j_0^*, z_0^*)$ produced by Algorithm 1

2: Compute the number $n_L$ of points with observed $X_{j_0^*}$ falling into $A_L$.

3: Compute the number $n_R$ of points with observed $X_{j_0^*}$ falling into $A_R$.

4: **for** all data with missing value along $j_0^*$ **do**

5: Send the data randomly to $A_L$ (resp. $A_R$) with probability $n_L/(n_L + n_R)$ (resp. $n_R/(n_L + n_R)$)

6: **end for**

**Algorithm 4 Block propagation**

1: **Input:** a node $A$, a sample $D_n$ of observations falling into $A$, a split $(j_0^*, z_0^*)$ produced by Algorithm 1

2: Consider sending all observations with missing values on $j_0^*$ into $A_L$. Compute the corresponding error (criterion on the right-hand side in 21).

3: Consider sending all observations with missing values on $j_0^*$ into $A_R$. Compute the corresponding error (criterion on the right-hand side in 21).

4: Choose the alternative with the lowest error and send all missing data on the same side accordingly.
<figure><img src="image_8.png" /><figcaption>Fig. 6: Theoretical risk of the splitting methods as a function of $p$, for three values of $\eta$ parameter that controls the amount of coupling between $X_1$ and $X_2$ in the model of Proposition 2.</figcaption></figure>

of MIA. Besides, surrogate split is better than MIA if the link between $X_1$ and $X_2$ is strong (small values of $\eta$) and worse if this link is weak (high values of $\eta$).

## B.5 Splitting algorithms

**Algorithm 3 Probabilistic split**

1: **Input:** a node $A$, a sample $\mathcal{D}_n$ of observations falling into $A$, a split $(j_0^*, z_0^*)$ produced by Algorithm 1

2: Compute the number $n_L$ of points with observed $X_{j_0^*}$ falling into $A_L$.

3: Compute the number $n_R$ of points with observed $X_{j_0^*}$ falling into $A_R$.

4: **for** all data with missing value along $j_0^*$ **do**

5: Send the data randomly to $A_L$ (resp. $A_R$) with probability $n_L/(n_L + n_R)$ (resp. $n_R/(n_L + n_R)$)

6: **end for**

**Algorithm 4 Block propagation**

1: **Input:** a node $A$, a sample $\mathcal{D}_n$ of observations falling into $A$, a split $(j_0^*, z_0^*)$ produced by Algorithm 1

2: Consider sending all observations with missing values on $j_0^*$ into $A_L$. Compute the corresponding error (criterion on the right-hand side in 21).

3: Consider sending all observations with missing values on $j_0^*$ into $A_R$. Compute the corresponding error (criterion on the right-hand side in 21).

4: Choose the alternative with the lowest error and send all missing data on the same side accordingly.
**Algorithm 5** Missing Incorporated in Attribute (MIA)

1: **Input:** a node A, a sample $\mathcal{D}_n$ of observations falling into A

2: **for** each split $(j, z)$ **do**

3: Send all observations with missing values on j on the left side. Compute the error $\epsilon_{j,z,L}$ (right-hand side in 21)

4: Send all observations with missing values on j on the right side. Compute the error $\epsilon_{j,z,R}$ (right-hand side in 21)

5: **end for**

6: **for** each $j \in \{1, \dots, d\}$ **do**

7: Compute the error $\epsilon_{j,-,\cdot}$ associated to separating observations with missing data on j from the remaining ones.

8: **end for**

9: Choose the split corresponding to the lowest error $\epsilon_{\cdot, \cdot, \cdot}$. Split the data accordingly.

## B.6 Proof of Proposition 1

**Cart splitting criterion.** Under the model given in Proposition 1, simple calculations show that

$$ 
\begin{align*}
\mathbb{E}[Y|X \in [0, s]] &= \frac{s}{2}, \quad \mathbb{E}[Y^2|X \in [0, s]] = \frac{s^2}{3} \\
\mathbb{E}[Y|X \in [s, 1]] &= \frac{1+s}{2}, \quad \mathbb{E}[Y^2|X \in [s, 1]] = \frac{1-s^3}{3(1-s)} \\
\mathbb{P}[X \in [0, s]] &= s, \quad \mathbb{P}[X \in [s, 1]] = 1-s.
\end{align*}
$$

Thus the CART splitting criterion can be written as

$$ 
\begin{align*}
C(1, s) &= \mathbb{E}[Y^2] - (\mathbb{P}[X \in [0, s]](\mathbb{E}[Y|X \in [0, s]])^2 + \mathbb{P}[X \in [s, 1]](\mathbb{E}[Y|X \in [s, 1]])^2) \\
&= \frac{1}{3} - \left(s\left(\frac{s}{2}\right)^2 + (1-s)\left(\frac{1+s}{2}\right)^2\right) \\
&= \frac{s(s-1)}{4} + \frac{1}{12}.
\end{align*}
$$

By definition,

$$ s^* = \underset{s \in [0,1]}{\operatorname{argmin}} \left( \frac{1}{4} s(s-1) + \frac{1}{12} \right) = 1/2, $$

and the criterion evaluated in $s = 1/2$ is equal to $1/48$. The calculations are exactly the same when a percentage of missing value is added if $M_1 \perp X_1$.

**MIA splitting criterion.** By symmetry, we can assume that missing values are sent left. It is equivalent to observing

$$ X' = 01_{M=1} + X1_{M=0}. $$
**Algorithm 5** Missing Incorporated in Attribute (MIA)

1: **Input:** a node A, a sample $\mathcal{D}_n$ of observations falling into A

2: **for** each split $(j, z)$ **do**

3: Send all observations with missing values on j on the left side. Compute the error $\epsilon_{j,z,L}$ (right-hand side in 21)

4: Send all observations with missing values on j on the right side. Compute the error $\epsilon_{j,z,R}$ (right-hand side in 21)

5: **end for**

6: **for** each $j \in \{1, \dots, d\}$ **do**

7: Compute the error $\epsilon_{j,-,\cdot}$ associated to separating observations with missing data on j from the remaining ones.

8: **end for**

9: Choose the split corresponding to the lowest error $\epsilon_{\cdot, \cdot, \cdot}$. Split the data accordingly.

## B.6 Proof of Proposition 1

**Cart splitting criterion.** Under the model given in Proposition 1, simple calculations show that

$$ 
\begin{aligned}
\mathbb{E}[Y|X \in [0, s]] &= \frac{s}{2}, \quad \mathbb{E}[Y^2|X \in [0, s]] = \frac{s^2}{3} \\
\mathbb{E}[Y|X \in [s, 1]] &= \frac{1+s}{2}, \quad \mathbb{E}[Y^2|X \in [s, 1]] = \frac{1-s^3}{3(1-s)} \\
\mathbb{P}[X \in [0, s]] &= s, \quad \mathbb{P}[X \in [s, 1]] = 1-s.
\end{aligned}
$$

Thus the CART splitting criterion can be written as

$$ 
\begin{aligned}
C(1, s) &= \mathbb{E}[Y^2] - (\mathbb{P}[X \in [0, s]](\mathbb{E}[Y|X \in [0, s]])^2 + \mathbb{P}[X \in [s, 1]](\mathbb{E}[Y|X \in [s, 1]])^2) \\
&= \frac{1}{3} - \left(s\left(\frac{s}{2}\right)^2 + (1-s)\left(\frac{1+s}{2}\right)^2\right) \\
&= \frac{s(s-1)}{4} + \frac{1}{12}.
\end{aligned}
$$

By definition,

$$ s^* = \underset{s \in [0,1]}{\operatorname{argmin}} \left( \frac{1}{4}s(s-1) + \frac{1}{12} \right) = 1/2, $$

and the criterion evaluated in $s = 1/2$ is equal to $1/48$. The calculations are exactly the same when a percentage of missing value is added if $M_1 \perp X_1$.

**MIA splitting criterion.** By symmetry, we can assume that missing values are sent left. It is equivalent to observing

$$ X' = 01_{M=1} + X1_{M=0}. $$
The MIA splitting criterion is then defined as

$$
\begin{align*}
s_{\text{MIA,L}}^* &= \underset{s \in [0,1]}{\operatorname{argmin}} \mathbb{E} \left[ (Y - \mathbb{E}[Y|X' \le s] \mathbf{1}_{X' \le s} - \mathbb{E}[Y|X' > s] \mathbf{1}_{X' > s})^2 \right] \\
&= \underset{s \in [0,1]}{\operatorname{argmin}} \mathbb{P}(X' \le s) \mathbb{E} \left[ (Y - \mathbb{E}[Y|X' \le s])^2 \mathbf{1}_{X' \le s} \right] \\
&\quad + \mathbb{P}(X' > s) \mathbb{E} \left[ (Y - \mathbb{E}[Y|X' > s])^2 \mathbf{1}_{X' > s} \right].
\end{align*}
$$

We have

$$
\begin{align*}
\mathbb{E}[Y'|X' \in [0, s]] &= \mathbb{E}[X|X' \in [0, s]] \\
&= \mathbb{E}[X \mathbf{1}_{M=1} + X \mathbf{1}_{M=0} | X' \in [0, s]] \\
&= \frac{1}{\mathbb{P}[X' \in [0, s]]} \mathbb{E}[X \mathbf{1}_{M=1, X' \in [0, s]} + X \mathbf{1}_{M=0, X' \in [0, s]}] \\
&= \frac{1}{p + (1-p)s} \left( \frac{p}{2} + \frac{(1-p)s^2}{2} \right).
\end{align*}
$$

Besides,

$$
\begin{align*}
\mathbb{E}[Y^2|X' \in [0, s]] &= \mathbb{E}[X^2|X' \in [0, s]] \\
&= \mathbb{E}[X^2 \mathbf{1}_{M=1} + X^2 \mathbf{1}_{M=0} | X' \in [0, s]] \\
&= \frac{1}{p + (1-p)s} \mathbb{E}[X^2 \mathbf{1}_{M=1, X' \in [0, s]} + X^2 \mathbf{1}_{M=0, X' \in [0, s]}] \\
&= \frac{1}{p + (1-p)s} \left( \frac{p}{3} + \frac{(1-p)s^3}{3} \right)
\end{align*}
$$

Thus the left-part of the criterion is given by

$$
\begin{align*}
& \mathbb{P}(X' \in [0, s]) \mathbb{E}[(Y - \mathbb{E}[Y|X' \in [0, s]])^2 | X' \in [0, s]] \\
&= \left( p + (1-p)s \right) \left( \mathbb{E}[Y^2|X' \in [0, s]] - \left( \mathbb{E}[Y|X' \in [0, s]] \right)^2 \right) \\
&= \left( p + (1-p)s \right) \left( \frac{1}{p + (1-p)s} \left( \frac{p}{3} + \frac{(1-p)s^3}{3} \right) \right. \\
&\qquad \left. - \left( \frac{1}{p + (1-p)s} \left( \frac{p}{2} + \frac{(1-p)s^2}{2} \right) \right)^2 \right) \\
&= \left( \frac{p}{3} + \frac{(1-p)s^3}{3} \right) - \frac{1}{p + (1-p)s} \left( \frac{p}{2} + \frac{(1-p)s^2}{2} \right)^2
\end{align*}
$$
The MIA splitting criterion is then defined as

$$
\begin{align*}
s_{\text{MIA,L}}^{*} &= \underset{s \in [0,1]}{\operatorname{argmin}} \mathbb{E} \left[ (Y - \mathbb{E}[Y|X' \le s] \mathbf{1}_{X' \le s} - \mathbb{E}[Y|X' > s] \mathbf{1}_{X' > s})^2 \right] \\
&= \underset{s \in [0,1]}{\operatorname{argmin}} \mathbb{P}(X' \le s) \mathbb{E} \left[ (Y - \mathbb{E}[Y|X' \le s])^2 \mathbf{1}_{X' \le s} \right] \\
&\quad + \mathbb{P}(X' > s) \mathbb{E} \left[ (Y - \mathbb{E}[Y|X' > s])^2 \mathbf{1}_{X' > s} \right].
\end{align*}
$$

We have

$$
\begin{align*}
\mathbb{E}[Y'|X' \in [0, s]] &= \mathbb{E}[X|X' \in [0, s]] \\
&= \mathbb{E}[X \mathbf{1}_{M=1} + X \mathbf{1}_{M=0} | X' \in [0, s]] \\
&= \frac{1}{\mathbb{P}[X' \in [0, s]]} \mathbb{E}[X \mathbf{1}_{M=1, X' \in [0, s]} + X \mathbf{1}_{M=0, X' \in [0, s]}] \\
&= \frac{1}{p + (1-p)s} \left( \frac{p}{2} + \frac{(1-p)s^2}{2} \right).
\end{align*}
$$

Besides,

$$
\begin{align*}
\mathbb{E}[Y^2|X' \in [0, s]] &= \mathbb{E}[X^2|X' \in [0, s]] \\
&= \mathbb{E}[X^2 \mathbf{1}_{M=1} + X^2 \mathbf{1}_{M=0} | X' \in [0, s]] \\
&= \frac{1}{p + (1-p)s} \mathbb{E}[X^2 \mathbf{1}_{M=1, X' \in [0, s]} + X^2 \mathbf{1}_{M=0, X' \in [0, s]}] \\
&= \frac{1}{p + (1-p)s} \left( \frac{p}{3} + \frac{(1-p)s^3}{3} \right)
\end{align*}
$$

Thus the left-part of the criterion is given by

$$
\begin{align*}
& \mathbb{P}(X' \in [0, s]) \mathbb{E}[(Y - \mathbb{E}[Y|X' \in [0, s]])^2 | X' \in [0, s]] \\
&= \left( p + (1-p)s \right) \left( \mathbb{E}[Y^2|X' \in [0, s]] - \left( \mathbb{E}[Y|X' \in [0, s]] \right)^2 \right) \\
&= \left( p + (1-p)s \right) \left( \frac{1}{p + (1-p)s} \left( \frac{p}{3} + \frac{(1-p)s^3}{3} \right) \right. \\
&\qquad \left. - \left( \frac{1}{p + (1-p)s} \left( \frac{p}{2} + \frac{(1-p)s^2}{2} \right) \right)^2 \right) \\
&= \left( \frac{p}{3} + \frac{(1-p)s^3}{3} \right) - \frac{1}{p + (1-p)s} \left( \frac{p}{2} + \frac{(1-p)s^2}{2} \right)^2
\end{align*}
$$
On the other hand, we have

$$
\begin{align*}
\mathbb{E}[Y|X' \in [s, 1]] &= \mathbb{E}[X|X' \in [s, 1]] \\
&= \mathbb{E}[X1_{M=1} + X1_{M=0}|X' \in [s, 1]] \\
&= \frac{1}{(1-p)(1-s)} \mathbb{E}[X1_{M=1, X' \in [s, 1]} + X1_{M=0, X' \in [s, 1]}] \\
&= \frac{1}{(1-p)(1-s)} \left( (1-p) \frac{1-s^2}{2} \right) \\
&= \frac{1+s}{2}.
\end{align*}
$$

Besides,

$$
\begin{align*}
\mathbb{E}[Y^2|X' \in [s, 1]] &= \mathbb{E}[X^2|X' \in [s, 1]] \\
&= \mathbb{E}[X^21_{M=1} + X^21_{M=0}|X' \in [s, 1]] \\
&= \frac{1}{(1-p)(1-s)} \mathbb{E}[X^21_{M=1, X' \in [s, 1]} + X^21_{M=0, X' \in [s, 1]}] \\
&= \frac{1}{(1-p)(1-s)} \left( (1-p) \frac{1-s^3}{3} \right) \\
&= \frac{1-s^3}{3(1-s)}.
\end{align*}
$$

Thus the right-part of the criterion is given by

$$
\begin{align*}
& \mathbb{P}(X' \in [s, 1]) \mathbb{E}[(Y - \mathbb{E}[Y|X' \in [s, 1]])^2 | X' \in [s, 1]] \\
&= \left( (1-p)(1-s) \right) \left( \mathbb{E}[Y^2 | X' \in [s, 1]] - \left( \mathbb{E}[Y | X' \in [s, 1]] \right)^2 \right) \\
&= \left( (1-p)(1-s) \right) \left( \frac{1-s^3}{3(1-s)} - \left( \frac{1+s}{2} \right)^2 \right) \\
&= (1-p) \frac{1-s^3}{3} - (1-p)(1-s) \left( \frac{1+s}{2} \right)^2.
\end{align*}
$$

Finally,

$$
s_{\text{MIA,L}}^* = \underset{s \in [0,1]}{\operatorname{argmin}} \left\{ \left( \frac{p}{3} + \frac{(1-p)s^3}{3} \right) - \frac{1}{p + (1-p)s} \left( \frac{p}{2} + \frac{(1-p)s^2}{2} \right)^2 \right. \\
\left. + (1-p) \frac{1-s^3}{3} - (1-p)(1-s) \left( \frac{1+s}{2} \right)^2 \right\},
$$

which concludes the proof.
On the other hand, we have

$$
\begin{align*}
\mathbb{E}[Y|X' \in [s, 1]] &= \mathbb{E}[X|X' \in [s, 1]] \\
&= \mathbb{E}[X1_{M=1} + X1_{M=0}|X' \in [s, 1]] \\
&= \frac{1}{(1-p)(1-s)} \mathbb{E}[X1_{M=1, X' \in [s, 1]} + X1_{M=0, X' \in [s, 1]}] \\
&= \frac{1}{(1-p)(1-s)} \left( (1-p) \frac{1-s^2}{2} \right) \\
&= \frac{1+s}{2}.
\end{align*}
$$

Besides,

$$
\begin{align*}
\mathbb{E}[Y^2|X' \in [s, 1]] &= \mathbb{E}[X^2|X' \in [s, 1]] \\
&= \mathbb{E}[X^21_{M=1} + X^21_{M=0}|X' \in [s, 1]] \\
&= \frac{1}{(1-p)(1-s)} \mathbb{E}[X^21_{M=1, X' \in [s, 1]} + X^21_{M=0, X' \in [s, 1]}] \\
&= \frac{1}{(1-p)(1-s)} \left( (1-p) \frac{1-s^3}{3} \right) \\
&= \frac{1-s^3}{3(1-s)}.
\end{align*}
$$

Thus the right-part of the criterion is given by

$$
\begin{align*}
& \mathbb{P}(X' \in [s, 1]) \mathbb{E}[(Y - \mathbb{E}[Y|X' \in [s, 1]])^2 | X' \in [s, 1]] \\
&= \left( (1-p)(1-s) \right) \left( \mathbb{E}[Y^2 | X' \in [s, 1]] - \left( \mathbb{E}[Y | X' \in [s, 1]] \right)^2 \right) \\
&= \left( (1-p)(1-s) \right) \left( \frac{1-s^3}{3(1-s)} - \left( \frac{1+s}{2} \right)^2 \right) \\
&= (1-p) \frac{1-s^3}{3} - (1-p)(1-s) \left( \frac{1+s}{2} \right)^2.
\end{align*}
$$

Finally,

$$
s_{\text{MIA,L}}^{*} = \underset{s \in [0,1]}{\operatorname{argmin}} \left\{ \left( \frac{p}{3} + \frac{(1-p)s^3}{3} \right) - \frac{1}{p + (1-p)s} \left( \frac{p}{2} + \frac{(1-p)s^2}{2} \right)^2 \right. \\
\left. + (1-p) \frac{1-s^3}{3} - (1-p)(1-s) \left( \frac{1+s}{2} \right)^2 \right\},
$$

which concludes the proof.
## B.7 Proof of proposition 2

**Probabilistic and block propagation.** First, note that the variable $X_2 = X_1 \mathbb{1}_{W=1}$ is similar to the variable studied for the computation of the MIA criterion in Proposition 1. Therefore, the value of the CART splitting criterion along the first variable is $C_{MIA}(1, 1/2, L, 0)$ and its value along the second variable is $C_{MIA}(2, s_{MIA,L}^*, L, \eta)$. Since the function

$$ \alpha \mapsto C_{MIA}(\cdot, s_{MIA,L}^*, L, \alpha) $$

is increasing, splitting along the first variable leads to the largest variance reduction. Thus, for probabilistic and block propagation, splits occur along the first variable. Let us now compare the value of these criteria. We have

$$ \mathbb{P}[X_1 \le 1/2] = \mathbb{P}[X_1 \ge 1/2] = 1/2. $$

The quantities related to the left cell are given by

$$ \mathbb{E}[Y|X_1 \le 1/2] = \frac{p+1}{4} \quad \text{and} \quad \mathbb{E}[Y^2|X_1 \le 1/2] = \frac{p}{4} + \frac{1}{12}. $$

The quantities related to the right cell are given by

$$ \mathbb{E}[Y|X_1 \ge 1/2] = \frac{3-p}{4} \quad \text{and} \quad \mathbb{E}[Y^2|X_1 \ge 1/2] = \frac{7}{12} - \frac{p}{4}. $$

Thus, the value of the criterion satisfies

$$ \begin{aligned} R(f_{\text{prob}}^*) &= \frac{1}{2}\left(\frac{p}{4} + \frac{1}{12} - \left(\frac{p+1}{4}\right)^2\right) \\ &\quad + \frac{1}{2}\left(\frac{7}{12} - \frac{p}{4} - \left(\frac{3-p}{4}\right)^2\right) \\ &= -\frac{p^2}{16} + \frac{p}{8} + \frac{1}{48}. \end{aligned} $$

Let, for all $p \in [0, 1]$,

$$ \begin{aligned} h(p) &= R(f_{\text{prob}}^*) - R(f_{\text{block}}^*) \\ &= -\frac{p^2}{16} + \frac{p}{8} + \frac{1}{48} - \left(-\frac{11}{48} + \frac{1}{8}\frac{3p+2}{2p+1}\right) \\ &= -\frac{p^2}{16} + \frac{p}{8} + \frac{1}{16} - \frac{1}{16}\frac{1}{2p+1}. \end{aligned} $$
## B.7 Proof of proposition 2

**Probabilistic and block propagation.** First, note that the variable $X_2 = X_1 \mathbb{1}_{W=1}$ is similar to the variable studied for the computation of the MIA criterion in Proposition 1. Therefore, the value of the CART splitting criterion along the first variable is $C_{MIA}(1, 1/2, L, 0)$ and its value along the second variable is $C_{MIA}(2, s_{MIA,L}^*, L, \eta)$. Since the function

$$ \alpha \mapsto C_{MIA}(\cdot, s_{MIA,L}^*, L, \alpha) $$

is increasing, splitting along the first variable leads to the largest variance reduction. Thus, for probabilistic and block propagation, splits occur along the first variable. Let us now compare the value of these criteria. We have

$$ \mathbb{P}[X_1 \le 1/2] = \mathbb{P}[X_1 \ge 1/2] = 1/2. $$

The quantities related to the left cell are given by

$$ \mathbb{E}[Y|X_1 \le 1/2] = \frac{p+1}{4} \quad \text{and} \quad \mathbb{E}[Y^2|X_1 \le 1/2] = \frac{p}{4} + \frac{1}{12}. $$

The quantities related to the right cell are given by

$$ \mathbb{E}[Y|X_1 \ge 1/2] = \frac{3-p}{4} \quad \text{and} \quad \mathbb{E}[Y^2|X_1 \ge 1/2] = \frac{7}{12} - \frac{p}{4}. $$

Thus, the value of the criterion satisfies

$$ \begin{aligned} R(f_{\text{prob}}^*) &= \frac{1}{2}\left(\frac{p}{4} + \frac{1}{12} - \left(\frac{p+1}{4}\right)^2\right) \\ &\quad + \frac{1}{2}\left(\frac{7}{12} - \frac{p}{4} - \left(\frac{3-p}{4}\right)^2\right) \\ &= -\frac{p^2}{16} + \frac{p}{8} + \frac{1}{48}. \end{aligned} $$

Let, for all $p \in [0, 1]$,

$$ \begin{aligned} h(p) &= R(f_{\text{prob}}^*) - R(f_{\text{block}}^*) \\ &= -\frac{p^2}{16} + \frac{p}{8} + \frac{1}{48} - \left(-\frac{11}{48} + \frac{1}{8}\frac{3p+2}{2p+1}\right) \\ &= -\frac{p^2}{16} + \frac{p}{8} + \frac{1}{16} - \frac{1}{16}\frac{1}{2p+1}. \end{aligned} $$
We have,

$$h'(p) = -\frac{p}{8} + \frac{1}{8} + \frac{1}{8(2p+1)^2},$$

and consequently,

$$h''(p) = -\frac{1}{8} - \frac{1}{2(2p+1)^3}.$$

An inspection of the variation of $h$ reveals that $h(p) \ge 0$ for all $p \in [0, 1]$, which concludes the first part of the proof.

**MIA.** As noticed above, the criterion computed along the second variable is given by

$$C_{\text{MIA}}(2, s_{\text{MIA},L}^*, L, \eta)$$

Since the function

$$\alpha \mapsto C_{\text{MIA}}(\cdot, s_{\text{MIA},L}^*, L, \alpha)$$

is increasing, MIA split will occur along the first variable if $p \le \eta$ and along the second variable if $p \ge \eta$. Therefore, the risk of the MIA splitting procedure is given by

$$R(f_{\text{MIA}}^*) = \min_{s \in [0,1]} C_{\text{MIA}}(1, s, L, p) \mathbf{1}_{p \le \eta} + \min_{s \in [0,1]} C_{\text{MIA}}(1, s, L, \eta) \mathbf{1}_{p > \eta}.$$

**Surrogate split.** Consider the model $Y = X_1$ and $X_2 = X_1 \mathbf{1}_{W=1}$, where $\mathbb{P}[W=0] = \eta$. Let us determine the best split along $X_2$ to predict $Z = \mathbf{1}_{X_1<0.5}$. Since $\{X_2 \le s\} = \{X_1 \le s, W=1\} \cup \{W=0\}$, and $\{X_2 > s\} = \{X_1 > s, W=1\}$,

$$\mathbb{P}[X_2 \le s] = s(1-p) + p \quad \text{and} \quad \mathbb{P}[X_2 > s] = (1-s)(1-p).$$

Consequently,

$$\begin{aligned} \mathbb{E}[Z | X_2 \le s] &= \frac{\mathbb{E}[\mathbf{1}_{X_1 \le 0.5, X_2 \le s}]}{\mathbb{P}[X_2 \le s]} \\ &= \frac{1}{s(1-p) + p} \mathbb{E}[\mathbf{1}_{X_1 \le 0.5, X_1 \le s, W=1} + \mathbf{1}_{X_1 \le 0.5, W=0}] \\ &= \frac{1}{s(1-p) + p} [(1-p) \min(0.5, s) + \frac{p}{2}]. \end{aligned}$$
We have,

$$h'(p) = -\frac{p}{8} + \frac{1}{8} + \frac{1}{8(2p+1)^2},$$

and consequently,

$$h''(p) = -\frac{1}{8} - \frac{1}{2(2p+1)^3}.$$

An inspection of the variation of $h$ reveals that $h(p) \ge 0$ for all $p \in [0, 1]$, which concludes the first part of the proof.

**MIA.** As noticed above, the criterion computed along the second variable is given by

$$C_{\text{MIA}}(2, s_{\text{MIA},L}^*, L, \eta)$$

Since the function

$$\alpha \mapsto C_{\text{MIA}}(\cdot, s_{\text{MIA},L}^*, L, \alpha)$$

is increasing, MIA split will occur along the first variable if $p \le \eta$ and along the second variable if $p \ge \eta$. Therefore, the risk of the MIA splitting procedure is given by

$$R(f_{\text{MIA}}^*) = \min_{s \in [0,1]} C_{\text{MIA}}(1, s, L, p) \mathbb{1}_{p \le \eta} + \min_{s \in [0,1]} C_{\text{MIA}}(1, s, L, \eta) \mathbb{1}_{p > \eta}.$$

**Surrogate split.** Consider the model $Y = X_1$ and $X_2 = X_1 \mathbb{1}_{W=1}$, where $\mathbb{P}[W=0] = \eta$. Let us determine the best split along $X_2$ to predict $Z = \mathbb{1}_{X_1<0.5}$. Since $\{X_2 \le s\} = \{X_1 \le s, W=1\} \cup \{W=0\}$, and $\{X_2 > s\} = \{X_1 > s, W=1\}$,

$$\mathbb{P}[X_2 \le s] = s(1-p) + p \quad \text{and} \quad \mathbb{P}[X_2 > s] = (1-s)(1-p).$$

Consequently,

$$\begin{aligned} \mathbb{E}[Z|X_2 \le s] &= \frac{\mathbb{E}[\mathbb{1}_{X_1 \le 0.5, X_2 \le s}]}{\mathbb{P}[X_2 \le s]} \\ &= \frac{1}{s(1-p) + p} \mathbb{E}[\mathbb{1}_{X_1 \le 0.5, X_1 \le s, W=1} + \mathbb{1}_{X_1 \le 0.5, W=0}] \\ &= \frac{1}{s(1-p) + p} [(1-p) \min(0.5, s) + \frac{p}{2}]. \end{aligned}$$
$$
\begin{align*}
\mathbb{E}[Z|X_2 \ge s] &= \frac{\mathbb{E}[1_{X_1 \le 0.5, X_2 > s}]}{\mathbb{P}[X_2 > s]} \\
&= \frac{1}{(1-s)(1-p)} \mathbb{P}[X_1 \le 0.5, W = 1, X_1 \ge s] \\
&= \frac{(0.5-s)_+}{1-s}.
\end{align*}
$$

Besides, note that $\mathbb{E}[Z^2] = \mathbb{P}[X_1 \le 0.5] = 0.5$. Therefore, the splitting criterion to predict $1_{X_1 \le 0.5}$ with $X_2$ is given by

$$
\begin{align*}
f(s) &= \frac{1}{2} - \mathbb{P}[X_2 \le s] (\mathbb{E}[Z|X_2 \le s])^2 - \mathbb{P}[X_2 > s] (\mathbb{E}[Z|X_2 > s])^2 \\
&= \frac{1}{2} - \frac{1}{s(1-p) + p} \left( (1-p) \min(0.5, s) + \frac{p}{2} \right)^2 - \frac{1-p}{1-s} ((0.5-s)_+)^2.
\end{align*}
$$

For $s \ge 1/2$,

$$
h(s) = \frac{1}{2} - \frac{1}{4(s(1-p) + p)},
$$

which is minimal for $s = 1/2$. For $s \le 1/2$,

$$
h(s) = \frac{1}{2} - \frac{1}{4} \left( \frac{p^2}{p + s(1-p)} + \frac{1-p}{1-s} \right).
$$

Hence,

$$
h'(s) = -\frac{1-p}{4} \frac{(1-2p)s^2 + 2ps}{(1-s)^2(s(1-p) + p)^2}.
$$

Let $g(s) = (1-2p)s^2 + 2ps$. If $p \le 1/2$, the solutions of $g(s) = 0$ are negative, thus, $g(s) \ge 0$ for all $s \in [0, 1/2]$ and thus the minimum of $h$ is reached at $s = 1/2$. If $p \ge 1/2$, one solution of $g(s) = 0$ is zero and the other is $s = 2p/(2p-1) > 1$. Thus, $g(s) \ge 0$ for all $s \in [0, 1/2]$ and the minimum of $h$ is reached at $s = 1/2$. Finally, the minimum of $h$ is reached at $s = 1/2$. The risk of the surrogate estimate is then given by

$$
\begin{align*}
R(f_{\text{surr}}^*) &= \mathbb{E}[(Y - f_{\text{surr}}^*(\mathbf{X}))^2] \\
&= \mathbb{E}[(Y - f_{\text{surr}}^*(\mathbf{X}))^2] \mathbb{1}_{M_1=0} + (Y - f_{\text{surr}}^*(\mathbf{X}))^2 \mathbb{1}_{M_1=1}.
\end{align*}
$$
$$ 
\begin{aligned}
\mathbb{E}[Z|X_2 \ge s] &= \frac{\mathbb{E}[1_{X_1 \le 0.5, X_2 > s}]}{\mathbb{P}[X_2 > s]} \\
&= \frac{1}{(1-s)(1-p)} \mathbb{P}[X_1 \le 0.5, W = 1, X_1 \ge s] \\
&= \frac{(0.5 - s)_+}{1 - s}.
\end{aligned}
$$

Besides, note that $\mathbb{E}[Z^2] = \mathbb{P}[X_1 \le 0.5] = 0.5$. Therefore, the splitting criterion to predict $1_{X_1 \le 0.5}$ with $X_2$ is given by

$$ 
\begin{aligned}
f(s) &= \frac{1}{2} - \mathbb{P}[X_2 \le s](\mathbb{E}[Z|X_2 \le s])^2 - \mathbb{P}[X_2 > s](\mathbb{E}[Z|X_2 > s])^2 \\
&= \frac{1}{2} - \frac{1}{s(1-p) + p}\left((1-p)\min(0.5, s) + \frac{p}{2}\right)^2 - \frac{1-p}{1-s}\left((0.5-s)_+\right)^2.
\end{aligned}
$$

For $s \ge 1/2$,

$$ h(s) = \frac{1}{2} - \frac{1}{4(s(1-p) + p)}, $$

which is minimal for $s = 1/2$. For $s \le 1/2$,

$$ h(s) = \frac{1}{2} - \frac{1}{4} \left( \frac{p^2}{p + s(1-p)} + \frac{1-p}{1-s} \right). $$

Hence,

$$ h'(s) = -\frac{1-p}{4} \frac{(1-2p)s^2 + 2ps}{(1-s)^2(s(1-p) + p)^2}. $$

Let $g(s) = (1-2p)s^2 + 2ps$. If $p \le 1/2$, the solutions of $g(s) = 0$ are negative, thus, $g(s) \ge 0$ for all $s \in [0, 1/2]$ and thus the minimum of $h$ is reached at $s = 1/2$. If $p \ge 1/2$, one solution of $g(s) = 0$ is zero and the other is $s = 2p/(2p-1) > 1$. Thus, $g(s) \ge 0$ for all $s \in [0, 1/2]$ and the minimum of $h$ is reached at $s = 1/2$. Finally, the minimum of $h$ is reached at $s = 1/2$. The risk of the surrogate estimate is then given by

$$ 
\begin{aligned}
R(f_{\text{surr}}^*) &= \mathbb{E}[(Y - f_{\text{surr}}^*(\mathbf{X}))^2] \\
&= \mathbb{E}[(Y - f_{\text{surr}}^*(\mathbf{X}))^2]1_{M_1=0} + (Y - f_{\text{surr}}^*(\mathbf{X}))^21_{M_1=1}.
\end{aligned}
$$
Here,

$$
\begin{align*}
& \mathbb{E}[(Y - f_{\text{surr}}^*(\mathbf{X}))^2 | M_1 = 1] \\
&= \mathbb{E}[(X_1 - 0.25)^2 \mathbf{1}_{X_2 < 0.5} + (X_1 - 0.75)^2 \mathbf{1}_{X_2 \ge 0.5}] \\
&= \eta \mathbb{E}[(X_1 - 0.25)^2] + (1 - \eta) \mathbb{E}[(X_1 - 0.25)^2 \mathbf{1}_{X_1 \le 0.5}] \\
&\quad + (1 - \eta) \mathbb{E}[(X_1 - 0.75)^2 \mathbf{1}_{X_1 > 0.5}] \\
&= \frac{1}{48} + \frac{6\eta}{48}.
\end{align*}
$$

Finally,

$$
R(f_{\text{surr}}^{*}) = \frac{1-p}{48} + p\left(\frac{1}{48} + \frac{6\eta}{48}\right) = \frac{1}{48} + \frac{6}{48}\eta p.
$$

C Miscellaneous

C.1 Variable selection properties of the tree methods with missing values

Decision trees based on the CART criterion (implemented in the R library rpart) and on
conditional trees (implemented in the R library partykit) lead to different ways of selecting
splitting variables. We illustrate this behaviour on the simple following model:

$$
\left\{
\begin{array}{l}
X_1 \perp X_2 \sim \mathcal{N}(0, 1) \\
\varepsilon \sim \mathcal{N}(0, 1) \\
Y = 0.25X_1 + \varepsilon.
\end{array}
\right.
$$

We insert MCAR values, either on the first variable or on both variables. Stumps (decision
trees of depth one) are fit on 500 Monte-Carlo repetitions. We vary the sample size and the
percentage of missing values. Figure 7 show that CART and conditional trees give similar
results when there are missing values on both variables. However, Figure 8 shows that CART
has a tendency to underselect X₁ when there are missing values only on X₁. For instance, for
a sample of size 50 with 75% missing values, CART selects the non-informative variable X₂
more frequently than X₁, while conditional trees keep selecting X₁ more often.
Here,

$$
\begin{align*}
& \mathbb{E}[(Y - f_{\text{surr}}^*(\mathbf{X}))^2 | M_1 = 1] \\
&= \mathbb{E}[(X_1 - 0.25)^2 \mathbf{1}_{X_2 < 0.5} + (X_1 - 0.75)^2 \mathbf{1}_{X_2 \ge 0.5}] \\
&= \eta \mathbb{E}[(X_1 - 0.25)^2] + (1 - \eta) \mathbb{E}[(X_1 - 0.25)^2 \mathbf{1}_{X_1 \le 0.5}] \\
&\quad + (1 - \eta) \mathbb{E}[(X_1 - 0.75)^2 \mathbf{1}_{X_1 > 0.5}] \\
&= \frac{1}{48} + \frac{6\eta}{48}.
\end{align*}
$$

Finally,

$$
R(f_{\text{surr}}^{*}) = \frac{1-p}{48} + p\left(\frac{1}{48} + \frac{6\eta}{48}\right) = \frac{1}{48} + \frac{6}{48}\eta p.
$$

C Miscellaneous

C.1 Variable selection properties of the tree methods with missing values

Decision trees based on the CART criterion (implemented in the R library rpart) and on
conditional trees (implemented in the R library partykit) lead to different ways of selecting
splitting variables. We illustrate this behaviour on the simple following model:

$$
\left\{
\begin{array}{l}
X_1 \perp X_2 \sim \mathcal{N}(0, 1) \\
\varepsilon \sim \mathcal{N}(0, 1) \\
Y = 0.25X_1 + \varepsilon.
\end{array}
\right.
$$

We insert MCAR values, either on the first variable or on both variables. Stumps (decision trees of depth one) are fit on 500 Monte-Carlo repetitions. We vary the sample size and the percentage of missing values. Figure 7 show that CART and conditional trees give similar results when there are missing values on both variables. However, Figure 8 shows that CART has a tendency to underselect X₁ when there are missing values only on X₁. For instance, for a sample of size 50 with 75% missing values, CART selects the non-informative variable X₂ more frequently than X₁, while conditional trees keep selecting X₁ more often.
<figure><img src="image_8.png" /><figcaption>Fig. 7: Frequency of selection of $X_1$ when there are missing values on $X_1$ and $X_2$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8: Frequency of selection of $X_1$ when there are missing values on $X_1$ only</figcaption></figure>

## C.2 Example of EM algorithm

Let us consider a simple case of $n$ observations $(\mathbf{x}_1, \mathbf{x}_2) = (x_{i1}, x_{i2})_{1 \le i \le n}$ sampled from the distribution of $(X_1, X_2)$, a bivariate Gaussian distribution with parameters $(\mu, \Sigma)$. We assume that $X_2$ is subjected to missing values and that only $r$ values are observed. The aim is to get the maximum likelihood estimates of $(\mu, \Sigma)$ from the incomplete data set. The algorithm described below can be straightforwardly extended to the multivariate case. Note that from $(\hat{\mu}, \hat{\Sigma})$, it is then possible to directly estimate the parameters of a linear regression model and thus to perform linear regression with missing values.

We denote by $f_{1,2}(\mathbf{x}_1, \mathbf{x}_2; \mu, \Sigma)$, $f_1(\mathbf{x}_1; \mu, \sigma_{11})$ and $f_{2|1}(\mathbf{x}_2|\mathbf{x}_1; \mu, \Sigma)$, respectively, the probability of joint distribution of $(X_1, X_2)$, marginal distribution of $X_1$ and conditional distribution of $X_2|X_1$. The joint distribution of observed data can be decomposed as:

$$f_{1,2}(\mathbf{x}_1, \mathbf{x}_2; \mu, \Sigma) = \prod_{i=1}^{n} f_1(x_{i1}; \mu_1, \sigma_{11}) \prod_{j=1}^{r} f_{2|1}(x_{j2}|x_{j1}; \mu, \Sigma),$$
<figure><img src="image_9.png" /><figcaption>Fig. 7: Frequency of selection of $X_1$ when there are missing values on $X_1$ and $X_2$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8: Frequency of selection of $X_1$ when there are missing values on $X_1$ only</figcaption></figure>

## C.2 Example of EM algorithm

Let us consider a simple case of $n$ observations $(\mathbf{x}_1, \mathbf{x}_2) = (x_{i1}, x_{i2})_{1 \le i \le n}$ sampled from the distribution of $(X_1, X_2)$, a bivariate Gaussian distribution with parameters $(\mu, \Sigma)$. We assume that $X_2$ is subjected to missing values and that only $r$ values are observed. The aim is to get the maximum likelihood estimates of $(\mu, \Sigma)$ from the incomplete data set. The algorithm described below can be straightforwardly extended to the multivariate case. Note that from $(\hat{\mu}, \hat{\Sigma})$, it is then possible to directly estimate the parameters of a linear regression model and thus to perform linear regression with missing values.

We denote by $f_{1,2}(\mathbf{x}_1, \mathbf{x}_2; \mu, \Sigma)$, $f_1(\mathbf{x}_1; \mu_1, \sigma_{11})$ and $f_{2|1}(\mathbf{x}_2|\mathbf{x}_1; \mu, \Sigma)$, respectively, the probability of joint distribution of $(X_1, X_2)$, marginal distribution of $X_1$ and conditional distribution of $X_2|X_1$. The joint distribution of observed data can be decomposed as:

$$f_{1,2}(\mathbf{x}_1, \mathbf{x}_2; \mu, \Sigma) = \prod_{i=1}^{n} f_1(x_{i1}; \mu_1, \sigma_{11}) \prod_{j=1}^{r} f_{2|1}(x_{j2}|x_{j1}; \mu, \Sigma),$$
and the observed log-likelihood is written (up to an additional constant that does not appear in the maximization and that we therefore drop):

$$ \ell(\mu, \Sigma; \mathbf{x}_1, \mathbf{x}_2) = -\frac{n}{2} \log(\sigma_{11}^2) - \frac{1}{2} \sum_{i=1}^{n} \frac{(x_{i1} - \mu_1)^2}{\sigma_{11}^2} - \frac{r}{2} \log \left( (\sigma_{22} - \frac{\sigma_{12}^2}{\sigma_{11}^2})^2 \right) \\ -\frac{1}{2} \sum_{i=1}^{r} \frac{\left( x_{i2} - \mu_2 - \frac{\sigma_{12}}{\sigma_{11}}(x_{i1} - \mu_1) \right)^2}{(\sigma_{22} - \frac{\sigma_{12}^2}{\sigma_{11}^2})^2} $$

We skip the computations and directly give the expression of the closed form maximum likelihood estimates of the mean:

$$ \hat{\mu}_1 = n^{-1} \sum_{i=1}^{n} x_{i1} \\ \hat{\mu}_2 = \hat{\beta}_{20.1} + \hat{\beta}_{21.1} \hat{\mu}_1, $$

where

$$ \hat{\beta}_{21.1} = s_{12}/s_{11}, \quad \hat{\beta}_{20.1} = \bar{x}_2 - \hat{\beta}_{21.1}\bar{x}_1, \\ \bar{x}_j = r^{-1} \sum_{i=1}^{r} x_{ij} \quad \text{and} \quad s_{jk} = r^{-1} \sum_{i=1}^{r} (x_{ij} - \bar{x}_j)(x_{ik} - \bar{x}_k), \quad j, k = 1, 2. $$

In this simple setting, we have an explicit expression of the maximum likelihood estimator despite missing values. However, this is not always the case but it is possible to use an EM algorithm to get the maximum likelihood estimators in the cases where data are missing.

The EM algorithm consists in maximizing the observed likelihood through successive maximization of the complete likelihood (if we had observed all n realizations of x₁ and x₂). Maximizing the complete likelihood

$$ \ell_c(\mu, \Sigma; \mathbf{x}_1, \mathbf{x}_2) = -\frac{n}{2} \log(\det(\Sigma)) - \frac{1}{2} \sum_{i=1}^{n} (x_{i1} - \mu_1)^T \Sigma^{-1} (x_{i1} - \mu_1) $$

would be straightforward if we had all the observations. However elements of this likelihood are not available. Therefore, we replace them by the conditional expectation given observed data and the parameters of the current iteration. These two steps of computation of the conditional expectation (E-step) and maximization of the completed likelihood (M step) are repeated until convergence. The update formulas for the E and M steps are as follows:

**E step:** The sufficient statistics of the likelihood are:

$$ s_1 = \sum_{i=1}^{n} x_{i1}, \quad s_2 = \sum_{i=1}^{n} x_{i2}, \quad s_{11} = \sum_{i=1}^{n} x_{i1}^2, \quad s_{22} = \sum_{i=1}^{n} x_{i2}^2, \quad s_{12} = \sum_{i=1}^{n} x_{i1}x_{i2}. $$
and the observed log-likelihood is written (up to an additional constant that does not appear in the maximization and that we therefore drop):

$$ \ell(\mu, \Sigma; \mathbf{x}_1, \mathbf{x}_2) = -\frac{n}{2} \log(\sigma_{11}^2) - \frac{1}{2} \sum_{i=1}^{n} \frac{(x_{i1} - \mu_1)^2}{\sigma_{11}^2} - \frac{r}{2} \log \left( (\sigma_{22} - \frac{\sigma_{12}^2}{\sigma_{11}^2})^2 \right) \\ -\frac{1}{2} \sum_{i=1}^{r} \frac{(x_{i2} - \mu_2 - \frac{\sigma_{12}}{\sigma_{11}}(x_{i1} - \mu_1))^2}{(\sigma_{22} - \frac{\sigma_{12}^2}{\sigma_{11}^2})^2} $$

We skip the computations and directly give the expression of the closed form maximum likelihood estimates of the mean:

$$ \hat{\mu}_1 = n^{-1} \sum_{i=1}^{n} x_{i1} \\ \hat{\mu}_2 = \hat{\beta}_{20.1} + \hat{\beta}_{21.1} \hat{\mu}_1, $$

where

$$ \hat{\beta}_{21.1} = s_{12}/s_{11}, \quad \hat{\beta}_{20.1} = \bar{x}_2 - \hat{\beta}_{21.1}\bar{x}_1, \\ \bar{x}_j = r^{-1} \sum_{i=1}^{r} x_{ij} \quad \text{and} \quad s_{jk} = r^{-1} \sum_{i=1}^{r} (x_{ij} - \bar{x}_j)(x_{ik} - \bar{x}_k), \quad j, k = 1, 2. $$

In this simple setting, we have an explicit expression of the maximum likelihood estimator despite missing values. However, this is not always the case but it is possible to use an EM algorithm to get the maximum likelihood estimators in the cases where data are missing.

The EM algorithm consists in maximizing the observed likelihood through successive maximization of the complete likelihood (if we had observed all n realizations of x₁ and x₂). Maximizing the complete likelihood

$$ \ell_c(\mu, \Sigma; \mathbf{x}_1, \mathbf{x}_2) = -\frac{n}{2} \log(\det(\Sigma)) - \frac{1}{2} \sum_{i=1}^{n} (x_{i1} - \mu_1)^T \Sigma^{-1} (x_{i1} - \mu_1) $$

would be straightforward if we had all the observations. However elements of this likelihood are not available. Therefore, we replace them by the conditional expectation given observed data and the parameters of the current iteration. These two steps of computation of the conditional expectation (E-step) and maximization of the completed likelihood (M step) are repeated until convergence. The update formulas for the E and M steps are as follows:

**E step:** The sufficient statistics of the likelihood are:

$$ s_1 = \sum_{i=1}^{n} x_{i1}, \quad s_2 = \sum_{i=1}^{n} x_{i2}, \quad s_{11} = \sum_{i=1}^{n} x_{i1}^2, \quad s_{22} = \sum_{i=1}^{n} x_{i2}^2, \quad s_{12} = \sum_{i=1}^{n} x_{i1}x_{i2}. $$
Since some values of $x_2$ are not available, we fill in the sufficient statistics with:

$$ \mathbb{E}[x_{i2}|x_{i1}; \mu, \Sigma] = \beta_{20.1} + \beta_{21.1}x_{i1} $$

$$ \mathbb{E}[x_{i2}^2|x_{i1}; \mu, \Sigma] = (\beta_{20.1} + \beta_{21.1}x_{i1})^2 + \sigma_{22.1} $$

$$ \mathbb{E}[x_{i2}x_{i2}|x_{i1}; \mu, \Sigma] = (\beta_{20.1} + \beta_{21.1}x_{i1})x_{i1}. $$

with, $\beta_{21.1} = \sigma_{12}/\sigma_{11}$, $\beta_{20.1} = \mu_2 - \beta_{21.1}\mu_1$, and $\sigma_{22.1} = \sigma_{22} - \sigma_{12}^2/\sigma_{11}$.

**M step:** The M step consists in computing the maximum likelihood estimates as usual. Given $s_1, s_2, s_{11}, s_{22}$, and $s_{12}$, update $\hat{\mu}$ and $\hat{\sigma}$ with

$$ \hat{\mu}_1 = s_1/n, \quad \hat{\mu}_2 = s_2/n, $$

$$ \hat{\sigma}_1 = s_{11}/n - \hat{\mu}_1^2, \quad \hat{\sigma}_2 = s_{22}/n - \hat{\mu}_2^2, \quad \hat{\sigma}_{12} = s_{12}/n - \hat{\mu}_1 \hat{\mu}_2 $$

Note that $s_1, s_{11}, \hat{\mu}_1$ and $\hat{\sigma}_1$ are constant across iterations since we do not have missing values on $x_1$.

*Remark 4* Note that EM imputes the sufficient statistics and not the data.
Since some values of $x_2$ are not available, we fill in the sufficient statistics with:

$$ \mathbb{E}[x_{i2}|x_{i1}; \mu, \Sigma] = \beta_{20.1} + \beta_{21.1}x_{i1} $$

$$ \mathbb{E}[x_{i2}^2|x_{i1}; \mu, \Sigma] = (\beta_{20.1} + \beta_{21.1}x_{i1})^2 + \sigma_{22.1} $$

$$ \mathbb{E}[x_{i2}x_{i2}|x_{i1}; \mu, \Sigma] = (\beta_{20.1} + \beta_{21.1}x_{i1})x_{i1}. $$

with, $\beta_{21.1} = \sigma_{12}/\sigma_{11}$, $\beta_{20.1} = \mu_2 - \beta_{21.1}\mu_1$, and $\sigma_{22.1} = \sigma_{22} - \sigma_{12}^2/\sigma_{11}$.

**M step:** The M step consists in computing the maximum likelihood estimates as usual. Given $s_1, s_2, s_{11}, s_{22}$, and $s_{12}$, update $\hat{\mu}$ and $\hat{\sigma}$ with

$$ \hat{\mu}_1 = s_1/n, \quad \hat{\mu}_2 = s_2/n, $$

$$ \hat{\sigma}_1 = s_{11}/n - \hat{\mu}_1^2, \quad \hat{\sigma}_2 = s_{22}/n - \hat{\mu}_2^2, \quad \hat{\sigma}_{12} = s_{12}/n - \hat{\mu}_1 \hat{\mu}_2 $$

Note that $s_1, s_{11}, \hat{\mu}_1$ and $\hat{\sigma}_1$ are constant across iterations since we do not have missing values on $x_1$.

*Remark 4* Note that EM imputes the sufficient statistics and not the data.
D Additional experiments

D.1 Varying the correlation strength and the missing rate

<figure><img src="image_10.png" /><figcaption>Fig. 9: *R<sup>2</sup>* scores on model 1 • Normalized explained variance for different mechanisms with 20% of missing values, *n* = 1000, *d* = 9 and *ρ* = 0.2.</figcaption></figure>
D Additional experiments

D.1 Varying the correlation strength and the missing rate

<figure><img src="image_11.png" /><figcaption>Fig. 9: *R<sup>2</sup>* scores on model 1 • Normalized explained variance for different mechanisms with 20% of missing values, *n* = 1000, *d* = 9 and *ρ* = 0.2.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 10: *R² scores on model 1* • Normalized explained variance for different mechanisms with 20% of missing values, *n* = 1000, *d* = 9 and *ρ* = 0.5.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 10: *R² scores on model 1* • Normalized explained variance for different mechanisms with 20% of missing values, *n* = 1000, *d* = 9 and *ρ* = 0.5.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 11: R² scores on model 1 • Normalized explained variance for different mechanisms with 20% of missing values, n = 1000, d = 9 and ρ = 0.8.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 11: R² scores on model 1 • Normalized explained variance for different mechanisms with 20% of missing values, n = 1000, d = 9 and ρ = 0.8.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 12: R² scores on model 1 • Normalized explained variance for different mechanisms with 40% of missing values, n = 1000, d = 9 and ρ = 0.5.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 12: R² scores on model 1 • Normalized explained variance for different mechanisms with 40% of missing values, n = 1000, d = 9 and ρ = 0.5.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 13: R² scores on model 1 • Normalized explained variance for different mechanisms with 60% of missing values, n = 1000, d = 9 and ρ = 0.5.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 13: R² scores on model 1 • Normalized explained variance for different mechanisms with 60% of missing values, n = 1000, d = 9 and ρ = 0.5.</figcaption></figure>
D.2 Statistical tests

<figure><img src="image_15.png" /><figcaption>Fig. 14: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MCAR model ($\rho$ = 0.2, missing rate of 20%)</figcaption></figure>
D.2 Statistical tests

<figure><img src="image_16.png" /><figcaption>Fig. 14: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MCAR model ($\rho$ = 0.2, missing rate of 20%)</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 15: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MNAR model ($\rho = 0.2$, missing rate of 20%)</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 15: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MNAR model ($\rho = 0.2$, missing rate of 20%)</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 16: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a Predictive M model ($\rho = 0.2$, missing rate of 20%)</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 16: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a Predictive M model ($\rho = 0.2$, missing rate of 20%)</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 17: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MCAR model ($\rho = 0.5$, missing rate of 20%)</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 17: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MCAR model ($\rho = 0.5$, missing rate of 20%)</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 18: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MNAR model ($\rho = 0.5$, missing rate of 20%)</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 18: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MNAR model ($\rho = 0.5$, missing rate of 20%)</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 19: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a Predictive M model ($\rho = 0.5$, missing rate of 20%)</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Fig. 19: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a Predictive M model ($\rho = 0.5$, missing rate of 20%)</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Fig. 20: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MCAR model ($\rho = 0.8$, missing rate of 20%)</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. 20: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MCAR model ($\rho = 0.8$, missing rate of 20%)</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. 21: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MNAR model ($\rho = 0.8$, missing rate of 20%)</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. 21: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MNAR model ($\rho = 0.8$, missing rate of 20%)</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. 22: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a Predictive M model ($\rho = 0.8$, missing rate of 20%)</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Fig. 22: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a Predictive M model ($\rho = 0.8$, missing rate of 20%)</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Fig. 23: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MCAR model ($\rho = 0.5$, missing rate of 40%)</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Fig. 23: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MCAR model ($\rho = 0.5$, missing rate of 40%)</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Fig. 24: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MNAR model ($\rho = 0.5$, missing rate of 40%)</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Fig. 24: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MNAR model ($\rho = 0.5$, missing rate of 40%)</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Fig. 25: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a Predictive M model ($\rho = 0.5$, missing rate of 40%)</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Fig. 25: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a Predictive M model ($\rho = 0.5$, missing rate of 40%)</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Fig. 26: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MCAR model ($\rho = 0.5$, missing rate of 60%)</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Fig. 26: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MCAR model ($\rho = 0.5$, missing rate of 60%)</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Fig. 27: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MNAR model ($\rho = 0.5$, missing rate of 60%)</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>Fig. 27: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a MNAR model ($\rho = 0.5$, missing rate of 60%)</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>Fig. 28: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a Predictive M model ($\rho = 0.5$, missing rate of 60%)</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Fig. 28: P-values of paired t-tests for each pair of imputation method and for each learning algorithm in a Predictive M model ($\rho = 0.5$, missing rate of 60%)</figcaption></figure>