# A second look at the Kurth solution in galactic dynamics

Markus Kunze1
1 Universität Köln, Institut für Mathematik, Weyertal 86-90,
D - 50931 Köln, Germany
email<EMAIL_ADDRESS>
Key words: Vlasov-Poisson system, Kurth solution, galactic dynamics

###### Abstract

The Kurth solution is a particular non-isotropic steady state solution to the
gravitational Vlasov-Poisson system. It has the property that by means of a
suitable time-dependent transformation it can be turned into a family of time-
dependent solutions. Therefore, for a general steady state
$Q(x,v)=\tilde{Q}(e_{Q},\beta)$, depending upon the particle energy $e_{Q}$
and $\beta=\ell^{2}=|x\wedge v|^{2}$, the question arises if solutions $f$
could be generated that are of the form

$f(t)=\tilde{Q}\Big{(}e_{Q}(R(t),P(t),B(t)),B(t)\Big{)}$

for suitable functions $R$, $P$ and $B$, all depending on $(t,r,p_{r},\beta)$
for $r=|x|$ and $p_{r}=\frac{x\cdot v}{|x|}$. We are going to show that, under
some mild assumptions, basically if $R$ and $P$ are independent of $\beta$,
and if $B=\beta$ is constant, then $Q$ already has to be the Kurth solution.

This paper is dedicated to the memory of Professor Robert Glassey.

##  1 Introduction

It is a remarkable fact that very few of Bob Glassey’s influential papers
concern the Vlasov-Poisson system. Certainly there is an in-depth treatment of
the existence of global solutions in his book [3], but apart from that only
[4, 5] comes to this author’s mind. Maybe this is due to Glassey’s
mathematical formation in the tradition of John, Nirenberg, Segal, Strauss …
that he liked better hyperbolic equations, and in particular the relativistic
Vlasov-Maxwell system.

For this reason, Glassey would have probably not paid much attention to the
present paper, but being a polite person, he would nevertheless have found
some friendly words for it. In addition, this paper has no hard analytic
proofs, which Glassey could do so well. Let us only mention [6] on global
existence for the ‘$2.5$’ dimensional relativistic Vlasov-Maxwell system,
jointly with Jack Schaeffer, which is not so well-known (in the sense that not
many people have read it in all detail), but which is a true masterpiece. One
has to use all kinds of structures in the system and is not allowed to loose
the tiniest part of an $\varepsilon$ to close the argument in the end.

Here we are going to consider the Vlasov-Poisson system in the gravitational
case, which is given by

$\partial_{t}f(t,x,v)+v\cdot\nabla_{x}f(t,x,v)-\nabla_{x}U_{f}(t,x)\cdot\nabla_{v}f(t,x,v)=0,$
(1.1)

where

$\Delta_{x}U_{f}(t,x)=4\pi\rho_{f}(t,x),\quad\lim_{|x|\to\infty}U_{f}(t,x)=0,\quad\rho_{f}(t,x)=\int_{\mathbb{R}^{3}}f(t,x,v)\,dv,$
(1.2)

for $(t,x,v)\in\mathbb{R}\times\mathbb{R}^{3}\times\mathbb{R}^{3}$. Therefore

$U_{f}(t,x)=-\int_{\mathbb{R}^{3}}\frac{\rho_{f}(t,y)}{|y-x|}\,dy.$ (1.3)

The system possesses an abundance of solutions $Q=Q(x,v)$ that are independent
of time. Let

$e_{Q}(x,v)=\frac{1}{2}\,|v|^{2}+U_{Q}(x)$

denote the particle energy and let

$\ell^{2}(x,v)=|L|^{2}=|x|^{2}|v|^{2}-(x\cdot v)^{2}$

be the square of the angular momentum $L=x\wedge v$, respectively. Then both
$e_{Q}$ and $\ell^{2}$ are conserved along solutions of the characteristic
equations $\ddot{X}(s)=-\nabla U_{Q}(X(s))$; note that also $U_{Q}$ is
independent of time. Next recall that a function $g=g(x,v)$ is said to be
spherically symmetric, if $g(Ax,Av)=g(x,v)$ for all $A\in{\rm SO}(3)$ and
$x,v\in\mathbb{R}^{3}$. Now it is the content of Jeans’s theorem that the
distribution function $Q$ of every spherically symmetric steady state solution
has to be of the form

$Q(x,v)=\tilde{Q}(e_{Q}(x,v),\ell^{2}(x,v))$

for a suitable function $\tilde{Q}$ of two variables; see [1, Section 2] for a
precise formulation. Such steady state solutions are called non-isotropic, in
contrast to the isotropic ones, which can be written as
$Q(x,v)=\tilde{Q}(e_{Q}(x,v))$; a solution of the latter form will necessarily
be spherically symmetric, [2, 11].

In this paper we will have a closer look at one particular and non-isotropic
steady state solution $Q$, which has been found by Kurth in 1978 and which
will be denoted by $Q_{{\rm K}}$ in the sequel; see [7]. It is surrounded by
time-periodic solutions $f_{\varepsilon}(t)$ such that $f_{\varepsilon}\to
Q_{{\rm K}}$ as $\varepsilon\to 0$. Since the $f_{\varepsilon}(t)$ are semi-
explicit, the Kurth solution is a good testing ground for all kinds of
questions, including some from numerics [10]. It is very degenerate in many
respects, so an important issue is to understand whether it reflects what
happens ‘generically’ close to steady states (in a sense to be made precise),
or on the contrary it is just a peculiarity. The Kurth solution is given by

$\displaystyle Q_{{\rm
K}}(x,v)=\frac{3}{4\pi^{3}}\,\frac{1}{(1-|x|^{2}-|v|^{2}+|x\wedge
v|^{2})^{1/2}}\quad\mbox{where}\quad(\ldots)>0\,\,\mbox{and}\,\,|x\wedge
v|<1,$ $\displaystyle\mbox{and}\quad Q_{{\rm K}}(x,v)=0\quad\mbox{else},$
(1.4)

for $x,v\in\mathbb{R}^{3}$. Then (see Lemma 3.1 below) its charge density

$\rho_{Q_{{\rm K}}}(x)=\int_{\mathbb{R}^{3}}Q_{{\rm
K}}(x,v)\,dv=\frac{3}{4\pi}\,{\bf 1}_{B_{1}(0)}(x)$

is, up to a factor, the characteristic function of the unit ball in
$\mathbb{R}^{3}$. The solution to $\Delta U_{Q_{{\rm K}}}=4\pi\rho_{Q_{{\rm
K}}}$ and $U_{Q_{{\rm K}}}(x)\to 0$ as $|x|\to\infty$ is given by

$U_{Q_{{\rm
K}}}(x)=\left\\{\begin{array}[]{c@{\quad:\quad}c}\frac{1}{2}\,|x|^{2}-\frac{3}{2}&|x|\leq
1\\\\[4.30554pt] -\frac{1}{|x|}&|x|>1\end{array}\right..$ (1.5)

Next consider the second-order ODE
$\ddot{\phi}=-\frac{1}{\phi^{2}}+\frac{1}{\phi^{3}}$ and let
$\phi_{\varepsilon}$ denote the solution such that $\phi_{\varepsilon}(0)=1$
and $\dot{\phi}_{\varepsilon}(0)=\varepsilon$. It follows that
$\phi_{\varepsilon}$ is periodic for $\varepsilon<1$ (in fact
$|\varepsilon|<1$), and its period is calculated to be
$T_{\varepsilon}=\frac{2\pi}{(1-\varepsilon^{2})^{3/2}}$. Defining

$f_{\varepsilon}(t,x,v)=Q_{{\rm
K}}\Big{(}\frac{x}{\phi_{\varepsilon}(t)},\phi_{\varepsilon}(t)v-\dot{\phi}_{\varepsilon}(t)x\Big{)},\quad
t\in\mathbb{R},\quad x,v\in\mathbb{R}^{3},$ (1.6)

the $f_{\varepsilon}(t)$ are $T_{\varepsilon}$-periodic and (formal) solutions
to the gravitational Vlasov-Poisson system; see Lemma 3.2 below. We may also
write

$f_{\varepsilon}(t)=Q_{{\rm
K}}\circ\Lambda_{\varepsilon}(t),\quad\Lambda_{\varepsilon}(t)(x,v)=\Big{(}\frac{x}{\phi_{\varepsilon}(t)},\phi_{\varepsilon}(t)v-\dot{\phi}_{\varepsilon}(t)x\Big{)}.$
(1.7)

The associated density is

$\rho_{\varepsilon}(t,x)=\int_{\mathbb{R}^{3}}f_{\varepsilon}(t,x,v)\,dv=\frac{3}{4\pi}\,\frac{1}{\phi_{\varepsilon}(t)^{3}}\,{\bf
1}_{\\{|x|<\phi_{\varepsilon}(t)\\}}=\frac{1}{\phi_{\varepsilon}(t)^{3}}\,\rho_{Q_{{\rm
K}}}\Big{(}\frac{x}{\phi_{\varepsilon}(t)}\Big{)},$

resulting in the potential

$U_{\varepsilon}(t,x)=\frac{1}{\phi_{\varepsilon}(t)}\,U_{Q_{{\rm
K}}}\Big{(}\frac{x}{\phi_{\varepsilon}(t)}\Big{)}.$ (1.8)

The function $Q_{{\rm K}}$ is spherically symmetric, hence so is
$f_{\varepsilon}(t)$, since

$f_{\varepsilon}(t,Ax,Av)=Q_{{\rm
K}}\Big{(}A\frac{x}{\phi_{\varepsilon}(t)},A[\phi_{\varepsilon}(t)v-\dot{\phi}_{\varepsilon}(t)x]\Big{)}=Q_{{\rm
K}}\Big{(}\frac{x}{\phi_{\varepsilon}(t)},\phi_{\varepsilon}(t)v-\dot{\phi}_{\varepsilon}(t)x\Big{)}=f_{\varepsilon}(t,x,v)$

for $t\in\mathbb{R}$, $A\in{\rm SO}(3)$, $x,v\in\mathbb{R}^{3}$. Therefore we
may re-express everything in the adapted spherically symmetric variables

$r=|x|,\quad p_{r}=\frac{x\cdot v}{|x|},\quad\ell=|x\wedge v|.$

To begin with,

$\displaystyle e_{Q_{K}}(r,p_{r},\ell^{2})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,|v|^{2}+U_{Q_{{\rm
K}}}(r)=\frac{1}{2}\,p_{r}^{2}+U_{{\rm eff,K}}(r,\ell^{2}),$ (1.9)
$\displaystyle\quad U_{{\rm eff,K}}(r,\ell^{2})$ $\displaystyle=$
$\displaystyle U_{Q_{{\rm
K}}}(r)+\frac{\ell^{2}}{2r^{2}}=\left\\{\begin{array}[]{c@{\quad:\quad}c}\frac{r^{2}}{2}-\frac{3}{2}+\frac{\ell^{2}}{2r^{2}}&r\leq
1\\\\[4.30554pt] -\frac{1}{r}+\frac{\ell^{2}}{2r^{2}}&r\geq
1\end{array}\right.,$ (1.12)

by (1.5). Here $U_{{\rm eff,K}}$ is called the effective potential, and
henceforth we will sometimes write $e$ instead of $e_{Q_{{\rm K}}}$. Also,

$1-|x|^{2}-|v|^{2}+|x\wedge
v|^{2}=1-r^{2}-p_{r}^{2}-\frac{\ell^{2}}{r^{2}}+\ell^{2}=-2(1+e)+\ell^{2},$
(1.13)

so that

$\displaystyle\tilde{Q}_{{\rm
K}}(e,\ell^{2})=\frac{3}{4\pi^{3}}\,\frac{1}{(-2(1+e)+\ell^{2})^{1/2}}\quad\mbox{where}\quad(\ldots)>0\,\,\mbox{and}\,\,\ell<1,$
$\displaystyle\mbox{and}\quad\tilde{Q}_{{\rm
K}}(e,\ell^{2})=0\quad\mbox{else}.$ (1.14)

In spherically symmetric coordinates, the $\Lambda_{\varepsilon}(t)$ from
(1.7) are identified with

$\Lambda_{\varepsilon}(t)(r,p_{r})=\Big{(}\frac{r}{\phi_{\varepsilon}(t)},\phi_{\varepsilon}(t)p_{r}-\dot{\phi}_{\varepsilon}(t)r\Big{)},$
(1.15)

since $f_{\varepsilon}(t)=Q_{{\rm K}}\circ\Lambda_{\varepsilon}(t)$ and

$\displaystyle
1-\Big{|}\frac{x}{\phi_{\varepsilon}(t)}\Big{|}^{2}-|\phi_{\varepsilon}(t)v-\dot{\phi}_{\varepsilon}(t)x|^{2}+\Big{|}\frac{x}{\phi_{\varepsilon}(t)}\wedge[\phi_{\varepsilon}(t)v-\dot{\phi}_{\varepsilon}(t)x]\Big{|}^{2}$
$\displaystyle=$ $\displaystyle
1-\frac{r^{2}}{\phi_{\varepsilon}(t)^{2}}-\phi_{\varepsilon}(t)^{2}\,\Big{(}p_{r}^{2}+\frac{\ell^{2}}{r^{2}}\Big{)}+2\phi_{\varepsilon}(t)\dot{\phi}_{\varepsilon}(t)\,rp_{r}-\dot{\phi}_{\varepsilon}(t)^{2}r^{2}+\ell^{2}$
$\displaystyle=$ $\displaystyle
1-\frac{r^{2}}{\phi_{\varepsilon}(t)^{2}}-(\phi_{\varepsilon}(t)p_{r}-\dot{\phi}_{\varepsilon}(t)r)^{2}-\phi_{\varepsilon}(t)^{2}\frac{\ell^{2}}{r^{2}}+\ell^{2}=F\Big{(}\frac{r}{\phi_{\varepsilon}(t)},\phi_{\varepsilon}(t)p_{r}-\dot{\phi}_{\varepsilon}(t)r,\ell^{2}\Big{)}$

for

$F(r,p_{r},\ell^{2})=1-r^{2}-p_{r}^{2}-\frac{\ell^{2}}{r^{2}}+\ell^{2}.$
(1.16)

Now we are in position to describe the main result of this paper. Writing
$\beta=\ell^{2}$, the Kurth solution is

$Q_{{\rm K}}=\tilde{Q}_{{\rm K}}(e_{Q},\beta)=\tilde{Q}_{{\rm
K}}\Big{(}e_{Q}(r,p_{r},\beta),\beta\Big{)},$

whereas the neighboring $f_{\varepsilon}(t)=f_{\varepsilon}(t,r,p_{r},\beta)$
can be expressed as

$f_{\varepsilon}(t)=\tilde{Q}_{{\rm
K}}\Big{(}e_{Q}(R_{\varepsilon}(t),P_{\varepsilon}(t),B_{\varepsilon}(t)),B_{\varepsilon}(t)\Big{)}$

for

$\displaystyle R_{\varepsilon}(t)$ $\displaystyle=$ $\displaystyle
R_{\varepsilon}(t,r,p_{r},\beta)=\frac{r}{\phi_{\varepsilon}(t)},$ (1.17)
$\displaystyle P_{\varepsilon}(t)$ $\displaystyle=$ $\displaystyle
P_{\varepsilon}(t,r,p_{r},\beta)=\phi_{\varepsilon}(t)p_{r}-\dot{\phi}_{\varepsilon}(t)r,$
(1.18) $\displaystyle B_{\varepsilon}(t)$ $\displaystyle=$ $\displaystyle
B_{\varepsilon}(t,r,p_{r},\beta)=\beta,$ (1.19)

according to (1.15), (1.13) and (1.16). It should be remarked (as is verified
in Lemma 3.3 below) that defining

$H_{\varepsilon,{\rm
K}}(t,r,p_{r},\beta)=-\frac{\dot{\phi}_{\varepsilon}(t)}{\phi_{\varepsilon}(t)}\,r\,p_{r}-\frac{1}{2}\,(\dot{\phi}_{\varepsilon}(t)^{2}-\phi_{\varepsilon}(t)\ddot{\phi}_{\varepsilon}(t))\,r^{2},$

then

$\frac{d}{dt}\,Z_{\varepsilon}(t)=J\nabla H_{\varepsilon,{\rm
K}}(t,Z_{\varepsilon}(t),B_{\varepsilon}(t)),$ (1.20)

where

$Z_{\varepsilon}=(R_{\varepsilon},P_{\varepsilon}),\quad
J=\left(\begin{array}[]{cc}0&1\\\ -1&0\end{array}\right).$

In other words, the time evolution of $Z_{\varepsilon}$ is governed by the
time-dependent Hamiltonian $H_{\varepsilon,{\rm K}}$. Also note that both
$R_{\varepsilon}$ and $P_{\varepsilon}$ are in fact independent of $\beta$,
and $B_{\varepsilon}=\beta$ is constant.

Thus, for a general steady state $Q(x,v)=\tilde{Q}(e_{Q},\beta)$, the question
arises if solutions $f$ could be found that are of the form

$f(t)=\tilde{Q}\Big{(}e_{Q}(R(t),P(t),B(t)),B(t)\Big{)}$ (1.21)

for suitable functions $R$, $P$ and $B$, all depending on $(t,r,p_{r},\beta)$,
such that the evolution of $Z=(R,P)$ is Hamiltonian. For the moment it will
play no role if the $f$ come in a family of $f_{\varepsilon}$ that is close to
$Q$ as $\varepsilon\to 0$, or if the function(s) are periodic or not.

We are going to show that, basically, if $R$ and $P$ are independent of
$\beta$, and if $B=\beta$ is constant, then $Q$ already has to be the Kurth
solution $Q_{{\rm K}}$.

###### Theorem 1.1

Suppose that the functions

$(R,P,B)(t)=(R(t,r,p_{r},\beta),P(t,r,p_{r},\beta),B(t,r,p_{r},\beta))$

are such that $f(t)=f(t,r,p_{r},\beta)$ is a solution to the gravitational
Vlasov-Poisson system and moreover $R$ and $P$ are independent of $\beta$, and
$B=\beta$ is constant:

$(R,P,B)(t)=(R(t,r,p_{r}),P(t,r,p_{r}),\beta).$

Let there exist a Hamiltonian $H=H(t,r,p_{r})$ such that
$\partial_{t}Z=J\nabla H(t,Z)$ is satisfied for $Z=(R,P)$. In addition, we
assume that

  * (a)

$\partial_{e}\tilde{Q}\neq 0$ on the support of $Q$ and $U^{\prime}_{Q}(0)=0$;

  * (b)

$\partial_{r}R>0$ and there is a function $\sigma(t)$ such that

$\lim_{\delta\to 0}\frac{R(t,\delta,p_{r})}{\delta}=\sigma(t);$ (1.22)

  * (c)

for the Jacobian of the map $(r,p_{r})\mapsto(R(0,r,p_{r}),P(0,r,p_{r}))$ we
have

$\det\Big{(}\frac{\partial(R(0,r,p_{r}),P(0,r,p_{r}))}{\partial(r,p_{r})}\Big{)}=1.$
(1.23)

Then, defining $\alpha=U^{\prime\prime}_{Q}(0)$, we must have

$U^{\prime}_{Q}(r)=\alpha r,\quad R(t,r)=\frac{r}{\phi(t)},\quad
P(t,r,p_{r})=\phi(t)p_{r}-\dot{\phi}(t)r,$

where $\phi$ solves

$\ddot{\phi}(t)=\alpha\,\Big{(}-\frac{1}{\phi(t)^{2}}+\frac{1}{\phi(t)^{3}}\Big{)}.$

###### Remark 1.2

(a) The proof in Section 2 is a physics-style calculation, and we are not very
precise about, for instance, the regularity of $\tilde{Q}$. However, this is
not the main focus of the paper and missing details could be filled in easily.

(b) The constant on the right-hand side of (1.23) needs not be $1$, any other
number $\neq 0$ would also work.

(c) Concerning hypothesis (b), it will turn out in the proof that
$\partial_{r}R\neq 0$, see (2.10) below. Thus we are going to assume
$\partial_{r}R>0$ without loss of generality. Also $R=R(t,r)$ will be shown to
be independent of $p_{r}$. Hence (1.22) means that in fact

$\sigma(t)=\lim_{\delta\to 0}\frac{R(t,\delta)}{\delta}=\partial_{r}R(t,0)$

is required to exist. This can be guaranteed for instance if we suppose that
$H\in C^{2}$.

(d) When we started to look into the question if in general solutions of the
form (1.21) could be found, this author was convinced that $\beta$ should play
no role for the argument, in the sense that everything will be constant in
$\beta$. Theorem 1.1 indicates that actually the situation is much more
complicated, and that, in a vague sense, some ‘phase mixing’ would be needed
in order that (1.21) could provide a time-dependent solution. $\diamondsuit$

##  2 Proof of Theorem 1.1

Since the solution $f$ is spherically symmetric by (1.21), its potential
$U_{f}$ satisfies

$\displaystyle\partial_{r}U_{f}(t,\tilde{r})$ $\displaystyle=$
$\displaystyle\frac{4\pi}{\tilde{r}^{2}}\int_{0}^{\tilde{r}}r^{2}\rho_{f}(t,r)\,dr=\frac{4\pi}{\tilde{r}^{2}}\int_{0}^{\tilde{r}}dr\,r^{2}\int
dvf(t,x,v)$ (2.1) $\displaystyle=$
$\displaystyle\frac{4\pi^{2}}{\tilde{r}^{2}}\int_{0}^{\tilde{r}}dr\int\int\,dp_{r}\,d\beta\,\tilde{Q}\Big{(}e_{Q}(R(t),P(t),B(t)),B(t)\Big{)},$

where the arguments of $(R,P,B)$ are $(t,r,p_{r},\beta)$ and we have used that
$dv=\frac{2\pi}{r^{2}}\,dp_{r}\,d\ell\,\ell=\frac{\pi}{r^{2}}\,dp_{r}\,d\beta$.
By definition,

$f(t,r,p_{r},\beta)=\tilde{Q}\Big{(}\frac{1}{2}\,P(t)^{2}+U(R(t))+\frac{B(t)}{2R(t)^{2}},B(t)\Big{)}$

for $U=U_{Q}$ denoting the potential generated by the steady state $Q$. From
the spherically symmetric version of the Vlasov equation (see [1]) we hence
obtain, for all $(t,r,p_{r},\beta)$,

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\partial_{t}f(t,r,p_{r},\beta)+p_{r}\,\partial_{r}f(t,r,p_{r},\beta)+\Big{(}\frac{\beta}{r^{3}}-\partial_{r}U_{f}(t,r)\Big{)}\,\partial_{p_{r}}f(t,r,p_{r},\beta)$
(2.2) $\displaystyle=$
$\displaystyle(\partial_{e}\tilde{Q})\Big{[}P(\partial_{t}P)+U^{\prime}(R)(\partial_{t}R)+\frac{1}{2R^{2}}\,(\partial_{t}B)-\frac{B}{R^{3}}\,(\partial_{t}R)\Big{]}+(\partial_{\beta}\tilde{Q})(\partial_{t}B)$
$\displaystyle+\,p_{r}(\partial_{e}\tilde{Q})\Big{[}P(\partial_{r}P)+U^{\prime}(R)(\partial_{r}R)+\frac{1}{2R^{2}}\,(\partial_{r}B)-\frac{B}{R^{3}}\,(\partial_{r}R)\Big{]}+p_{r}(\partial_{\beta}\tilde{Q})(\partial_{r}B)$
$\displaystyle+\,\Big{(}\frac{\beta}{r^{3}}-\partial_{r}U_{f}(t,r)\Big{)}(\partial_{e}\tilde{Q})\Big{[}P(\partial_{p_{r}}P)+U^{\prime}(R)(\partial_{p_{r}}R)+\frac{1}{2R^{2}}\,(\partial_{p_{r}}B)-\frac{B}{R^{3}}\,(\partial_{p_{r}}R)\Big{]}$
$\displaystyle+\,\Big{(}\frac{\beta}{r^{3}}-\partial_{r}U_{f}(t,r)\Big{)}(\partial_{\beta}\tilde{Q})(\partial_{p_{r}}B)$
$\displaystyle=$
$\displaystyle(\partial_{e}\tilde{Q})P\,\Big{[}\partial_{t}P+p_{r}(\partial_{r}P)+\Big{(}\frac{\beta}{r^{3}}-\partial_{r}U_{f}(t,r)\Big{)}(\partial_{p_{r}}P)\Big{]}$
$\displaystyle+(\partial_{e}\tilde{Q})\Big{(}U^{\prime}(R)-\frac{B}{R^{3}}\Big{)}\,\Big{[}\partial_{t}R+p_{r}(\partial_{r}R)+\Big{(}\frac{\beta}{r^{3}}-\partial_{r}U_{f}(t,r)\Big{)}\,(\partial_{p_{r}}R)\Big{]}$
$\displaystyle+\Big{(}(\partial_{e}\tilde{Q})\frac{1}{2R^{2}}+\partial_{\beta}\tilde{Q}\Big{)}\,\Big{[}\partial_{t}B+p_{r}(\partial_{r}B)+\Big{(}\frac{\beta}{r^{3}}-\partial_{r}U_{f}(t,r)\Big{)}(\partial_{p_{r}}B)\Big{]}.$

Since $B=\beta$ is constant by hypothesis, the last line drops out. Also
$\partial_{e}\tilde{Q}\neq 0$ on the support of $Q$, whence (2.2) reduces to

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
P\,\Big{[}\partial_{t}P+p_{r}(\partial_{r}P)+\Big{(}\frac{\beta}{r^{3}}-\partial_{r}U_{f}(t,r)\Big{)}(\partial_{p_{r}}P)\Big{]}$
(2.3)
$\displaystyle+\,\Big{(}U^{\prime}(R)-\frac{\beta}{R^{3}}\Big{)}\,\Big{[}\partial_{t}R+p_{r}(\partial_{r}R)+\Big{(}\frac{\beta}{r^{3}}-\partial_{r}U_{f}(t,r)\Big{)}\,(\partial_{p_{r}}R)\Big{]}.$

As $R$ and $P$ are assumed to be independent of $\beta$, we can compare the
coefficients of the powers $\beta^{0},\beta^{1},\beta^{2}$ in $\beta$ to
deduce that

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
P\Big{[}\partial_{t}P+p_{r}(\partial_{r}P)-\partial_{r}U_{f}(t,r)(\partial_{p_{r}}P)\Big{]}$
(2.4)
$\displaystyle+\,U^{\prime}(R)\Big{[}\partial_{t}R+p_{r}(\partial_{r}R)-\partial_{r}U_{f}(t,r)(\partial_{p_{r}}R)\Big{]},$
$\displaystyle 0$ $\displaystyle=$ $\displaystyle
P\frac{1}{r^{3}}(\partial_{p_{r}}P)+U^{\prime}(R)\frac{1}{r^{3}}(\partial_{p_{r}}R)-\frac{1}{R^{3}}\,\Big{[}\partial_{t}R+p_{r}(\partial_{r}R)-\partial_{r}U_{f}(t,r)(\partial_{p_{r}}R)\Big{]},$
(2.5) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\frac{1}{r^{3}R^{3}}(\partial_{p_{r}}R).$

Thus $\partial_{p_{r}}R=0$ and (2.4), (2.5) simplify to

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
P\Big{[}\partial_{t}P+p_{r}(\partial_{r}P)-\partial_{r}U_{f}(t,r)(\partial_{p_{r}}P)\Big{]}+U^{\prime}(R)\Big{[}\partial_{t}R+p_{r}(\partial_{r}R)\Big{]},$
(2.6) $\displaystyle 0$ $\displaystyle=$ $\displaystyle
P\frac{1}{r^{3}}(\partial_{p_{r}}P)-\frac{1}{R^{3}}\,\Big{[}\partial_{t}R+p_{r}(\partial_{r}R)\Big{]},$
(2.7)

which is a PDE system for $(R(t,r),P(t,r,p_{r}))$. Coming back to (2.1), we
have

$\displaystyle\partial_{r}U_{f}(t,\tilde{r})$ $\displaystyle=$
$\displaystyle\frac{4\pi^{2}}{\tilde{r}^{2}}\int_{0}^{\tilde{r}}dr\int\int\,dp_{r}\,d\beta\,\tilde{Q}\Big{(}\frac{1}{2}\,P(t)^{2}+U(R(t))+\frac{\beta}{2R(t)^{2}},\beta\Big{)}$
(2.8) $\displaystyle=$
$\displaystyle\frac{4\pi^{2}}{\tilde{r}^{2}}\int_{0}^{\infty}d\beta\int_{0}^{\tilde{r}}dr\int\,dp_{r}\,\tilde{Q}\Big{(}\frac{1}{2}\,P(t)^{2}+U(R(t))+\frac{\beta}{2R(t)^{2}},\beta\Big{)}.$

Let $\phi_{t}$ denote the solution map that is associated to
$\partial_{t}z=J\nabla H(t,z)$, i.e., $z(t)=\phi_{t}(r,p_{r})$ solves the
equation and satisfies $z(0)=(r,p_{r})$. Thus with
$Z(t,r,p_{r})=(R(t,r),P(t,r,p_{r}))$ we get

$Z(t)=\phi_{t}(Z(0)),$ (2.9)

since $Z$ is assumed to be a solution. Note that we do not suppose that
$Z(0)=(r,p_{r})$, since this is also not satisfied for the Kurth solution:
from (1.17), (1.18) we have
$(R_{\varepsilon},P_{\varepsilon})(0)=(r,p_{r}-\varepsilon r)$ in this case.
Since the system is Hamiltonian, each map $\phi_{t}$ is a symplectomorphism
[8, Lemma 1.10], and hence in particular $\det D\phi_{t}=1$ holds for its
Jacobian determinant. By (1.23) from assumption (c) we also know that $\det
DZ(0)=1$. Therefore (2.9) shows that $\det DZ(t)=1$, and hence

$(\partial_{r}R)(t,r)(\partial_{p_{r}}P)(t,r,p_{r})-(\partial_{p_{r}}R)(t,r)(\partial_{r}P)(t,r,p_{r})=1,$

which in our case is

$(\partial_{r}R)(t,r)(\partial_{p_{r}}P)(t,r,p_{r})=1$ (2.10)

for all $(t,r,p_{r})$. At fixed $t$ we are going to apply the change of
variables

$(r,p_{r})\mapsto(R(t,r),P(t,r,p_{r}))=Z(t,r,p_{r})=(R,P)$

to (2.8), which has $\det DZ(t,r,p_{r})=1$ by (2.10) and $R(t,0)=0$ due to
assumption (b). Then we get

$\displaystyle\partial_{r}U_{f}(t,\tilde{r})$ $\displaystyle=$
$\displaystyle\frac{4\pi^{2}}{\tilde{r}^{2}}\int_{0}^{\infty}d\beta\int_{0}^{R(t,\tilde{r})}dR\int\,dP\,\tilde{Q}\Big{(}\frac{1}{2}\,P^{2}+U(R)+\frac{\beta}{2R^{2}},\beta\Big{)}$
(2.11) $\displaystyle=$
$\displaystyle\frac{4\pi^{2}}{\tilde{r}^{2}}\int_{0}^{R(t,\tilde{r})}dR\int\,dP\int\,d\beta\,\tilde{Q}\Big{(}\frac{1}{2}\,P^{2}+U(R)+\frac{\beta}{2R^{2}},\beta\Big{)}.$

On the other hand,

$\displaystyle U^{\prime}(\tilde{r})$ $\displaystyle=$
$\displaystyle\frac{4\pi}{\tilde{r}^{2}}\int_{0}^{\tilde{r}}r^{2}\rho_{Q}(r)\,dr=\frac{4\pi}{\tilde{r}^{2}}\int_{0}^{\tilde{r}}dr\,r^{2}\int
dv\,Q(x,v)$ (2.12) $\displaystyle=$
$\displaystyle\frac{4\pi^{2}}{\tilde{r}^{2}}\int_{0}^{\tilde{r}}dr\int
dp_{r}\int
d\beta\,\tilde{Q}\Big{(}\frac{1}{2}\,p_{r}^{2}+U(r)+\frac{\beta}{2r^{2}},\beta\Big{)}.$

Comparing (2.11) to (2.12), we have shown that

$\partial_{r}U_{f}(t,r)=\frac{R(t,r)^{2}}{r^{2}}\,U^{\prime}(R(t,r))$ (2.13)

is verified. From (2.6), (2.13) and (2.7) it follows that

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
P(\partial_{t}P+p_{r}(\partial_{r}P))-P\,\partial_{r}U_{f}(t,r)(\partial_{p_{r}}P)+U^{\prime}(R)(\partial_{t}R+p_{r}(\partial_{r}R))$
$\displaystyle=$ $\displaystyle
P(\partial_{t}P+p_{r}(\partial_{r}P))-P\,\frac{R^{2}}{r^{2}}\,U^{\prime}(R)\,(\partial_{p_{r}}P)+U^{\prime}(R)(\partial_{t}R+p_{r}(\partial_{r}R))$
$\displaystyle=$ $\displaystyle
P(\partial_{t}P+p_{r}(\partial_{r}P))-P\,\frac{R^{2}}{r^{2}}\,U^{\prime}(R)\,(\partial_{p_{r}}P)+R^{3}P\frac{1}{r^{3}}(\partial_{p_{r}}P)\,U^{\prime}(R)$
$\displaystyle=$
$\displaystyle\frac{P}{r^{3}}\Big{[}r^{3}(\partial_{t}P+p_{r}(\partial_{r}P))-rR^{2}\,U^{\prime}(R)\,(\partial_{p_{r}}P)+R^{3}(\partial_{p_{r}}P)\,U^{\prime}(R)\Big{]},$

so that

$r^{3}(\partial_{t}P+p_{r}(\partial_{r}P))+R^{2}(R-r)\,U^{\prime}(R)\,(\partial_{p_{r}}P)=0.$
(2.14)

Furthermore, by (2.7),

$\partial_{p_{r}}P^{2}=2P(\partial_{p_{r}}P)=2\,\frac{r^{3}}{R^{3}}\,(\partial_{t}R+p_{r}(\partial_{r}R)),$

and since $R$ is independent of $p_{r}$, integration $\int dp_{r}$ yields

$P^{2}(t,r,p_{r})-P^{2}(t,r,0)=2\,\frac{r^{3}}{R^{3}}\,p_{r}\,(\partial_{t}R)+\frac{r^{3}}{R^{3}}\,p_{r}^{2}\,(\partial_{r}R).$
(2.15)

In addition, integration $\int dp_{r}$ of (2.10) leads to

$(\partial_{r}R)(P(t,r,p_{r})-P(t,r,0))=p_{r}.$

Combining this relation with (2.15), we get

$P(t,r,p_{r})+P(t,r,0)=2\,(\partial_{r}R)\frac{r^{3}}{R^{3}}\,(\partial_{t}R)+(\partial_{r}R)\frac{r^{3}}{R^{3}}\,p_{r}\,(\partial_{r}R),$

which in turn implies that

$\displaystyle P(t,r,p_{r})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Big{(}P(t,r,p_{r})-P(t,r,0)+P(t,r,p_{r})+P(t,r,0)\Big{)}$
(2.16) $\displaystyle=$
$\displaystyle\frac{p_{r}}{2\,\partial_{r}R}+(\partial_{r}R)\,\frac{r^{3}}{R^{3}}\,\Big{(}\partial_{t}R+\frac{1}{2}\,p_{r}\,(\partial_{r}R)\Big{)};$

there are only $R$’s on the right-hand side, which are independent of $p_{r}$.
If we take the derivative w.r. to $t$, we obtain

$\displaystyle\partial_{t}P$ $\displaystyle=$
$\displaystyle-\frac{p_{r}}{2\,(\partial_{r}R)^{2}}\,(\partial^{2}_{tr}R)+(\partial^{2}_{tr}R)\,\frac{r^{3}}{R^{3}}\,\Big{(}\partial_{t}R+\frac{1}{2}\,p_{r}\,(\partial_{r}R)\Big{)}$
$\displaystyle-\,(\partial_{r}R)\,(\partial_{t}R)\,\frac{3r^{3}}{R^{4}}\,\Big{(}\partial_{t}R+\frac{1}{2}\,p_{r}\,(\partial_{r}R)\Big{)}+(\partial_{r}R)\,\frac{r^{3}}{R^{3}}\,\Big{(}\partial^{2}_{tt}R+\frac{1}{2}\,p_{r}\,(\partial^{2}_{tr}R)\Big{)}.$

Due to (2.14) and (2.10), we get

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
r^{3}\partial_{t}P+r^{3}p_{r}(\partial_{r}P)+R^{2}(R-r)\,U^{\prime}(R)\,(\partial_{p_{r}}P)$
$\displaystyle=$ $\displaystyle
r^{3}\Big{[}-\frac{p_{r}}{2\,(\partial_{r}R)^{2}}\,(\partial^{2}_{tr}R)+(\partial^{2}_{tr}R)\,\frac{r^{3}}{R^{3}}\,\Big{(}\partial_{t}R+\frac{1}{2}\,p_{r}\,(\partial_{r}R)\Big{)}$
$\displaystyle\hskip
15.00002pt-\,(\partial_{r}R)\,(\partial_{t}R)\,\frac{3r^{3}}{R^{4}}\,\Big{(}\partial_{t}R+\frac{1}{2}\,p_{r}\,(\partial_{r}R)\Big{)}+(\partial_{r}R)\,\frac{r^{3}}{R^{3}}\,\Big{(}\partial^{2}_{tt}R+\frac{1}{2}\,p_{r}\,(\partial^{2}_{tr}R)\Big{)}+p_{r}(\partial_{r}P)\Big{]}$
$\displaystyle+R^{2}(R-r)\,U^{\prime}(R)\,\frac{1}{(\partial_{r}R)}.$

Taking $p_{r}=0$, this yields

$0=\frac{r^{6}}{R^{3}}\Big{[}(\partial^{2}_{tr}R)\,(\partial_{t}R)-(\partial_{r}R)\,(\partial_{t}R)^{2}\,\frac{3}{R}+(\partial_{r}R)\,(\partial^{2}_{tt}R)\Big{]}+R^{2}(R-r)\,U^{\prime}(R)\,\frac{1}{(\partial_{r}R)}.$
(2.17)

If $t$ is fixed, then the map
$(r,p_{r})\mapsto(R(t,r),P(t,r,p_{r}))=Z(t,r,p_{r})$ is symplectic, owing to
(2.10), and it is generated by the ‘point transformation’ $r\mapsto R(t,r)$.
Thus, using [9, equ. (1.44)], there is a scalar function $v=v(t,r)$ such that

$P(t,r,p_{r})=\frac{1}{\partial_{r}R(t,r)}\,(p_{r}-\partial_{r}v(t,r)).$
(2.18)

Therefore, from (2.16),

$\frac{1}{\partial_{r}R}\,(p_{r}-\partial_{r}v(t,r))=P=\frac{p_{r}}{2\,\partial_{r}R}+(\partial_{r}R)\,\frac{r^{3}}{R^{3}}\,\Big{(}\partial_{t}R+\frac{1}{2}\,p_{r}\,(\partial_{r}R)\Big{)}.$
(2.19)

Taking once again $p_{r}=0$, we see that

$\partial_{r}v(t,r)=-(\partial_{r}R)^{2}\,(\partial_{t}R)\,\frac{r^{3}}{R^{3}}.$
(2.20)

If we plug this relation back to (2.19), it follows that

$(\partial_{r}R)^{3}\,\frac{r^{3}}{R^{3}}=1,$

or

$\partial_{r}R=\frac{R}{r}.$

For $\delta>0$ integration yields

$R(t,r)=r\,\frac{R(t,\delta)}{\delta}\,e^{C(t)}$

for a suitable function $C(t)$. By assumption, taking the limit $\delta\to 0$,
we get

$R(t,r)=a(t)r,$

where $a(t)=\sigma(t)e^{C(t)}$. Thus (2.17) simplifies to

$0=\frac{r^{3}}{a(t)^{3}}\Big{[}\dot{a}(t)^{2}r-a(t)\,\dot{a}(t)^{2}r^{2}\,\frac{3}{a(t)r}+a(t)\ddot{a}(t)r\Big{]}+a(t)^{2}r^{3}(a(t)-1)\,U^{\prime}(a(t)r)\,\frac{1}{a(t)},$

which is

$U^{\prime}(a(t)r)=\frac{2\dot{a}(t)^{2}-a(t)\ddot{a}(t)}{a(t)^{4}(a(t)-1)}\,r.$

Replacing $a(t)r$ by $r$, this leads to

$U^{\prime}(r)=\frac{2\dot{a}(t)^{2}-a(t)\ddot{a}(t)}{a(t)^{5}(a(t)-1)}\,r.$

Since the variables are separated, we deduce that there is
$\alpha\in\mathbb{R}$ such that

$\frac{U^{\prime}(r)}{r}=\alpha=\frac{2\dot{a}(t)^{2}-a(t)\ddot{a}(t)}{a(t)^{5}(a(t)-1)}$
(2.21)

for all $(t,r)$. Thus if we set $\phi(t)=\frac{1}{a(t)}$, then
$\dot{\phi}(t)=-\frac{\dot{a}(t)}{a(t)^{2}}$ and, by (2.21),

$\ddot{\phi}(t)=-\frac{\ddot{a}(t)}{a(t)^{2}}+2\,\frac{\dot{a}(t)^{2}}{a(t)^{3}}=\frac{1}{a(t)^{3}}\,(2\dot{a}(t)^{2}-a(t)\ddot{a}(t))=\alpha\,a(t)^{2}(a(t)-1)=\alpha\,\Big{(}-\frac{1}{\phi(t)^{2}}+\frac{1}{\phi(t)^{3}}\Big{)}.$

Moreover, we also have

$R(t,r)=\frac{r}{\phi(t)}.$

Next, using (2.20),

$\partial_{r}v(t,r)=-a(t)^{2}\,\dot{a}(t)r\,\frac{r^{3}}{a(t)^{3}r^{3}}=-\frac{\dot{a}(t)}{a(t)}\,r=\frac{\dot{\phi}(t)}{\phi(t)}\,r,$

and hence upon integration

$v(t,r)=\frac{1}{2}\,\frac{\dot{\phi}(t)}{\phi(t)}\,r^{2}+\gamma(t)$

for a suitable function $\gamma(t)$. Thus (2.18) implies that

$P(t,r,p_{r})=\frac{1}{\partial_{r}R(t,r)}\,(p_{r}-\partial_{r}v(t,r))=\phi(t)\,\Big{(}p_{r}-\frac{\dot{\phi}(t)}{\phi(t)}\,r\Big{)}=\phi(t)\,p_{r}-\dot{\phi}(t)\,r.$

Lastly, $U^{\prime}(0)=0$ in conjunction with (2.21) shows that

$\alpha=\lim_{r\to 0}\frac{U^{\prime}(r)}{r}=U^{\prime\prime}(0).$

This completes the proof of Theorem 1.1. $\Box$

##  3 Some technical results

###### Lemma 3.1

One has

$\int_{\mathbb{R}^{3}}Q_{{\rm K}}(x,v)\,dv=\frac{3}{4\pi}\,{\bf
1}_{B_{1}(0)}(x)$

Proof : We only consider $r=|x|<1$. First note that

$1-|x|^{2}-|v|^{2}+|x\wedge
v|^{2}=1-r^{2}-p_{r}^{2}-\frac{\ell^{2}}{r^{2}}+\ell^{2}>0$

means that $p_{r}^{2}<(1-r^{2})(1-\frac{\ell^{2}}{r^{2}})\leq 1-r^{2}$ and
also

$\ell^{2}<r^{2}-\frac{r^{2}p_{r}^{2}}{1-r^{2}}=:\ell_{0}^{2}.$

In spherical symmetry we have $dv=\frac{2\pi}{r^{2}}\,dp_{r}\,d\ell\,\ell$.
Therefore by (1),

$\displaystyle\int_{\mathbb{R}^{3}}Q_{{\rm K}}(x,v)\,dv$ $\displaystyle=$
$\displaystyle\frac{2\pi}{r^{2}}\int_{|p_{r}|\leq\sqrt{1-r^{2}}}dp_{r}\int_{0}^{\ell_{0}}d\ell\,\ell\,Q_{{\rm
K}}(r,p_{r},\ell)$ $\displaystyle=$
$\displaystyle\frac{3}{2\pi^{2}r^{2}}\int_{|p_{r}|\leq\sqrt{1-r^{2}}}dp_{r}\int_{0}^{\ell_{0}}d\ell\,\ell\,\frac{1}{(1-r^{2}-p_{r}^{2}-\frac{\ell^{2}}{r^{2}}+\ell^{2})^{1/2}}$
$\displaystyle=$
$\displaystyle\frac{3}{2\pi^{2}}\,\frac{1}{1-r^{2}}\int_{|p_{r}|\leq\sqrt{1-r^{2}}}dp_{r}\int_{0}^{\ell_{0}}d\ell\,(-1)\frac{d}{d\ell}\,\Big{(}1-r^{2}-p_{r}^{2}-\frac{\ell^{2}}{r^{2}}+\ell^{2}\Big{)}^{1/2}$
$\displaystyle=$
$\displaystyle\frac{3}{2\pi^{2}}\,\frac{1}{1-r^{2}}\int_{|p_{r}|\leq\sqrt{1-r^{2}}}dp_{r}\,(1-r^{2}-p_{r}^{2})^{1/2}$
$\displaystyle=$
$\displaystyle\frac{3}{\pi^{2}}\,\int_{0}^{1}ds\,(1-s^{2})^{1/2}$
$\displaystyle=$ $\displaystyle\frac{3}{4\pi},$

as was to be shown. $\Box$

###### Lemma 3.2

The function $f_{\varepsilon}$ from (1.6) is a (formal) solution to the
gravitational Vlasov-Poisson system.

Proof : To begin with, (1.6) yields

$\partial_{t}f_{\varepsilon}=-\frac{\dot{\phi}_{\varepsilon}}{\phi_{\varepsilon}^{2}}\,x\cdot((\nabla_{x}Q_{{\rm
K}})\circ\Lambda_{\varepsilon})-\ddot{\phi}_{\varepsilon}\,x\cdot((\nabla_{v}Q_{{\rm
K}})\circ\Lambda_{\varepsilon})+\dot{\phi}_{\varepsilon}\,v\cdot((\nabla_{v}Q_{{\rm
K}})\circ\Lambda_{\varepsilon}),$ (3.1)

where $\Lambda_{\varepsilon}$ is defined in (1.7). Observe that on $\\{Q_{{\rm
K}}\neq 0\\}$:

$\displaystyle\nabla_{x}Q_{{\rm K}}(x,v)$ $\displaystyle=$
$\displaystyle\Big{(}\frac{4\pi^{3}}{3}\Big{)}^{2}Q_{{\rm
K}}^{3}(x,v)\,((1-|v|^{2})x+\langle x,v\rangle v),$ (3.2)
$\displaystyle\nabla_{v}Q_{{\rm K}}(x,v)$ $\displaystyle=$
$\displaystyle\Big{(}\frac{4\pi^{3}}{3}\Big{)}^{2}Q_{{\rm
K}}^{3}(x,v)\,((1-|x|^{2})v+\langle x,v\rangle x).$ (3.3)

Furthermore,

$\displaystyle\nabla_{x}f_{\varepsilon}$ $\displaystyle=$
$\displaystyle\frac{1}{\phi_{\varepsilon}}\,((\nabla_{x}Q_{{\rm
K}})\circ\Lambda_{\varepsilon})-\dot{\phi}_{\varepsilon}\,((\nabla_{v}Q_{{\rm
K}})\circ\Lambda_{\varepsilon}),$ $\displaystyle\nabla_{v}f_{\varepsilon}$
$\displaystyle=$ $\displaystyle\phi_{\varepsilon}\,((\nabla_{v}Q_{{\rm
K}})\circ\Lambda_{\varepsilon}).$

On $\\{f_{\varepsilon}\neq 0\\}$ we have

$U_{\varepsilon}(t,x)=\frac{1}{\phi_{\varepsilon}(t)}\,\bigg{(}\frac{|x|^{2}}{2\phi_{\varepsilon}(t)^{2}}-\frac{3}{2}\bigg{)},\quad\nabla_{x}U_{\varepsilon}(t,x)=\frac{x}{\phi_{\varepsilon}(t)^{3}}.$

Thus we obtain from
$\ddot{\phi}_{\varepsilon}=-\frac{1}{\phi_{\varepsilon}^{2}}+\frac{1}{\phi_{\varepsilon}^{3}}$
that

$\displaystyle\partial_{t}f_{\varepsilon}+v\cdot\nabla_{x}f_{\varepsilon}-\nabla_{x}U_{\varepsilon}\cdot\nabla_{v}f_{\varepsilon}$
$\displaystyle=$
$\displaystyle-\frac{\dot{\phi}_{\varepsilon}}{\phi_{\varepsilon}^{2}}\,x\cdot((\nabla_{x}Q_{{\rm
K}})\circ\Lambda_{\varepsilon})-\ddot{\phi}_{\varepsilon}\,x\cdot((\nabla_{v}Q_{{\rm
K}})\circ\Lambda_{\varepsilon})+\dot{\phi}_{\varepsilon}\,v\cdot((\nabla_{v}Q_{{\rm
K}})\circ\Lambda_{\varepsilon})$
$\displaystyle+\,\frac{1}{\phi_{\varepsilon}}\,v\cdot((\nabla_{x}Q_{{\rm
K}})\circ\Lambda_{\varepsilon})-\dot{\phi}_{\varepsilon}\,v\cdot((\nabla_{v}Q_{{\rm
K}})\circ\Lambda_{\varepsilon})-\frac{x}{\phi_{\varepsilon}^{2}}\cdot((\nabla_{v}Q_{{\rm
K}})\circ\Lambda_{\varepsilon})$ $\displaystyle=$
$\displaystyle-\frac{\dot{\phi}_{\varepsilon}}{\phi_{\varepsilon}^{2}}\,x\cdot((\nabla_{x}Q_{{\rm
K}})\circ\Lambda_{\varepsilon})-\frac{1}{\phi_{\varepsilon}^{3}}\,x\cdot((\nabla_{v}Q_{{\rm
K}})\circ\Lambda_{\varepsilon})+\,\frac{1}{\phi_{\varepsilon}}\,v\cdot((\nabla_{x}Q_{{\rm
K}})\circ\Lambda_{\varepsilon})$ $\displaystyle=$
$\displaystyle\frac{1}{\phi_{\varepsilon}^{3}}\Big{(}\frac{4\pi^{3}}{3}\Big{)}^{2}(Q_{{\rm
K}}^{3}\circ\Lambda_{\varepsilon})\bigg{[}-\phi_{\varepsilon}\dot{\phi}_{\varepsilon}\,x\cdot\Big{(}(1-|\phi_{\varepsilon}v-\dot{\phi}_{\varepsilon}x|^{2})\frac{x}{\phi_{\varepsilon}}+\langle\frac{x}{\phi_{\varepsilon}},\phi_{\varepsilon}v-\dot{\phi}_{\varepsilon}x\rangle[\phi_{\varepsilon}v-\dot{\phi}_{\varepsilon}x]\Big{)}$
$\displaystyle\hskip
100.00015pt-\,x\cdot\Big{(}(1-\Big{|}\frac{x}{\phi_{\varepsilon}}\Big{|}^{2})[\phi_{\varepsilon}v-\dot{\phi}_{\varepsilon}x]+\langle\frac{x}{\phi_{\varepsilon}},\phi_{\varepsilon}v-\dot{\phi_{\varepsilon}}x\rangle\frac{x}{\phi_{\varepsilon}}\Big{)}$
$\displaystyle\hskip
100.00015pt+\,\phi_{\varepsilon}^{2}\,v\cdot\Big{(}(1-|\phi_{\varepsilon}v-\dot{\phi}_{\varepsilon}x|^{2})\frac{x}{\phi_{\varepsilon}}+\langle\frac{x}{\phi_{\varepsilon}},\phi_{\varepsilon}v-\dot{\phi}_{\varepsilon}x\rangle[\phi_{\varepsilon}v-\dot{\phi}_{\varepsilon}x]\Big{)}\bigg{]}$
$\displaystyle=$ $\displaystyle 0,$

where the last step requires some calculation. Apart from that,

$\Delta U_{\varepsilon}(t,x)=\frac{1}{\phi_{\varepsilon}(t)^{3}}\,\Delta
U_{Q_{{\rm
K}}}\Big{(}\frac{x}{\phi_{\varepsilon}(t)}\Big{)}=\frac{4\pi}{\phi_{\varepsilon}(t)^{3}}\,\rho_{Q_{{\rm
K}}}\Big{(}\frac{x}{\phi_{\varepsilon}(t)}\Big{)}=4\pi\rho_{\varepsilon}(t,x),$

which completes the somewhat formal argument. $\Box$

###### Lemma 3.3

The functions $R_{\varepsilon}$, $P_{\varepsilon}$ and $B_{\varepsilon}$ from
(1.17), (1.18) and (1.19), respectively, provide a solution to (1.20).

Proof : Since $H_{\varepsilon,{\rm K}}$ is independent of $\beta$, we drop
this variable. Then

$\partial_{r}H_{\varepsilon,{\rm
K}}(t,r,p_{r})=-\frac{\dot{\phi}_{\varepsilon}(t)}{\phi_{\varepsilon}(t)}\,p_{r}-(\dot{\phi}_{\varepsilon}^{2}(t)-\phi_{\varepsilon}(t)\ddot{\phi}_{\varepsilon}(t))\,r,\quad\partial_{p_{r}}H_{\varepsilon,{\rm
K}}(t,r,p_{r})=-\frac{\dot{\phi}_{\varepsilon}(t)}{\phi_{\varepsilon}(t)}\,r.$

This yields

$\dot{R}_{\varepsilon}(t)=-\frac{r}{\phi_{\varepsilon}^{2}(t)}\,\dot{\phi}_{\varepsilon}(t)=\partial_{p_{r}}H_{\varepsilon,{\rm
K}}(t,R_{\varepsilon}(t),P_{\varepsilon}(t))$

as well as

$\displaystyle\dot{P}_{\varepsilon}(t)$ $\displaystyle=$
$\displaystyle\dot{\phi}_{\varepsilon}(t)p_{r}-\ddot{\phi}_{\varepsilon}(t)r$
$\displaystyle=$
$\displaystyle\frac{\dot{\phi}_{\varepsilon}(t)}{\phi_{\varepsilon}(t)}\,\Big{(}\phi_{\varepsilon}(t)p_{r}-\dot{\phi}_{\varepsilon}(t)r\Big{)}+(\dot{\phi}_{\varepsilon}^{2}(t)-\phi_{\varepsilon}(t)\ddot{\phi}_{\varepsilon}(t))\,\frac{r}{\phi_{\varepsilon}(t)}$
$\displaystyle=$
$\displaystyle\frac{\dot{\phi}_{\varepsilon}(t)}{\phi_{\varepsilon}(t)}\,P_{\varepsilon}(t)+(\dot{\phi}_{\varepsilon}^{2}(t)-\phi_{\varepsilon}(t)\ddot{\phi}_{\varepsilon}(t))\,R_{\varepsilon}(t)$
$\displaystyle=$ $\displaystyle-\partial_{r}H_{\varepsilon,{\rm
K}}(t,R_{\varepsilon}(t),P_{\varepsilon}(t)),$

and altogether this yields (1.20). $\Box$

## References

  * [1] Batt J., Faltenbacher W. & Horst E.: Stationary spherically symmetric models in stellar dynamics, Arch. Rational Mech. Anal. 93, 159-183 (1986)
  * [2] Gidas B., Ni W.M. & Nirenberg L.: Symmetry and related properties via the maximum principle, Comm. Math. Phys. 68, 209-243 (1979)
  * [3] Glassey R.: The Cauchy Problem in Kinetic Theory, Society for Industrial and Applied Mathematics (SIAM), Philadelphia/PA 1996
  * [4] Glassey R. & Schaeffer J.: On symmetric solutions of the relativistic Vlasov-Poisson system, Comm. Math. Phys. 101, 459-473 (1985)
  * [5] Glassey R. & Schaeffer J.: Time decay for solutions to the linearized Vlasov equation, Transport Theory Statist. Phys. 23, 411-453 (1994)
  * [6] Glassey R. & Schaeffer J.: The “two and one-half-dimensional” relativistic Vlasov Maxwell system, Comm. Math. Phys. 185, 257-284 (1997)
  * [7] Kurth R.: A global particular solution to the initial value problem of stellar dynamics, Quart. Appl. Math. 36, 325-329 (1978/79)
  * [8] McDuff D. & Salamon D.: Introduction to Symplectic Topology, 2nd edition, Oxford University Press, Oxford 1998
  * [9] Moser J. & Zehnder E.: Notes on Dynamical Systems, American Mathematical Society, Providence/RI 2005
  * [10] Ramming T. & Rein G.: Oscillating solutions of the Vlasov-Poisson system–a numerical investigation, Phys. D 365, 72-79 (2018)
  * [11] Rein G.: Stability of spherically symmetric steady states in galactic dynamics against general perturbations, Arch. Ration. Mech. Anal. 161, 27-42 (2002)


# On the Relation between Prediction and Imputation Accuracy under Missing
Covariates

Burim Ramosaj∗, Justus Tulowietzki, Markus Pauly  Faculty of Statistics
Institute of Mathematical Statistics and Applications in Industry
Technical University of Dortmund
44227 Dortmund, Germany

###### Abstract

Missing covariates in regression or classification problems can prohibit the
direct use of advanced tools for further analysis. Recent research has
realized an increasing trend towards the usage of modern Machine Learning
algorithms for imputation. It originates from their capability of showing
favourable prediction accuracy in different learning problems. In this work,
we analyze through simulation the interaction between imputation accuracy and
prediction accuracy in regression learning problems with missing covariates
when Machine Learning based methods for both, imputation and prediction are
used. In addition, we explore imputation performance when using statistical
inference procedures in prediction settings, such as coverage rates of (valid)
prediction intervals. Our analysis is based on empirical datasets provided by
the UCI Machine Learning repository and an extensive simulation study.

###### keywords:

Missing Covariates , Imputation Accuracy , Prediction Accuracy , Prediction
Intervals , Bagging , Boosting

††journal: Journal of LaTeX Templates

##  1 Introduction

The presence of missing values in data preparation and data analysis makes the
usage of state-of-the art statistical methods difficult to apply. Seeking for
a universal answer to such problems was the main idea of [1], who introduced
multiple imputation. Through imputation, one provides data analysts
(sequences) of completed datasets, based on which, various data analysis
procedures can be conducted. An alternative to imputation is the usage of so-
called data adjustment methods: statistical methods, that directly treat
missing instances during training or parameter estimation such as the full-
information-maximum-likelihood method (see e.g. [2]) or the expectation-
maximization algorithm (c.f. [3]). A big disadvantage of these methods is the
expertise knowledge on theoretical model construction, where the likelihood
function of parameters of interest, for example, needs to be adopted
appropriately in order to account for missing information. Such examples can
be found in [4, 5, 6], e.g., where whole statistical testing procedures were
adjusted to account for missing values. It is already well-known that more
naive methods such as list-wise deletion or mean imputation can lead to severe
estimation bias, see e.g. [7, 8, 9, 10, 1]. Therefore, we do not discuss these
approaches further.

In the current paper we focus on regression problems, where we do not have
complete information on the set of covariates. Missing covariates in
supervised regression learning have been part in a variety of theoretical and
applicative research fields. In [10], for example, a theoretical analysis
based on maximum semiparametric likelihood for constructing consistent
regression estimates was conducted. While in [11, 12] or [13], for example,
multiple imputation is used as a tool in medical research for variable
selection or bias-reduction in parameter estimation. More recent research has
been focusing on Machine Learning (ML) based imputation. In [14, 15, 16, 17],
for example, the Random Forest method is used to impute missing values in
various datasets with mixed variable scales. Other ML-based methods such as
the $k$-nearest neighbor method, boosting machines or Bayesian Regression in
combination with classification and regression trees have been part of
(multiple) imputation, see e.g. [18, 19, 20, 21, 22].

Choosing an appropriate imputation method for missing data problems usually
depends on several aspects such as the capability to treat mixed-type data, to
impute reasonable values in variables of e.g. bounded support and to provide a
fast imputation solution. Imputation accuracy can usually be assessed through
the consideration of performance measures. Here, depending on the subsequent
application, one may focus on data reproducibility measured through normalized
root mean squared error ($NRMSE$) and proportion of false classification
($PFC$) or distribution preserving measures such as Kolmogorov-Smirnov based
statistics, see e.g. [14, 21, 23]. It is important to realize that these two
classes of performance measures for the evaluation of imputation do not always
agree, see, e.g., [24, 23]. In fact, one could be provided by an imputation
scheme with comparably low $NRMSE$ values, for which subsequent statistical
inference can have highly-inflated type-I error rates. Therefore, choosing
imputation methods by solely focusing on data reproducibility will not always
lead to correct inference procedures. To control for the latter, [1] coined
the term proper imputation: a property, that guarantees correct statistical
inference under the (multiple) imputation scheme.

While imputation accuracy can be accessed through data reproducibility or
distributional recovery, prediction performance of subsequently applied ML
methods (i.e. after imputation) are often evaluated using the mean-squared
error ($MSE$) or the misclassification error ($MCE$). Under missing
covariates, however, the sole focus on these measures is not sufficient, as
the disagreement between data reproducibility and statistical inference has
shown. In fact, beyond point-prediction we may also be interested in
uncertainty quantifiation in form of prediction intervals. The effect of
missing covariates on the latter remains mostly unknown. In this paper, we aim
to close this gap in empirical and simulation-based scenarios.

We thereby put a special focus on the relation between data reproducibility
and correct statistical inference for post-imputation prediction. While the
latter could also be measured by distributional discrepancies like in [23], we
aim to put a special focus on coverage rates and prediction interval lengths
in post-imputation prediction instead. The reason for this is the raise in ML-
based imputation methods and their competitive predictive performance for
both, imputation and prediction. Taking into account that data reproducibility
and statistical inference are not always in harmony under the missing
framework, an interesting question remains whether appropriate data
reproducible imputation schemes lead to favourable prediction results after
imputation. Furthermore, it is unknown whether imputation schemes with
comparably low $NRMSE$ values will lead to accurate predictive machines in
terms of delivering appropriate point predictions for future outcomes while
also correctly quantifying their uncertainty. Therefore, based on different
ML-methods in supervised regression learning problems, we (i) aim to clarify
the interaction between $NRMSE$ and $MSE$ as measures for data reproducibility
and predictive post-imputation ability. Furthermore, we (ii) aim to enlighten
the impact of accurate data-reproducible imputation schemes on predictive,
post-imputation statistical inference in terms of correct uncertainty
quantification. To measure the latter, we take into account correct coverage
rates and (narrow) interval lenghts of pointwise prediction intervals in post-
imputation settings obtained through ML-methods.

##  2 Measuring Accuracy

Measuring imputation accuracy can happen in many ways. A lot of research has
been focused on the general idea of reconstructing missing instances and being
as close as possible to the true underlying data (c.f. [14], [21]). Although
this approach seems reasonable, several disadvantages have been discovered
when using data reproducible measures such as squared error loss, especially
for later statistical inference; see e.g. [1, 24, 23]. Therefore, we discuss
several suitable measures for assessing prediction accuracy in regression
learning problems with missing covariates. In the sequel, we assume that we
have access to iid data collected in
$\mathcal{D}_{n}=\\{[\boldsymbol{X}_{i}^{\top},Y_{i}]^{\top}\in{\mathbb{R}}^{p+1}:i=1,\dots,n\\}$,
where

$\displaystyle Y_{i}=m(\boldsymbol{X}_{i})+\epsilon_{i}.$ (1)

Here, $m(\boldsymbol{x})=\mathbb{E}[Y_{1}|\boldsymbol{X}_{1}=\boldsymbol{x}]$
is the regression function, $\\{\epsilon_{i}\\}_{i=1}^{n}$ is a sequence of
iid random variables with $\mathbb{E}[\epsilon_{1}]=0$ and
$Var(\epsilon_{1})=\sigma^{2}\in(0,\infty)$ and we assume continuos covariates
$X_{i}$. Missing positions in the features
$\\{\boldsymbol{X}_{i}\\}_{i=1}^{n}$ are modeled by the indicator matrix
$\boldsymbol{R}=\\{R_{ij}\\}_{ij}\in\\{0,1\\}^{n\times p}$, where $R_{ij}=0$
indicates that the $i$-th observation of feature $j\in\\{1,\dots,p\\}$ is not
observable. Focusing on the general issue of predicting outcomes in regression
learning for new feature outcomes, we restrict our attention to data
reproducible accuracy measures and model prediction accuracy measures. In
order to cover statistical inference correctness for prediction, we use ML-
based prediction intervals as proposed in [25, 26, 27] to account for coverage
rates and interval lengths.

###  2.1 Imputation and Prediction Accuracy

In our setting, (missing) covariates are continuously scaled leading to the
usage of accuracy measures for continuous random variables. Regarding
imputation accuracy, we consider the $NRMSE$ formally given by

$\displaystyle NRMSE$ $\displaystyle=\frac{\sqrt{\sum\limits_{(i,j)\in
N_{mis}}(X_{ij}^{imp}-X_{ij}^{mis})^{2}}}{\sqrt{\sum\limits_{(i,j)\in
N_{mis}}(X_{ij}^{imp}-\bar{X}_{\cdot\cdot}^{mis})^{2}}},$ (2)

where $N_{mis}=\\{(i,j)\in\\{1,\dots,n\\}\times\\{1,\dots,p\\}|R_{ij}=0\\}$ is
the set of all observations and features with missing entries in those
positions. Here, $X_{ij}^{imp}$ denotes the imputed value of observation $i$
for variable $j$, while $X_{ij}^{mis}$ is the true, unobserved component of
those positions. $X_{\cdot\cdot}^{mis}$ is the mean of the sequence
$\\{X_{ij}:R_{ij}=0\\}$. Regarding the overall model performance on
prediction, we make use of the mean squared error

$\displaystyle MSE$
$\displaystyle=\mathbb{E}[(Y-\widehat{m}_{n}(\boldsymbol{X}))^{2}],$ (3)

where $\widehat{m}_{n}$ is an ML-based estimator of $m$ on $\mathcal{D}_{n}$
and $[\boldsymbol{X}^{\top},Y]^{\top}$ is an independent copy of
$[\boldsymbol{X}_{1}^{\top},Y_{1}]^{\top}$. Note that in the missing
framework, $m$ is estimated on the imputed dataset $\mathcal{D}_{n}^{imp}$,
while the $MSE$ is (usually) estimated using cross-validation procedures.

###  2.2 Prediction Intervals

Based on the methods for uncertainty quantification proposed in [25, 28, 26,
27], we make use of Random Forest based prediction intervals. In an extensive
simulation study in [27], it could be seen that other, ML-based prediction
intervals such as the (stochastic) gradient tree boosting (c.f. [29]) or the
XGBoost method (c.f. [30]) did not perform well under completely observed
covariates. Therefore, we restrict our attention to those already indicating
accurate coverage rates in completely observed settings. Meinshausen’s
Quantile Regression Forest (see [25]), for example, delivers a pointwise
prediction interval for an unseen feature point
$\boldsymbol{X}=\boldsymbol{x}$, which is formally given by

$\displaystyle PI_{QRF,n}$
$\displaystyle=[\widehat{Q}_{n,\alpha/2}(\boldsymbol{x});\quad\widehat{Q}_{n,1-\alpha/2}(\boldsymbol{x})],$
(4)

where
$\widehat{Q}_{n,\alpha/2}(\boldsymbol{x})=\inf\\{y|\widehat{F}_{n}(y|\boldsymbol{x})\geq\alpha/2\\}$
and $\widehat{F}_{n}(y|\boldsymbol{x})$ is a Random Forest based estimator for
the conditional distribution function $F(y|\boldsymbol{x})$ of
$Y|\boldsymbol{X}=\boldsymbol{x}$. Other prediction intervals based on the
Random Forest are, e.g., given in [26] and [27]. Following the same notation
as in [27], we refer with $\widehat{m}_{n,M}(\boldsymbol{x})$ to a Random
Forest prediction at $\boldsymbol{x}$, trained on $\mathcal{D}_{n}$ using $M$
decision trees, while $z_{1-\alpha}$ is the corresponding quantile of the
standard normal distribution. We consider the same residual variance
estimators as in [27], where $\widehat{\sigma}_{n,M}$ is the trivial residual
variance estimate, $\widehat{\sigma}_{n,Mcorrect}$ is the residual variance
estimate with finite-$M$ bias correction and $\widehat{\sigma}_{n,M;W}$ is a
weighted residual variance estimator, see also [28]. Moreover, we denote with
$\widehat{D}_{n,\alpha/2}^{RF}$ the empirical quantile of the Random Forest
Out-of-Bag residuals. With this notation we obtain four more prediction
intervals:

$\displaystyle PI_{n,empQ}(\boldsymbol{x})$
$\displaystyle=[\widehat{m}_{n,M}(\boldsymbol{x})+\widehat{D}_{n,\alpha/2}^{RF};\quad\widehat{m}_{n,M}(\boldsymbol{x})+\widehat{D}_{n,1-\alpha/2}^{RF}],$
(5) $\displaystyle PI_{n,ResVar}(\boldsymbol{x})$
$\displaystyle=[\widehat{m}_{n,M}(\boldsymbol{x})-z_{1-\alpha/2}\cdot\widehat{\sigma}_{n,M};\quad\widehat{m}_{n,M}(\boldsymbol{x})+z_{1-\alpha/2}\cdot\widehat{\sigma}_{n,M}],$
(6) $\displaystyle PI_{n,Mcorrect}(\boldsymbol{x})$
$\displaystyle=[\widehat{m}_{n,M}(\boldsymbol{x})-z_{1-\alpha/2}\cdot\widehat{\sigma}_{n,Mcorrect};\quad\widehat{m}_{n,M}(\boldsymbol{x})+z_{1-\alpha/2}\cdot\widehat{\sigma}_{n,Mcorrect}],$
(7) $\displaystyle PI_{n,weighted}(\boldsymbol{x})$
$\displaystyle=[\widehat{m}_{n,M}(\boldsymbol{x})-z_{1-\alpha/2}\cdot\widehat{\sigma}_{n,M;W};\quad\widehat{m}_{n,M}(\boldsymbol{x})+z_{1-\alpha/2}\cdot\widehat{\sigma}_{n,M;W}].$
(8)

For benchmarking we additionally consider a prediction interval that is build
under the linear model assumption. Imputation accuracy in inferential
prediction under missing covariates is then assessed by considering Monte-
Carlo estimated coverage rates and interval lengths.

##  3 Imputation and Prediction Models

We made use of the following state-of-the-art ML regression models for
prediction

  * •

the Random Forest as implemented in the R-package ranger (c.f [17]),

  * •

the (stochastic) gradient tree boosting (SGB) method from the R-package gbm
(c.f. [29]) and

  * •

the XGBoost method, also known as Queen of ML (c.f. [31]), as implemented in
the R-package xgboost.

For each of them, we fit a prediction model to the (imputed) data. Both
boosting methods rely on additive regression trees that are fitted
sequentially using the principles of gradient descent for loss minimization.
XGBoost, however, is slightly different by introducing extra randomization in
tree construction, a proportional shrinkage on the leaf nodes and a clever
penalization of trees. We refer to [29, 30, 32] for details on the concrete
algorithms. For benchmarking, a linear model is trained as well.

Although several imputation models are available on various (statistical)
software packages, we put a special focus on Random Forest based imputation
schemes and the multivariate imputation using chained equations (MICE)
procedure (c.f. [11], [14], [21]). The reasons for this are twofold, but both
have roots in the same theoretical issue called congeniality, see [33] for a
formal definition. First, using the same model for prediction and imputation
yields to inferential congeniality; a potential disagreement of both models,
however, can result in uncongenial (multiple) imputation methods. The latter
yields to invalid (multiple) imputation inference, as can be seen in [34] or
[33], for example. Secondly, focusing on Bayesian models for imputation such
as the MICE procedure, is in line with the general framework of congeniality
and the idea of (multiple) imputation. Although we do not directly compute
point-estimates during the analysis phase, interesting quantities in our
framework are Random Forest based prediction intervals and estimators of the
$MSE$.

missForest in R is an iterative algorithm developed by [14], that imputes
continuous and discrete random variables using trained Random Forests on
complete subsets of the data, and imputes missing values through prediction
with the trained Random Forest model. The process iterates in imputing missing
values until a pre-defined stopping criterion is met. Similar to the
missForest algorithm, we substituted the core learning method by other ML-
based methods such as the SGB method (in the sequel referred to as the gbm for
the algorithmic implementation) and the XGBoost method (in the sequel referred
to with xgboost for the algorithmic implementation). Both methods are
implemented in R using the same algorithmic framework as missForest, while
substituting the Random Forest method with the SGB resp. XGBoost. That means
that we train the SGB resp. XGBoost on (complete) subsets of the data and
impute missing values through the prediction of the trained model in an
iterative fashion.

MICE is a family of Bayesian imputation models developed in [35, 36]. Under
the normality assumption (i.e. MICE NORM), the method assumes a (Bayesian)
linear regression model, where every parameter in that model is drawn from
suitable priors. The predictive mean matching approach (MICE PMM) is similar
to MICE NORM, but does not impute missing values through the prediction of
those points using the Bayesian linear model, but randomly selects among
observed points that are closest to the same model prediction as MICE NORM. In
addition to these methods, MICE enables the implementation of Random Forest
based methods, referred to as MICE RF, see, e.g., [37]. The latter assumes a
modified Random Forest, where additional randomization is applied compared to
the missForest. For example, instead of simply predicting missing values
through averaging observations in leaf nodes, the method randomly selects
them. In addition, in the complete subset of the data determined for training
the Random Forest, potential missing values are not initially imputed by mean
or mode values, but by random draws among observed values. In the sequel, we
refer to the algorithmic implementation in R of all these method using the
terms mice_norm, mice_pmm and mice_rf.

##  4 Simulation Design

Our simulation design is separated in two parts. In the first part, empirical
data from the UCI Machine Learning Repository covering regression learning
problems are considered for the purpose of measuring imputation and prediction
accuracy. The following five datasets are considered:

  1. 1.

The Airfoile Data consists of $(p+1)=6$ variables measured in $n=1,503$
observations, where the target variable is the scaled sound pressure level
measured in decibel. The aim of this study conducted by NASA is to detect the
impact of physical shapes of airfoils on the produced noise.

  2. 2.

In the Concrete Data, $(p+1)=7$ variables are measured in $n=1,030$
observations. The target variable is the concrete compressive strength
measured in MPa units.

  3. 3.

The aim in the QSAR Data is to predict aquatic toxicity for a certain fish
species. It consists of $(p+1)=9$ variables measured in $n=546$ observations.

  4. 4.

The Real Estate Data has $(p+1)=6$ variables and $n=413$ observations. The aim
is to build a prediction model for house price developments in the area of New
Taipei City in Taiwan.

  5. 5.

The Power Plant Data consists of $n=9,568$ observations with $(p+1)=5$
variables. The actual dataset is much larger in terms of observations, but
only the first $9,568$ are selected to speed the computations. The aim of this
dataset is to predict the electric power generation of a water power-plant in
Turkey.

For each dataset, missing values under the MCAR scheme have been inserted with
$r\in\\{0.1,0.2,0.3,0.5\\}$ missing rates. Then, missing values have been
(once) imputed with the imputation methods mentioned in Section 3. The whole
process is iterated using $MC_{imp}=500$ Monte-Carlo iterates. Based on each
imputed dataset, all of the above mentioned prediction models are trained and
their prediction accuracy is measured using a five-fold cross-validated $MSE$.
Regarding hyper-parameter tuning of the various prediction models, we
conducted a grid-search using a ten-fold cross-validation procedure with ten
replications on the completely observed data, prior to the generation of
missing values. This has been conducted using the R-function trainControl of
the caret-package [38].

In the second part of our simulation study, synthetic data has been generated
with missing covariates to detect the effect of imputation accuracy on
prediction interval coverage rates. Here, we have focused on pointwise
prediction intervals. For sample sizes $n\in\\{100,500,1000\\}$, regression
learning problems of the form
$\\{[\boldsymbol{X}_{i}^{\top},Y_{i}]^{\top}\\}_{i=1}^{n}$ have been generated
using a $p=10$ dimensional covariate space and model $(\ref{RegModel})$, where
$X_{i}\overset{\text{iid}}{\sim}\mathcal{N}_{p=10}(0,\Sigma)$ and
$\epsilon_{i}\overset{\text{iid}}{\sim}\mathcal{N}(0,\sigma^{2})$ are
simulated independent of each other. Missing values have been inserted under
the MCAR scheme using various missing rates $r_{PI}\in\\{0.1,0.2,0.3\\}$.
Regarding the functional relationship between features and response, different
regression functions with coefficient
$\boldsymbol{\beta}_{0}=[2,4,2,-3,1,7,-4,0,0,0]^{\top}$ are used such as

  1. 1.

a linear model:
$m(\mathbf{x}_{i})=\mathbf{x}_{i}^{\top}\boldsymbol{\beta}_{0}$,

  2. 2.

a polynomial model:
$m(\mathbf{x}_{i})=\sum\limits_{j=1}^{p}\beta_{0,j}x_{i,j}^{j}$,

  3. 3.

a trigonometric model:
$m(\mathbf{x}_{i})=2\cdot\sin(\mathbf{x}_{i}^{\top}\boldsymbol{\beta}_{0}+2)$,

  4. 4.

and a non-continuous model:

$\displaystyle m(\mathbf{x}_{i})$
$\displaystyle=\begin{cases}\beta_{0,1}x_{i,1}+\beta_{0,2}x_{i,2}+\beta_{0,3}x_{i,3},&\text{
if }x_{i,3}>0.5,\\\ \beta_{0,4}x_{i,4}+\beta_{0,5}x_{i,5}+3&\text{ if
}x_{i,3}\leq 0.5.\end{cases}$

In order to capture potential dependencies among the features, various choices
for the covariance matrix $\Sigma$ are considered: a positive auto-regressive,
negative auto-regressive, compound symmetric, Toeplitz and the scaled identity
structure. In addition, we aim to take care of the systematic variation
originating from $m(\boldsymbol{X}_{1})$, and the noise $\epsilon_{1}$, by
choosing $\sigma^{2}$ in such a way that the signal-to-noise ratio
$SN:=Var(m(\boldsymbol{X}_{1}))/\sigma^{2}=1$. Finally, using $MC_{PI}=1,000$
Monte-Carlo iterations, prediction interval performance of the intervals
proposed in Section 2 are evaluated by approximating coverage rates and
(average) interval lengths over the Monte-Carlo iterates.

##  5 Simulation Results

We present the results for each parts of the simulation study separately,
starting with the empirical data sets with artificially inserted missing
values.

###  5.1 Results on Imputation Accuracy and Model Prediction Accuracy

In this section, we present the results for the empirical data analysis based
on the Airfoile dataset using the imputation and prediction accuracy measures
described in Section 2 for evaluation. We thereby focus on the Random Forest
and the XGBoost prediction model. The results of the linear and the SGB model
as well as the results for all other datasets are given in the supplement (see
Figures 1 - 19) and summarized at the end of this section.

Random Forest as Prediction Model. Figure 5.1 summarizes for each imputation
method the imputation error ($NRMSE$) and the model prediction error ($MSE$)
over $MC_{imp}=500$ Monte-Carlo iterates using the Random Forest method for
prediction on the imputed dataset. On average, the smallest imputation error
measured with the $NRMSE$ could be attained when using missForest and the gbm
imputation method. In addition, these methods yielded low variations in
$NRMSE$ across the Monte-Carlo iterates. In contrast, the mice_norm, mice_pmm
and mice_rf behaved similarly resulting into largest $NRMSE$ values across the
different imputation schemes with an increased variation in $NRMSE$ values.
The xgboost method performed slightly worse than missForest and gbm, when
focusing on imputation accuracy. In addition, all methods seem to be more or
less robust towards an increased missing rate. Interesting is the fact that
volatility decreases, as missing rates increase for the MICE procedures.
Prediction accuracy measured in terms of cross-validated $MSE$ using the
Random Forest model is lowest under the missForest, xgboost and gbm, which
corresponds with the $NRMSE$ results. As expected, the estimated $MSE$ suffers
from missing covariates and the effect gets worse with an increased missing
rate. Hence, model prediction accuracy suffers from an increased amount of
missing values, independent of the used imputation scheme. Although
congeniality was defined for valid statistical inference procedures, the
effect of using the same method for imputation and prediction seem to have
also a positive effect on model prediction accuracy.

Figure 1:  Imputation (left column) and prediction accuracy (right column)
measured by $NRMSE$ and $MSE$ using the Random Forest method for predicting
scaled sound pressure in the Airfoile dataset under various missing rates
(listed row-wise). The $MSE$ is estimated based on a five-fold cross-
validation procedure on the imputed dataset. The boxplots are over $500$
Monte-Carlo iterates. True refers to the Airfoile dataset without any missing
values and a fitted Random Forest on the complete dataset.

XGBoost as Prediction Model. Switching the prediction method to XGBoost, one
realizes an increase in model prediction accuracy for missing rates up to
$20\%$ as can be seen in Figure 5.1. In addition, for those missing rates, the
xgboost imputation is competitive to the missForest method but looses in
accuracy for larger missing rates compared to the missForest. Different to the
Random Forest, the XGBoost prediction method is more sensitive towards an
increased missing rate. Although under the completely observed framework the
XGBoost method performed best in terms of estimated $MSE$, the results
indicate that missing covariates can disturb the ranking. In fact, for missing
rates $r\geq 30\%$, the Random Forest exhibits a better prediction accuracy.

Figure 2:  Imputation (left column) and prediction accuracy (right column)
measured by $NRMSE$ and $MSE$ using the XGBoost method for predicting scaled
sound pressure in the Airfoile dataset under various missing rates (listed
row-wise). The $MSE$ is estimated based on a five-fold cross-validation
procedure on the imputed dataset. The boxplots are over $500$ Monte-Carlo
iterates. True refers to the Airfoile dataset without any missing values and a
fitted XGBoost on the complete dataset. Internal is the internal missing value
treatment of the XGBoost method without prior imputation.

Other Prediction models. Using the linear model as the prediction model
resulted in worse prediction accuracy with $MSE$ values ranging from $25$
($r=10\%$) to $45$ ($r=50\%$). For all missing scenarios, using the missForest
or the gbm method for imputation before prediction with the linear model
resulted in the lowest $MSE$. The results for the SGB method were even worse
with $MSE$ values between $80$ and $99$. As a surprising result prediction
accuracy measured in terms of cross-validated $MSE$ decreased with increasing
missing rate. A potential source of this effect could be the general weakness
of the SGB in the Airfoile dataset without any missing values. After inserting
and imputing missing values, which can yield to distributional changes of the
data, it seems that the SGB method benefits from these effects. However, model
prediction accuracy is still not satisfactory, see Figure 2 in the supplement.

Other Datasets. For the other datasets, similar effects could be obtained. The
Random Forest and the XGBoost showed the best prediction accuracy, see Figures
4 - 19 in the supplement. Again, larger missing rates affected model
prediction accuracy for the XGBoost method, but the Random Forest was more
robust to them overcoming XGBoost prediction performance measured in cross-
validated $MSE$ for larger missing rates. Overall, $NRMSE$ and cross-validated
$MSE$ seem to be positively associated to each other. Hence, more accurate
imputation models seem to yield better model prediction measured by $MSE$.

###  5.2 Results on Prediction Coverage and Length

Using the prediction intervals in Section 2, we present coverage rates and
interval lengths of pointwise prediction intervals in simulated data. Both
quantities have been computed using $1,000$ Monte-Carlo iterations with sample
sizes $n\in\\{100,500,1000\\}$. The boxplots presented here (see Figures 5.2,
5.2, 5.2 and 5.2) and in the supplement (see Figures 19 - 30) spread over the
different covariance structures used during the simulation. Every row
corresponds to one of the simulated missing rates $r\in\\{0.1,0.2,0.3\\}$,
while the columns reflect the different Random Forest based prediction
intervals. The left column summarizes the results for the Random Forest based
prediction interval using empirical quantiles ($PI_{n,empQ}$), the center
column reflects the Random Forest prediction interval using the simple
residual variance estimator on Out-of-Bag errors ($PI_{n,ResVar}$), while the
right column summarizes the Random Forest based prediction interval using the
weighted residual variance estimator ($PI_{n,weighted}$). We shifted the
results of $PI_{QRF,n}$, $PI_{n,MCorrect}$ and the prediction interval based
on the linear model to the supplement (see Figures 19, 20, 21, 22, 24, 26, 28
and 30 in the supplement). Under the complete case scenarios, the latter
methods did not show comparably well coverage rates as $PI_{n,emQ}$,
$PI_{n,ResVar}$ and $PI_{n,weighted}$. For imputed missing covariates, the
methods performed less accurate in terms of correct coverage rates, when
comparing them with $PI_{n,emQ}$, $PI_{n,ResVar}$ and $PI_{n,weighted}$.
Although the interval lengths of $PI_{QRF,n}$, $PI_{n,MCorrect}$ and the
linear model were, on average smaller, coverage rate was not sufficient to
make them competitive with $PI_{n,emQ}$, $PI_{n,ResVar}$ and
$PI_{n,weighted}$. For prediction intervals that underestimated the $0.95$
threshold in the complete case scenario, we could observe more accurate
coverage rates for larger missing rates. It seems that larger missing rates
increase coverage rates for the $PI_{QRF,n}$ and $PI_{n,MCorrect}$ methods,
independent of the used imputation scheme.

In Figure 5.2 the boxplots of the linear regression model are presented. In
general the usage of Random Forest based prediction intervals with empirical
quantiles ($PI_{n,empQ}$) or simple variance estimation ($PI_{n,ResVar}$) show
competitive behaviour in the complete case scenario. When considering the
various imputation schemes, under the different missing rates, it can be seen
that coverage rate slightly suffers compared to the complete case. To be more
precise, larger missing rates lead to slightly larger coverage rates for
$PI_{n,emQ}$, $PI_{n,ResVar}$ and $PI_{n,weighted}$. For the Random Forest
based prediction interval with weighted residual variance, this effect seems
to be positive, i.e. larger missing rates will lead to better coverage rates
for $PI_{n,weighted}$. Comparing the results with the previous findings, we
see that the xgboost yields, on average, the best coverage results across the
different imputation schemes. While the MICE procedures did not reveal
competitive performance in model prediction accuracy, the mice_norm method
under the linear model performed similar to the missForest procedure when
comparing coverage rates.

Figure 3:  Boxplots of prediction coverage rates under the linear model. The
variation is over the different covariance structures of the features. Each
row corresponds to one of the missing rates $r\in\\{0.1,0.2,0.3\\}$, while
each column to the following prediction intervals: $PI_{n,empQ}$,
$PI_{n,ResVar}$ and $PI_{n,weighted}$. The tripple (red, green and blue)
correspond to the sample sizes $n\in(100,500,1000)$.

Figure 5.2 summarizes coverage rates of pointwise prediction intervals under
the trigonometric model. Similar to the linear case, all three methods
$PI_{n,empQ}$, $PI_{n,ResVar}$ and $PI_{n,weighted}$ yield accurate coverage
rates showing better approximation to the $0.95$ threshold when the sample
size increase under the complete observation case. On average, the xgboost
imputation method remains competitive compared to the other imputation
methods. Slightly different to the linear case, the mice_norm approach gains
in correct coverage rate approximation compared to the missForest, together
with the mice_pmm approach. Nevertheless, the approximations between
mice_norm,mice_pmm and missForest are close to each other. As mentioned
earlier, $PI_{n,weighted}$ turns more accurate in terms of correct coverage
rates, when the missing rate increases.

Figure 4:  Boxplots of prediction coverage rates under the trigonometric
model. The variation is over the different covariance structures of the
features. Each row corresponds to one of the missing rates
$r\in\\{0.1,0.2,0.3\\}$, while each column to the following prediction
intervals: $PI_{n,empQ}$, $PI_{n,ResVar}$ and $PI_{n,weighted}$. The tripple
(red, green and blue) correspond to the sample sizes $n\in(100,500,1000)$.

Similar results compared to the linear and trigonometric case could be
obtained for the polynomial model and the non-continuous model. Boxplots of
the coverage rates can be found in Figures 24 and 28 of the supplement.

Regarding the length of the intervals for the linear model (see Figure 5.2),
the prediction interval $PI_{n,empQ}$ and the parametric interval
$PI_{n,ResVar}$ yield similar interval lengths. Under imputed missing
covariates, however, the $PI_{n,empQ}$ interval leads to slightly smaller
intervals than $PI_{n,ResVar}$. Nevertheless, the prediction interval based on
the weighted residual variance estimator $PI_{n,wighted}$ had the smallest
intervals on average. This comes with the cost of less accurate coverage rates
as can be seen in Figure 5.2. In addition, independent of the used prediction
interval, an increased missing rate yields to larger intervals making the
learning method such as the Random Forest more insecure about future
predictions. Regarding the used imputation method, almost all imputation
methods result into similar interval lengths. On average, the missForest
method has slightly smaller intervals comparable to the xgboost imputation.

Figure 5: Boxplots of prediction interval lengths under the linear model. The
variation is over the different covariance structures of the features. Each
row corresponds to one of the missing rates $r\in\\{0.1,0.2,0.3\\}$, while
each column to the following prediction intervals: $PI_{n,empQ}$,
$PI_{n,ResVar}$ and $PI_{n,weighted}$. The tripple (red, green and blue)
correspond to the sample sizes $n\in(100,500,1000)$.

Similar results on prediction lengths are obtained with other models.
Considering the trigonometric function as in Figure 5.2, it can be seen that
$PI_{n,empQ}$ results into slightly smaller intervals than $PI_{n,ResVar}$.
However, the interval length for the empirical quantiles under the
trigonometric model are more robust towards dependent covariates. Similar to
the linear case, $PI_{n,weighted}$ results into the smallest interval lengths,
but suffers from less accurate coverage. Furthermore, all imputation methods
behave similar with respect to prediction interval lengths under the
trigonometric case and other models (see Figures 21 and 22 in the supplement)
. It can be seen that Random Forest based prediction intervals are, more or
less, universally applicable to the different imputation schemes used in this
scenario yielding similar interval lengths.

Figure 6: Boxplots of prediction interval lengths under the trigonometric
model. The variation is over the different covariance structures of the
features. Each row corresponds to one of the missing rates
$r\in\\{0.1,0.2,0.3\\}$, while each column to the following prediction
intervals: $PI_{n,empQ}$, $PI_{n,ResVar}$ and $PI_{n,weighted}$. The tripple
(red, green and blue) correspond to the sample sizes $n\in(100,500,1000)$.

In summary, Random Forest based prediction intervals with imputed missing
covariates yield slightly wider intervals compared to the regression framework
without missing values. For prediction intervals that are underestimating the
true coverage rate, such as $PI_{QRF,n}$, $PI_{n,MCorrect}$ and
$PI_{n,weighted}$, an increased missing rate had positive effects on the
coverage rate. Overall, missForest and xgboost were competitive imputation
schemes when considering accurate coverage rates and interval lengths using
the $PI_{n,empQ}$ and $PI_{n,ResVar}$ intervals. mice_norm resulted into
similar, but slightly less accurate coverages compared to missForest and
xgboost

##  6 Conclusions

Missing covariates in regression learning problems are common issues in data
analysis problems. Providing a general approach for enabling the application
of various analysis models is often obtained through imputation. The use of
ML-based methods in this framework has obtained an increased attention during
the last decade, since fast and easy to use ML methods such as the Random
Forest may provide us with quick and accurate fixes of data analysis problems.
In our work, we put a special focus on variants of ML-based methods for
imputation, that mainly rely on decision trees as base learners and their
aggregation is conducted in a Random Forest based fashion or a boosting
approach. We aimed to shed light into the general issue when and which
imputation method should be used for missing covariates in regression learning
problems that provide accurate point predictions with correct uncertainty
ranges. To provide an answer to this, we conducted empirical analyses and
simulations which are led by the following questions: Does an imputation
scheme with a low imputation error (measured with the $NRMSE$) automatically
provide us with accurate model prediction performance (in terms of cross-
validated $MSE$)? How do ML-based imputation methods perform in estimating
uncertainty ranges for future prediction points in form of pointwise
prediction intervals? Are the results in harmony; that is, does an accurate
imputation method with a low $NRMSE$ provide us with good model prediction
accuracy measured in $MSE$ while delivering accurate and narrow prediction
interval lengths?

By analyzing empirical data from the UCI Machine Learning repository, we found
out that imputation methods with low imputation error measured with the
$NRMSE$ yield better model prediction measured by cross-validated $MSE$.
Especially for larger missing rates, the usage of the same ML method for both,
imputation and prediction was beneficial. This is in line with the
congeniality assumption in [33]; a theoretical term that (partly) guarantees
correct inference after (multiple) imputation. In particular the missForest
and our modified xgboost method for imputation yielded preferable results in
terms of low imputation error and good model prediction. It is expected that
ML methods with accurate model prediction capabilities measured in $MSE$ can
be transformed to be used as an imputation method yielding low imputation
errors, too. Regarding statistical inference procedures in prediction settings
such as the construction of valid prediction intervals, Random Forest based
imputation schemes such as the missForest and the xgboost method yielded
competitive coverage rates and interval lengths. In addition, the MICE
procedure with a Bayesian linear regression and normal assumption was under
the aspect of correct coverage rates and interval lengths, competitive too.
However, the method did not reveal low imputation error and overall good model
prediction. Hence, based on our findings, the missForest and the xgboost
method in combination with Random Forest based prediction intervals using
empirical quantiles resp. Out-of-Bag estimated residual variances are
competitive in three aspects: providing low imputation errors measured with
the $NRMSE$, yielding comparably low model prediction errors measured by
cross-validated $MSE$ and providing comparably accurate prediction interval
coverage rates and narrow widths using Random Forest based intervals.
Regarding the latter, our results also indicate that these intervals are
competitively applicable to a wide range of imputation schemes.
Future research will focus on a theoretical exploration of the interaction
between the $NRMSE$ and $MSE$ and the effect of the considered imputation
methods on uncertainty estimators. Insights into their theoretical interaction
will provide additional information to the general issue that imputation is
not just prediction.

## Acknowledgment

Burim Ramosaj’s work is funded by the Ministry of Culture and Science of the
state of NRW (MKW NRW) through the research grand programme KI-Starter. The
authors gratefully acknowledge the computing time provided on the Linux HPC
cluster at Technical University Dortmund (LiDO3), partially funded in the
course of the Large-Scale Equipment Initiative by the German Research
Foundation (DFG) as project 271512359.

## References

  * [1] D. B. Rubin, Multiple imputation for nonresponse in surveys, Vol. 81, John Wiley & Sons, (2004).
  * [2] C. K. Enders, The Performance of the Full Information Maximum Likelihood Estimator in Multiple Regression Models with Missing Data, Educational and Psychological Measurement 61 (5) (2001) 713–740.
  * [3] N. J. Horton, N. M. Laird, Maximum Likelihood Analysis of Generalized Linear models with Missing Covariates, Statistical Methods in Medical Research 8 (1) (1999) 37–50.
  * [4] L. Amro, M. Pauly, Permuting incomplete paired data: a novel exact and asymptotic correct randomization test, Journal of Statistical Computation and Simulation 87 (6) (2017) 1148–1159.
  * [5] L. Amro, F. Konietschke, M. Pauly, Multiplication-combination tests for incomplete paired data, Statistics in medicine 38 (17) (2019) 3243–3255.
  * [6] L. Amro, M. Pauly, B. Ramosaj, Asymptotic-based bootstrap approach for matched pairs with missingness in a single arm, Biometrical Journal 63 (7) (2021) 1389–1405.
  * [7] S. Greenland, W. D. Finkle, A Critical Look at Methods for Handling Missing Covariates in Epidemiologic Regression Analyses, American Journal of Epidemiology 142 (12) (1995) 1255–1264.
  * [8] J. W. Graham, S. M. Hofer, D. P. MacKinnon, Maximizing the Usefulness of Data Obtained with Planned Missing Value Patterns: An Application of Maximum Likelihood Procedures, Multivariate Behavioral Research 31 (2) (1996) 197–218.
  * [9] M. P. Jones, Indicator and Stratification Methods for Missing Explanatory Variables in Multiple Linear Regression, Journal of the American Statistical Association 91 (433) (1996) 222–230.
  * [10] H. Y. Chen, Nonparametric and Semiparametric Models for Missing Covariates in Parametric Regression, Journal of the American Statistical Association 99 (468) (2004) 1176–1189.
  * [11] S. Van Buuren, H. C. Boshuizen, D. L. Knook, Multiple Imputation of Missing Blood Pressure Covariates in Survival Analysis, Statistics in Medicine 18 (6) (1999) 681–694.
  * [12] X. Yang, T. R. Belin, W. J. Boscardin, Imputation and Variable Selection in Linear Regression Models with Missing Covariates, Biometrics 61 (2) (2005) 498–506.
  * [13] J. A. Sterne, I. R. White, J. B. Carlin, M. Spratt, P. Royston, M. G. Kenward, A. M. Wood, J. R. Carpenter, Multiple imputation for missing data in epidemiological and clinical research: potential and pitfalls, BMJ 338 (2009).
  * [14] D. J. Stekhoven, P. Bühlmann, MissForest—non-parametric missing value imputation for mixed-type data, Bioinformatics 28 (1) (2012) 112–118.
  * [15] A. D. Shah, J. W. Bartlett, J. Carpenter, O. Nicholas, H. Hemingway, Comparison of Random Forest and Parametric Imputation Models for Imputing Missing Data using MICE: A CALIBER Study, American Journal of Epidemiology 179 (6) (2014) 764–774.
  * [16] F. Tang, H. Ishwaran, Random forest missing data algorithms, Statistical Analysis and Data Mining: The ASA Data Science Journal 10 (6) (2017) 363–377.
  * [17] M. Mayer, M. M. Mayer, Package ‘missranger’ (2018).
  * [18] J. Chen, J. Shao, Nearest Neighbor Imputation for Survey Data, Journal of Official Statistics 16 (2) (2000) 113.
  * [19] D. Xu, M. J. Daniels, A. G. Winterstein, Sequential BART for imputation of missing covariates, Biostatistics 17 (3) (2016) 589–602.
  * [20] D. Dobler, S. Friedrich, M. Pauly, Nonparametric manova in mann-whitney effects, arXiv preprint arXiv:1712.06983 (2017).
  * [21] B. Ramosaj, M. Pauly, Predicting missing values: a comparative study on non-parametric approaches for imputation, Computational Statistics 34 (4) (2019) 1741–1764.
  * [22] X. Zhang, C. Yan, C. Gao, B. Malin, Y. Chen, XGBoost Imputation for Time Series Data, in: 2019 IEEE International Conference on Healthcare Informatics (ICHI), IEEE, (2019), pp. 1–3.
  * [23] M. Thurow, F. Dumpert, B. Ramosaj, M. Pauly, Goodness (of fit) of Imputation Accuracy: The GoodImpact Analysis, arXiv preprint arXiv:2101.07532 (2021).
  * [24] B. Ramosaj, L. Amro, M. Pauly, A cautionary tale on using imputation methods for inference in matched-pairs design, Bioinformatics 36 (10) (2020) 3099–3106.
  * [25] N. Meinshausen, Quantile Regression Forests, Journal of Machine Learning Research 7 (6) (2006).
  * [26] H. Zhang, J. Zimmerman, D. Nettleton, D. J. Nordman, Random Forest Prediction Intervals, The American Statistician (2019).
  * [27] B. Ramosaj, Interpretable Machines: Constructing Valid Prediction Intervals with Random Forests, arXiv preprint arXiv:2103.05766 (2021).
  * [28] B. Ramosaj, M. Pauly, Consistent estimation of residual variance with random forest out-of-bag errors, Statistics & Probability Letters 151 (2019) 49–57.
  * [29] J. H. Friedman, Stochastic gradient boosting, Computational Statistics & Data analysis 38 (4) (2002) 367–378.
  * [30] T. Chen, C. Guestrin, Xgboost: A scalable tree boosting system, in: Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining, 2016, pp. 785–794.
  * [31] V. Morde, XgBoost Algorithm: Long May She Reign!, https://towardsdatasc1ence.com/https-medium-com-vishalmorde-XgBoost-algorithmlong-shemay-rein-edd9f99be63d, [Online; accessed 03-December-2021] (2019).
  * [32] J. H. Friedman, The elements of statistical learning: Data mining, inference, and prediction, springer open, 2017.
  * [33] X.-L. Meng, Multiple-imputation Inferences with Uncongenial Sources of Input, Statistical Science (1994) 538–558.
  * [34] R. E. Fay, When Are Inferences from Multiple Imputation Valid?, US Census Bureau, (1992).
  * [35] S. Van Buuren, K. Groothuis-Oudshoorn, mice: Multivariate Imputation by Chained Equations in R, Journal of Statistical Software 45 (2011) 1–67.
  * [36] S. Van Buuren, Flexible Imputation of Missing Data, CRC press, (2018).
  * [37] L. L. Doove, S. Van Buuren, E. Dusseldorp, Recursive partitioning for missing data imputation in the presence of interaction effects, Computational Statistics & Data Analysis 72 (2014) 92–104.
  * [38] M. Kuhn, A Short Introduction to the caret Package, R Foundation for Statistical Computing 1 (2015) 1–10.