# The Infinitesimal Jackknife and Combinations of Models

Indrayudh Ghosal111Denotes equal lead contributions. $\;{}^{1}$, Yunzhe
Zhou11footnotemark: 1 $\;{}^{2}$, and Giles Hooker$\;{}^{3}$
1 Department of Statistics and Data Science, Cornell University
2 Department of Biostatistics, UC Berkeley
3 Department of Statistics, UC Berkeley

###### Abstract

The Infinitesimal Jackknife is a highly general method for estimating
variances, both in parametric models, and more recently for ensemble methods.
In this paper we extend the Infinitesimal Jackknife to estimate the covariance
between any two models. This can be used to provide uncertainty estimates for
combinations of models, or to construct test statistics for comparing
different models or combinations of models estimated from the same dataset. We
give the specific examples including boosted combination of models, the
comparisons of random forests and M-estimators, directional derivatives of
neural networks and ensemble XGBoost. We illustrate the efficacy of the test
through extensive simulations and Beijing Housing Data, and provide
theoretical guarantee for the consistency of the Infinitesimal Jackknife
covariance estimate.

The Infinitesimal Jackknife is a general method for estimating variances of
parametric models, and more recently also for some ensemble methods. In this
paper we extend the Infinitesimal Jackknife to estimate the covariance between
any two models. This can be used to quantify uncertainty for combinations of
models, or to construct test statistics for comparing different models or
ensembles of models estimated using the same training dataset. Specific
examples in this paper use boosted combinations of models like random forests
and M-estimators. We also use neural networks and XGBoost. We illustrate the
efficacy of the test through extensive simulations and real data examples
(Beijing Housing), and show theoretical consistency of the Infinitesimal
Jackknife covariance estimate.

Key Words: Infinitesimal Jackknife; M-estimator; Ensemble; Kernel;
V-statistics; Random Forests.

##  1 Introduction

This paper focuses on methods that enable uncertainty quantification for
combinations of models that are not immediately comparable. The jackknife is a
resampling technique that can be used for bias and variance estimation
[tukey1958bias]. Here the estimator is The estimator can be constructed by
repeatedly calculating the statistic, each time leaving one observation from
the sample out and averaging all estimates. As opposed to the ordinary
jackknife, the Infinitesimal Jackknife (IJ) looks at the behavior of the
statistic after giving an infinitesimally more weight to one observation (and
down-weighting the rest to keep the total weight constant) a statistic after
down-weighting each observation by an infinitesimal amount
[jaeckel1972infinitesimal]. In recent work, the IJ has been widely used to
quantify uncertainty for ensemble methods, most particularly for variance
estimation for the predictions of random forests [wager2014confidence,
athey2019generalized]. In addition, recent work has examined estimating the
covariance between two random forests - specifically when they are used as
base learners in a boosted setting [ghosal2020boosting]. However, the
definition of the IJ as proposed in efron1982jackknife is much more general
and can also be applied to a wide variety of parametric models.

This paper focuses on the application of IJ methods to combinations of models.
Within machine learning, combinations of models arise from mixture models and
from boosting methods [as explored in ghosal2020boosting]. They may also arise
if we wish to compare model predictions and provide uncertainty estimates for
their difference. This may be used for carried out for the purposes of
assessing goodness of fit, or when assessing the impact of a model or hyper-
parameter choices.

In this paper we demonstrate the use of the IJ for defining a general
covariance between the predictions of any two models for which IJ estimates
are available individually. Our demonstration makes unique and useful
contributions in several ways (do we need this sentence?) (jeff: we can remove
it.). First, we use the predictions of models as a means of incorporating non-
parametric functions into inference and employ IJ covariance definitions to
compare different models or combinations of models estimated from the same
training dataset. Secondly, apart from random forests, more general types of
models can be incorporated into our framework, including Generalized Linear
Models (GLMs) and kernel methods. Note that these models are incorporated in
very different ways. IJ can be inherently applied to random forests or general
M-estimators (in particular, GLMs and small neural networks). However, these
methods can also be employed with generic models – we use XGboost in
experiments – when they are combined with subsampling to produce an ensemble.
(I didn’t understand what the last three sentences are trying to say. Starting
from ”Secondly, …”) These variance estimates can then be used to create
intervals to describe the variability of the underlying model, although their
inferential validity will depend on the bias of that model.

The rest of the paper is organized as follows. Because we will be concerned
with the sample properties of confidence intervals and tests, we first
introduce a common simulation framework and terminology for all our
experiments in section 2. Then we give the definition of the IJ in section 3,
with the specific examples of M-estimators, random forests, and kernel
methods. In section 4, we show how to use the IJ to calculate the confidence
interval for the general boosting methods and local model modifications
(mostly used for bias corrections). This is a natural extensions of
ghosal2020boosting and ghosal2021generalised. Then we demonstrate how it can
be used to compare models in section 5 with generalised linear models as a
special case of M-estimators. In section 6, we also explore how to derive the
IJ for small neural network by regarding it as a M-estimator. This can make it
possible to compare neural network with other models. We next discuss the
general V-statistics framework to calculate the IJ for ensemble models and use
an example of the ensemble XGBoost for illustrations in section 7. Finally, we
apply these methods to Beijing House Data in section 8. All technical theorems
and proofs are provided in the Supplementary Appendix.

##  2 Simulation Framework

A key component of our contribution is to examine the many ways in which
models may be combined. In order to assess the finite sample performance of
our methods, we will use the following common simulation framework for all our
experiments. One or more of the definitions from this section will be used in
all the subsequent sections. In our studies, we consider only the regression
problems with response $Y|X=x\sim N(\eta(x),1)$ and $X\sim U[-1,1]^{6}$ with
the following generative models:

  * •

Friedman: $\eta(x)=10\sin(\pi
x_{1}x_{2})+20(x_{3}-\frac{1}{2})^{2}+10x_{4}+5x_{5}$.

  * •

Linear: $\eta(x)=x_{1}+x_{2}+x_{3}+x_{4}$

  * •

Constant: $\eta(x)=2$

but we note that the theory we produce readily extends to an exponential
family response. For each simulation, we fix a set of 5 (for model
comparisons) or 100 (for coverage calculations) query points at which our
estimated function will be evaluated. We then generate 200 data sets, each of
1000 points uniformly distributed on $[-1,1]^{6}$. This allows us to obtain a
distribution of predictions at each of the query points, along with the
covariance of predictions between them. The following metrics will be used
throughout:

  * •

Coverage of Expectation CoE (CoE):Let $\hat{f}_{ij}$ denote the prediction at
the $i$th query point and $j$th replication. $\hat{V}_{ij}$ is its variance
estimate. Then we construct 95% confidence intervals by $\hat{f}_{ij}\pm
1.96\times\sqrt{\hat{V}_{ij}}$. For $j=1,2,...,200$ in each fixed query point
$i$, we calculate a coverage probability by checking whether the expected
prediction value (approximated by $\frac{1}{200}\sum_{j=1}^{200}\hat{f}_{ij}$)
falls into this interval. We generate the boxplots for the coverage
probability of all the query points. We use this to assess the accuracy of our
estimates of variance and to distinguish this from model bias which also
enters the Coverage of Truth below.

  * •

Coverage of Truth (CoT): Borrowing the notations above, we instead calculate
the coverage probability by checking whether the true expectation $E(Y|X)$
value falls into the constructed 95% confidence intervals.

  * •

Coverage of Reproduction (CoR): In section 4 we use and evaluate intervals
that cover the value of an alternative $\hat{f}_{ij}$ obtained from
independent data set across 95% of replications.

We also summarize all the models that is considered in the simulations:

  * •

Random Forests: implemented by using sklearn package in python [scikit-learn].
We vary the number of maximum depth in the range of $\\{3,5,7,9\\}$ and also
consider a random forest with trees grown to full depth. We set the number of
trees to be 1000 or 5000 and use a subsample size of 200, taken with
replacement.

  * •

Generalized Linear Models(GLM): We consider the linear regression as an
example in this paper.

  * •

Neural Networks: implemented by using Tensorflow2.0 package in python
[abadi2016tensorflow]. We vary the number of hidden units in the range of
$\\{1,3,5,10,20\\}$ and consider ReLU and Sigmoid activation functions.

  * •

XGBoost: We use the xgboost package for implementation and use the default
hyperparameters () [chen2016xgboost].

##  3 Infinitesimal Jackknife

The IJ is a general-purpose framework for estimating the variance of any
statistic. efron1982jackknife defines the IJ for an estimate of the form
$\hat{\theta}(P^{0})$, where $P^{0}$ is the uniform probability distribution
over the empirical dataset. We re-write our estimate more generally as
$\hat{\theta}(P^{*})$ where $P^{*}$ is any re-weighting of the empirical
distribution. This is then approximated by the hyperplane tangent to the
surface $\hat{\theta}(P^{*})$ at the point $P^{*}=P^{0}$, i.e.
$\hat{\theta}(P^{*})\approx\hat{\theta}_{\text{TAN}}(P^{*})=\hat{\theta}(P^{0})+(P^{*}-P^{0})^{\top}U$,
where $U$ is a vector of the directional derivatives given by

$U_{i}=\lim_{\epsilon\to
0}\frac{\hat{\theta}(P^{0}+\epsilon(\delta_{i}-P^{0}))-\hat{\theta}(P^{0})}{\epsilon},i=1,\ldots,n$

where $\delta_{i}$ the $i$th coordinate vector. Under a suitable asymptotic
normal distribution for $P^{*}-P^{0}$ we can obtain the variance of
$\hat{\theta}_{\text{TAN}}(P^{*})$ to be

$\text{Var}[\hat{\theta}_{\text{TAN}}(P^{*})]=\frac{1}{n^{2}}\sum_{i=1}^{n}U_{i}^{2}$

This is the IJ variance estimate for the estimator $\hat{\theta}(P^{0})$.
Below we derive the form of IJ estimators for three specific models.

###  3.1 IJ for M-estimators

We first consider the general M-estimators as an example. M-estimators are a
broad class of extremum estimators for which the objective function is a
sample average [hayashi2000extremum]. Consider a data space
$\mathscr{Z}\subseteq\mathds{R}$ and the training data
$Z_{1},Z_{2},\cdots,Z_{n}$ as i.i.d copies of $Z\in\mathscr{Z}$. Then we
define the M-estimator by

$\hat{\theta}=\operatorname*{arg\,max}_{\theta\in\Theta}\mathds{E}_{\hat{\mathscr{D}}}[m(\theta,Z)]$

where $\mathds{E}_{\hat{\mathscr{D}}}$ is the expectation with regard to the
empirical distribution $\hat{\mathscr{D}}$ over the training data, that is
$\mathds{E}_{\hat{\mathscr{D}}}[m(\theta,Z)]=\frac{1}{n}\sum_{i=1}^{n}m(\theta,Z_{i})$.
$\Theta\subseteq\mathds{R}^{p}$ is the parametric space and $m(\theta,Z)$ is a
well-behaved function from $\Theta\times\mathscr{Z}$ to $\mathds{R}$. Since we
are interested in the predictions of M-estimators, we denote by
$\eta(\theta,x)$ the prediction function corresponding to any parameter
$\theta\in\Theta$ and query point $x\in\mathscr{Z}^{d}$. Then the directional
derivatives of the model predictions are given by

$\displaystyle
U_{i}(x)=-\nabla_{\theta}\eta(\hat{\theta},x)^{\top}\left[\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta},Z)]\right]^{-1}\nabla_{\theta}m(\hat{\theta},Z_{i})$
(3.1)

Thus, the variance of the model predictions can be given by

$\displaystyle\text{Var}(\eta(\hat{\theta},x))=\frac{1}{n^{2}}\sum_{i=1}^{n}U^{2}_{i}(x)$
$\displaystyle=\frac{1}{n^{2}}\sum_{i=1}^{n}\nabla_{\theta}\eta(\hat{\theta},x)^{\top}\left[\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta},Z)]\right]^{-1}\nabla_{\theta}m(\hat{\theta},Z_{i})\nabla_{\theta}m(\hat{\theta},Z_{i})^{\top}\left[\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta},Z)]\right]^{-1}\nabla_{\theta}\eta(\hat{\theta},x)^{\top}$
$\displaystyle=\nabla_{\theta}\eta(\hat{\theta},x)^{\top}\left[\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta},Z)]\right]^{-1}\nabla_{\theta}\eta(\hat{\theta},x)^{\top}$
$\displaystyle\approx\nabla_{\theta}\eta(\theta^{*},x)^{\top}\left[\mathds{E}_{\mathscr{D}}[\nabla_{\theta}^{2}m(\theta^{*},Z)]\right]^{-1}\nabla_{\theta}\eta(\theta^{*},x)^{\top}$

where $\mathds{E}_{\mathscr{D}}$ is the expectation with regard to the
distribution of $Z$ and
$\theta^{*}=\operatorname*{arg\,max}_{\theta\in\Theta}\mathds{E}[m(\theta,Z)]$
(the approximation holds true for large samples due to the law of large
numbers, the consistency of M-estimators and the continuity of $\eta$). This
results in the usual plug-in values for the sandwich form of the asymptotic
variance of $\eta(\hat{\theta},x)$.

The derivation of these directional derivatives and consistency of
corresponding variance estimate is shown in appendix B and appendix C. In
section 6, we will further discuss how to formalize small neural network model
as an M-estimator and calculate its directional derivatives.

###  3.2 IJ for Random Forests

We next consider the example of random forests. In efron2014estimation,
directional derivatives for a random forest are given by

$U_{i}(x)=n\cdot cov_{b}(N_{i,b},T_{b}(x))$

where the covariance is over $b=1,\dots,B$ trees, $N_{i,b}$ is the number of
times the $i$th training data point is included in the sample used to train
the $b$th tree and $T_{b}$ is the $b$th tree kernel.

(The derivation of the directional derivatives could probably be moved to the
appendix. It might even be moved to the XGBoost section. I don’t think
repeating the derivation here serves any purpose.) The calculation of these
directional derivatives is general for any ensemble model. Specifically, given
the observations $z_{1},z_{2},\cdots,z_{n}$, we define the reweighted
empirical probability distribution as

$\hat{F}^{*}:\text{mass}\;P^{*}_{i}\;\text{on}\;z_{i},\quad i=1,2,\cdots,n$

where the resampling vector
$P^{*}\in\Big{\\{}\Big{(}P^{*}_{1},P^{*}_{2},\cdots,P^{*}_{n}\Big{)}\Big{|}P^{*}_{i}\geq
0,\sum_{i=1}^{n}P^{*}_{i}=1\Big{\\}}$. In particular,
$P^{0}=(\frac{1}{n},\frac{1}{n},\cdots,\frac{1}{n})$. We define the
multinomial expectation of $\theta(P^{*})$ as

$\displaystyle S(P)=\mathbb{E}_{*}{[\hat{\theta}(P^{*})]},\quad
P^{*}\sim\frac{\text{Mult}_{n}(n,P)}{n},$

for any resampling vector $P$ and $\frac{\text{Mult}_{n}(n,P)}{n}$ is a
rescaled multinomial distribution, that is $n$ independent draws on $n$
categories each having probability $1/n$. $\hat{\theta}(P^{*})$ is the
prediction of the random forest and $\mathbb{E}_{*}$ indicates expectation
under the multinomial distribution of $P^{*}$. So $S(P^{0})$ is the ideal
smoothed bootstrap estimate. Define the directional derivative

$\displaystyle U_{i}=\lim_{\begin{subarray}{c}\epsilon\to
0\end{subarray}}\frac{S(P^{0}+\epsilon(\delta_{i}-P^{0}))-S(P^{0})}{\epsilon},\quad
i=1,\cdots,n$

where $\delta_{i}$ the $i$th coordinate vector. Formula (3.20) of
efron2014estimation gives

$\displaystyle
S(P^{0}+\epsilon(\delta_{i}-P^{0}))=S(P^{0})+n\cdot(cov_{*}(N_{i,*},T_{*}(x)))^{2}+o(\epsilon)$

where $cov_{*}$ represents the ideal covariance under the bootstrap sampling.
This yields

$U_{i}(x)=n\cdot cov_{*}(N_{i,*},T_{*}(x))$

In applications, we generate finite number of bootstrap replications and use

$U_{i}(x)=n\cdot cov_{b}(N_{i,b},T_{b}(x))$

as the estimate instead. In section 7, we will give an example of puting
XGBoost into an ensemble model and using the above framework to calculate its
directional derivatives.

Derivation of directional derivatives ends here

efron2014estimation proved that the variance estimator using these directional
derivatives is consistent for random forests where the training data for each
tree is a full bootstrap resample of the original training data, whereas
wager2017estimation did the same where each tree uses a random subsample
(without replacement) of the original training data.

Following similar arguments as in Appendix B.1 of ghosal2020boosting the
theoretical variance for a random forest prediction at a query point $x$ can
be given by
$\frac{k_{n}^{2}}{n}\zeta_{1,k_{n}}+\frac{1}{B_{n}}\zeta_{k_{n},k_{n}}$, where
$\frac{k_{n}^{2}}{n}\zeta_{1,k_{n}}$ is estimated with the usual IJ
directional derivatives, and $\zeta_{k_{n},k_{n}}=var_{*}(T_{*}(x))$, i.e.,
the theoretical variance between the tree kernels $T_{*}$, can be estimated by
$var_{b}(T_{b}(x))$.

####  3.2.1 Bias Corrections for IJ Estimates with Ensembles

Estimating $\frac{k_{n}^{2}}{n}\zeta_{1,k_{n}}$ with the infinitesimal
Jackknife method introduces an upward bias (zhou2021v, wager2017estimation,
wager2014confidence). The source of this bias is due to the fact that random
forests are incomplete U-statistics, i.e, the number of trees used ($B_{n}$)
are very small compared to the total number of possible trees that could be
used (the theoretical $U_{i}$ is derived for a model in which every subsample
is used. Hence we obtain a Monte Carlo estimate of it. assumes a model using
every possible subsample whereas the in practise we use a Monte Carlo
estimate). In practice $B_{n}\ll\binom{n}{k_{n}}$ for subsampling without
replacement and $B_{n}\ll(k_{n})^{n}$ for subsampling with replacement;
requiring a correction term as discussed in below.

  * •

Ranger’s Correction: Using the IJ directional derivates defined above, the
uncorrected variance estimate for a random forest is given by
$\sum_{i=1}^{n}(cov_{b}(N_{i,b},T_{b}(x)))^{2}$. But this is actually an
estimate for the population quantity
$\sum_{i=1}^{n}(cov_{*}(N_{i,*},T_{*}(x)))^{2}$, where the covariance is over
an infinite number of all possible trees. Focusing on the $i$th term we are
estimating $(cov_{*}(N_{i,*},T_{*}(x)))^{2}$ with
$(cov_{b}(N_{i,b},T_{b}(x)))^{2}$. If we define $N^{(c)}$ and $T^{(c)}$ to be
the mean-centered versions of the inbag count and trees respectively then
$\left(\mathds{E}[N^{(c)}_{i,*}T^{(c)}_{*}(x)]\right)^{2}$ is estimated by
$\left(\frac{1}{B_{n}}\displaystyle\sum_{b=1}^{B_{n}}N^{(c)}_{i,b}T^{(c)}_{b}(x)\right)^{2}$.
When using a finite ensemble, we write

$\frac{1}{B_{n}}\displaystyle\sum_{b=1}^{B_{n}}N^{(c)}_{i,b}T^{(c)}_{b}(x)=\mathds{E}[N^{(c)}_{i,*}T^{(c)}_{*}(x)]+e,$

where $e$ is the Monte Carlo error. We assume that this error $e$ is
independent of the population quantity
$\mathds{E}[N^{(c)}_{i,*}T^{(c)}_{*}(x)]$, and thus
$\left(\mathds{E}[N^{(c)}_{i,*}T^{(c)}_{*}(x)]+e\right)^{2}$ has an expected
value of
$\left(\mathds{E}[N^{(c)}_{i,*}T^{(c)}_{*}(x)]\right)^{2}+\mathds{E}[e^{2}]$,
where $\mathds{E}(e^{2})$ can be estimated by

$\hat{\mathds{E}}[e^{2}]=\frac{var_{b}(N^{(c)}_{i,b}T^{(c)}_{b}(x))}{B_{n}}\approx\frac{var_{b}(N^{(c)}_{i,b})\cdot
var_{b}(T^{(c)}_{b}(x))}{B_{n}}$

in which $N^{(c)}_{i,b}$ and $T^{(c)}_{b}(x)$ are assumed to be approximately
independent. Finally note that $var_{b}(N^{(c)}_{i,b})=var_{b}(N_{i,b})$ and
$var_{b}(T^{(c)}_{b}(x))=var_{b}(T_{b}(x))$. Hence we can estimate
$\sum_{i=1}^{n}(cov_{*}{N_{i,*},T_{*}(x)})^{2}$ by

$\sum_{i=1}^{n}(cov_{b}(N_{i,b},T_{b}(x)))^{2}-\frac{1}{B_{n}}\sum_{i=1}^{n}var_{b}(N_{i,b})\cdot
var_{b}(T_{b}(x)).$

The final corrected variance estimate for a random forest prediction at a
query point $x$ is thus given by

$\sum_{i=1}^{n}(cov_{b}(N_{i,b},T_{b}(x)))^{2}-\frac{1}{B_{n}}\left(\sum_{i=1}^{n}var_{b}(N_{i,b})-1\right)\cdot
var_{b}(T_{b}(x))$

We can extend this correction to the covariance terms as well. For two query
points $x_{1}$ and $x_{2}$ the original uncorrected covariance is given by

$\sum_{i=1}^{n}cov_{b}(N_{i,b},T_{b}(x_{1}))\cdot
cov_{b}(N_{i,b},T_{b}(x_{2}))$

and which is corrected by subtracting

$\frac{1}{B_{n}}\left(\sum_{i=1}^{n}var_{b}(N_{i,b})-1\right)\cdot
cov_{b}(T_{b}(x_{1}),T_{b}(x_{2})).$

We label this ”ranger’s correction” since it is the correction implemented in
the ranger package implementing random forests [rangerRpackage]. While
computationally simple, it relies on the strong underlying assumption that
$N^{(c)}_{i,b}$ and $T^{(c)}_{b}(x)$ are independent, which introduces bias
for the corrected variance estimator. The constructed confidence interval
based on this estimator could suffer from the under-coverage.

  * •

V-statistics based Correction: Another method for variance estimation
correction is through a ”V-statistics based correction” method. When the
ensemble is obtained using subsamples taken with replacement, the random
forests estimator can be regarded as a V-statistic [zhou2021v]. Similar to
sun2011efficient, an ANOVA-like estimator of variance can then be constructed.
Let $N_{i}=\sum_{b=1}^{B_{n}}N_{i,b}$,
$m_{i}(x)=\sum_{b=1}^{B_{n}}\frac{N_{i,b}}{N_{i}}T_{b}(x)$ and
$\bar{m}(x)=\frac{1}{n}\sum_{i=1}^{n}m_{i}(x)$. Define

$\displaystyle\text{SS}_{\tau}(x)=\sum_{i=1}^{n}N_{i}(m_{i}(x)-\bar{m}(x))^{2}\quad\text{and}\quad\text{SS}_{\epsilon}(x)=\sum_{i=1}^{n}\sum_{b=1}^{B_{n}}N_{i,b}(T_{b}(x)-m_{i}(x))^{2}$

A bias-corrected estimate can be given by

$\frac{\text{SS}_{\tau}(x)-(n-1){\hat{\sigma}}^{2}_{\epsilon}(x)}{C-\sum_{i=1}^{n}N^{2}_{i}/C},$

where $C=\sum_{i=1}^{n}N_{i}$ and
${\hat{\sigma}}^{2}_{\epsilon}(x)=\frac{\text{SS}_{\epsilon}(x)}{C-n}$.

Similarly, this can be also extended to the case of covariance estimation as
well. For two query points $x_{1}$ and $x_{2}$, define

$\displaystyle\text{SS}_{\tau}(x_{1},x_{2})$
$\displaystyle=\sum_{i=1}^{n}N_{i}(m_{i}(x_{1})-\bar{m}(x_{1}))(m_{i}(x_{2})-\bar{m}(x_{2}))$
$\displaystyle\text{SS}_{\epsilon}(x_{1},x_{2})$
$\displaystyle=\sum_{i=1}^{n}\sum_{b=1}^{B_{n}}N_{i,b}(T_{b}(x_{1})-m_{i}(x_{1}))(T_{b}(x_{2})-m_{i}(x_{2}))$

The bias-corrected estimate can be constructed by

$\frac{\text{SS}_{\tau}(x_{1},x_{2})-(n-1){\hat{\sigma}}^{2}_{\epsilon}(x_{1},x_{2})}{C-\sum_{i=1}^{n}N^{2}_{i}/C},$

where
${\hat{\sigma}}^{2}_{\epsilon}(x_{1},x_{2})=\frac{\text{SS}_{\epsilon}(x_{1},x_{2})}{C-n}$
for the natural extensions of $\text{SS}_{\epsilon}$ and $\text{SS}_{\tau}$ to
multiple query points.

This estimate is an unbiased (over the choice of subsamples) for the IJ
estimate that results from using an infinite number of all possible trees
[searle2009variance]. This framework can be also naturally extended to general
ensemble models rather than just random forests and provides a unified
treatment of variance estimation. In section 4, we will also incorporate the
ensemble XGBoost models into this framework. We point out that this method
does not work for U-statistics for the reason that its sampling schema is not
equivalent to sampling from the empirical distribution. There is also work
that develops an unbiased variance estimator based on incomplete U-statistics
for variance estimation of random forests, where the tree size is allowed to
be comparable with the sample size [xu2022variance]. However, we only focus on
the case when the subsampling is taken with replacement in this paper.

In order to illustrate the advantage of V-statistics based correction over
ranger’s correction, we present a small simulation to compare their CoE for
predictions by using the experiment settings mentioned in section 2. We
consider three data generating processes. We use the random forests with trees
grown to full depth and set the number of trees B to be 1000 or 5000. The
result is presented in Figure 1. We draw the boxplots for the CoE over 100
query points. From the plots, we observe that ranger’s correction suffers from
undercoverage, especially when the number of trees is not very large. In
contrast, V-statistics based correction can produce coverage probability close
to or above 0.95 with reasonable number of trees. In the rest of the sections
of this paper, we will use V-statistics based correction for estimating the
covariance matrix of random forests.

Figure 1: Comparison of the CoE (95% level) between ranger’s correction (red)
and V-stat based correction (blue). We consider three data generating
processes. The number of trees B is set to be 1000 or 5000. We draw the
boxplots for the CoE over 100 query points.

###  3.3 IJ for Kernel Methods

Kernel regression is a non-parametric technique to estimate the conditional
expectation of a random variable $Y$ relative to a variable $X$
[nadaraya1964estimating; watson1964smooth]. Given the dataset
$\\{X_{k},Y_{k}\\}_{k=1}^{n}$, the Nadaraya–Watson estimator for any query
point $x$ is defined as

$\displaystyle\hat{m}_{h}(x)=\frac{\sum_{k=1}^{n}K_{n}(x-X_{k})Y_{k}}{\sum_{k=1}^{n}K_{n}(x-X_{k})}$

We define
$\hat{m}_{h}(x)=\frac{\sum_{k=1}^{n}C_{k}Y_{k}}{\sum_{k=1}^{n}C_{k}}$, where
$C_{k}=K_{n}(x-X_{k})$. Then the directional derivatives of $\hat{m}_{h}(x)$
are given by

$U_{i}=\lim_{\epsilon\to
0}\frac{\hat{\theta}(P^{0}+\epsilon(\delta_{i}-P^{0}))-\hat{\theta}(P^{0})}{\epsilon},\text{
where
}\hat{\theta}(P)=\frac{\sum_{k=1}^{n}P_{k}C_{k}Y_{k}}{\sum_{k=1}^{n}P_{k}C_{k}},\sum_{k=1}^{n}P_{k}=1$

We let $p=\sum_{k=1}^{n}C_{k}Y_{k}$, $q=\sum_{k=1}^{n}C_{k}$, $r=nC_{i}Y_{i}$
and $s=nC_{i}$ then we get

$\displaystyle U_{i}$ $\displaystyle=\lim_{\epsilon\to
0}\frac{1}{\epsilon}\left[\frac{\frac{1}{n}(1-\epsilon+n\epsilon)C_{i}Y_{i}+\sum_{k\neq
i}\frac{1}{n}(1-\epsilon)C_{k}Y_{k}}{\frac{1}{n}(1-\epsilon+n\epsilon)C_{i}+\sum_{k\neq
i}\frac{1}{n}(1-\epsilon)C_{k}}-\frac{\sum_{k=1}^{n}C_{k}Y_{k}}{\sum_{k=1}^{n}C_{k}}\right]$
$\displaystyle=\lim_{\epsilon\to 0}\frac{1}{\epsilon}\left[\frac{n\epsilon
C_{i}Y_{i}+(1-\epsilon)\sum_{k=1}^{n}C_{k}Y_{k}}{n\epsilon
C_{i}+(1-\epsilon)\sum_{k=1}^{n}C_{k}}-\frac{\sum_{k=1}^{n}C_{k}Y_{k}}{\sum_{k=1}^{n}C_{k}}\right]$
$\displaystyle=\lim_{\epsilon\to
0}\frac{1}{\epsilon}\left[\frac{(1-\epsilon)p+\epsilon
r}{(1-\epsilon)q+\epsilon s}-\frac{p}{q}\right]$
$\displaystyle=\lim_{\epsilon\to 0}\frac{(1-\epsilon)pq+\epsilon
rq-(1-\epsilon)pq-\epsilon ps}{\epsilon q((1-\epsilon)q+\epsilon s)}$
$\displaystyle=\frac{rq-ps}{q^{2}}$

With the directional derivatives above, we can calculate variance estimates
and produce confidence intervals for kernel regression predictions. In section
4.3, we will show that the local model modifications proposed in
[lu2021unified] share similar structures with kernel regression so we can
follow the calculations above to obtain their directional derivatives.

##  4 Boosting Methods

If a model is inadequate to capture the underlying relationships between $y$
and $x$ in a data set, another model can be added to pick up the rest of the
signal. We define this process to be be boosting in reference to gradient
boosting [friedman2001greedy] and since it is sequential and the second model
is used to boost the signal. Previous work in ghosal2020boosting has shown
that boosting a random forest with another one achieves significant reduction
in bias and test-set MSE with a small increase in variance. An extension of
this can be applied to non-Gaussian signals where we initialise with an MLE
type estimator and then boost with (one or two) random forests with
appropriately defined “pseudo-residuals” as the training signals
[ghosal2021generalised]. We further extend the idea in this paper to boost any
model with another one and obtaining uncertainty quantification via the IJ
covariance estimate.

A common problem shared by this boosting technique is the inheritance of
variability from one stage of modelling to the next. This does not impede the
modelling process on empirical data but asymptotic analysis of these boosted
models are easier if we assume the following condition.

######  Condition 1.

We’re given an initial dataset $D=(Z_{i})_{i=1}^{n}=(Y_{i},X_{i})_{i=1}^{n}$.
Suppose $g$ is a fixed function of $X$ (the predictor) and $\hat{g}$ is an
unbiased estimator of $g$ based on this dataset. Assume a model $\hat{f}$ is
fit with training signal $\big{(}h(\hat{g}(X_{i}),Z_{i})\big{)}_{i=1}^{n}$. If
$\check{f}$ is the model fitted the same way as $\hat{f}$ but with training
signals $\big{(}h(g(X_{i}),Z_{i})\big{)}_{i=1}^{n}$ then for any query point
$x$

$\frac{\hat{f}(x)-\check{f}(x)}{\sqrt{var(\hat{f}(x))}}\xrightarrow{p}0.$

The function $h$ may include training weights if needed. Special cases of this
condition has been used in previous literature:

  * •

In ghosal2020boosting $\hat{g}$ and $\hat{f}$ are the first and second (base
and boost) stage random forests respectively.

  * •

In ghosal2021generalised $\hat{g}$ could be the MLE-type estimator
($\hat{\eta}_{MLE}^{(0)}$) or the sum of $\hat{\eta}_{MLE}^{(0)}$ and the
first random forest ($\hat{f}_{1}$) - in those cases $h$ would be the training
signals and weights for the random forests $\hat{f}_{1}$ and $\hat{f}_{2}$
respectively.

For two-step boosted models our estimators will be of the form
$\hat{f}_{1}=\hat{F}_{1}$ and $\hat{f}_{2}=\hat{F}_{1}+\hat{F}_{2}$. Suppose
in this case the directional derivatives for $\hat{F}_{1}$ and $\hat{F}_{2}$
are given by $U_{i}^{(1)}(x)$ and $U_{i}^{(2)}(x)$ respectively for any query
point $x$ and $i=1,\dots,n$ and also define

$\hat{\Sigma}^{(pq)}_{ij}=\frac{1}{n^{2}}\sum_{k=1}^{n}U_{k}^{(p)}(x_{i})U_{k}^{(q)}(x_{j}),\;p,q=1,2,\;i,j=1,\dots,m$

to be the covariance estimate for
$\Sigma^{(pq)}_{ij}:=cov(\hat{F}_{p}(x_{i}),\hat{F}_{q}(x_{j}))$. Then we can
construct valid confidence interval for $\hat{f}_{2}$. Specifically, the
covariance matrix of $\hat{f}_{2}$ over the query points can be calculated by
$\hat{\Sigma}=\hat{\Sigma}^{(11)}+\hat{\Sigma}^{(22)}+\hat{\Sigma}^{(12)}+\hat{\Sigma}^{(21)}$
which can then be used to construct confidence intervals.

###  4.1 Boosted Forests

ghosal2020boosting discussed a version of boosting models for Gaussian
responses. There the first stage model (base learner) was a random forest and
the subsequent boosted stage model was another random forest whose training
signals were the residuals of the first forest. Here, boosted forests provided
lower test set MSE and higher confidence interval coverage. With the
techniques of this paper, we compare the CoT in prediction between random
forests and boosted forest models.

We use the same simulation setup as in section 2. Figure 2 presents boxplots
of the CoT for 95% confidence intervals. We use a random forest with the
maximum depth of 3 as the initial estimator at the first step. In the second
step, we use random forests with different maximum depth (3,5,7,9, ”Full” on
the x axis) for boosting, where ”Full” represents the random forests with
trees grown to full depth. We also include the case when there is no second
step of boosting (red boxplot in each subplot). For comparison purpose, we
also use green star symbols to present the average CoE with 95% level over 100
query points.

Figure 2: CoT (95% level) for boosting of random forests. Three scenarios of
model signals are considered and the number of trees $B$ is set to be 1000 or
5000. We use a random forest with the maximum depth of 3 as the initial
estimator at the first step. In each subplot, the x axis represents the model
for the second boosting step. The red boxplot with x axis label ”Without
Boosting” is for the case when there is no second step of boosting. The number
”3,5,7,9” denotes the different values of maximum depth when using random
forests for boosting. And ”Full” represents the random forests with trees
grown to full depth for boosting. For comparison purpose, we also use the
green star symbols to present the average CoE with 95% level over 100 query
points.

A first observation is that under Linear and Friedman signals, the initial
random forests model suffers from low coverage of the truth and the boosting
process can efficiently improve it. This tells us that the initial model has
serious bias and this bias can be further reduced by using boosting methods.
Particularly for the Friedman signal, we can see that increasing the maximum
depth in the second boosting step can lead to a higher coverage rate. By
contrast, in the Constant signal, we can see that the initial model has
already achieved a good coverage rate. As a result, boosting with a second
model is unnecessary. The CoE is consistently good across all the settings
indicating that our IJ methods capture the variance well.

###  4.2 Forest Refinements of a Linear Model

An important observation within this paper is that the generality of the IJ
estimate means that we need not restrict ourselves to boosting using only
ensembles. Here we work with the GLM (generalized linear model) class of
models to which we add a random forest to capture remaining signal and compare
the performance of such a boosted model with the base GLM. We use linear model
as a specific example of GLM for the demonstrations.

We consider the exact same setting in section 4.1, except that we use simple
linear regression model for the initial step of boosting. The resulting
coverage is shown in Figure 3. Clearly in case of the Linear and Constant
signals the linear model performs optimally making the boosting step
unnecessary. Also as expected the linear model doesn’t perform well for the
Friedman signal and it’s forest refinement can achieve better coverage. In
section 5, we will demonstrate that we can also look at tests to check whether
boosting gives a significant difference for the model predictions.

Figure 3: CoT (95% level) for boosting linear model with random forests.
Three scenarios of model signals are considered and the number of trees $B$ is
set to be 1000 or 5000. We use simple linear regression model as the initial
estimator at the first step. In each subplot, the x axis represents the model
for the second boosting step. The red boxplot with x axis label ”Without
Boosting” is for the case when there is no second step of boosting. The number
”3,5,7,9” denotes the different values of maximum depth when using random
forests for boosting. And ”Full” represents the random forests with trees
grown to full depth for boosting. For comparison purpose, we also use green
star symbols to present the average CoE with 95% level over 100 query points.
The rest is the same as Figure 2.

###  4.3 Local Model Modifications of Random Forests

The approach of boosting models above can be cast within the framework of
smoothing residuals to test for goodness of fit [hart2013nonparametric], and
this can also be thought of as providing a bias correction by boosting with a
second-stage kernel smooth. We illustrate this procedure with the following
example method.

In lu2021unified the a local linear modification is defined by a bias term is
given by $\widehat{Bias}(x)=\sum_{k=1}^{n}v_{k}(x)(Y_{k}-\hat{f}(X_{k}))$,
where $(Y_{k},X_{k})$ is the $k$th training data, $x$ is a query point and
$v_{k}(x)$ is the out-of-bag in-leaf proportion given by

$v_{k}(x)=\frac{\sum_{b=1}^{B_{n}}\mathds{1}\left\\{Z_{k}\notin I_{b},X_{k}\in
L_{b}(x)\right\\}}{\sum_{\ell=1}^{n}\sum_{b=1}^{B_{n}}\mathds{1}\left\\{Z_{\ell}\notin
I_{b},X_{\ell}\in L_{b}(x)\right\\}}$

in which $X_{k}\in L_{b}(x)$ indicates that $X_{k}$ falls into the same leaf
as $x$ in the $b$th tree, and we weight co-inclusion only over out-of-bag
samples. In that paper, prediction intervals are produced from the
distribution of out-of-bag errors weighted by the $v_{k}$. Here we show that
we use the IJ variance estimates to provide confidence intervals of the bias
corrected estimator. We rewrite
$\widehat{Bias}(x)=\frac{\sum_{k=1}^{n}C_{k}D_{k}}{\sum_{k=1}^{n}C_{k}}$,
where $C_{k}=\sum_{b=1}^{B_{n}}\mathds{1}\left\\{Z_{k}\notin I_{b},X_{k}\in
L_{b}(x)\right\\}$ and $D_{k}=Y_{k}-\hat{f}(X_{k})$. If we denote
$p=\sum_{k=1}^{n}C_{k}D_{k}$, $q=\sum_{k=1}^{n}C_{k}$, $r=nC_{i}D_{i}$ and
$s=nC_{i}$, then following the same calculations in section 3.3, the
directional derivatives of $\widehat{Bias}(x)$ are given by

$\displaystyle U_{i}=\frac{rq-ps}{q^{2}}$

Then the final estimate given by $\hat{f}(x)+\widehat{Bias}(x)$ will have a
variance that is consistently estimated by

$V(x)=\frac{1}{n^{2}}\sum_{i=1}^{n}(U^{\prime}_{i}+U_{i})^{2},\qquad
U^{\prime}_{i}=n\cdot cov_{b}(N_{i,b},T_{b}(x))$

Using this variance estimate we can thus compare the models $\hat{f}$ vs
$\hat{f}+\widehat{Bias}$ and also produce confidence intervals for the final
estimator $\hat{f}+\widehat{Bias}$.

Using the same simulation setup in section 2, we compare the CoT between
original model and the model after modification. The resulting plot is in
Figure 4. We observe that the local linear modifications does not affect
coverage for the Constant signal since the random forest with a simple
structure is already enough to fit the signal. As for Linear and Friedman
signal, we can see that local model modification can efficiently improve the
coverage rate by reducing the bias of original model estimator.

We notice that Figure 2, 3 and 4 share similar patterns. Under the Constant
signal, the original model (without boosting or local linear modifications)
already achieves good CoT. So boosting or model modifications doesn’t change
the coverage much at all. However, under the Friedman signal, the original
model suffers from bias and the boosted or modified model can efficiently
improve the coverage. Under Linear signal, the original linear regression
model of Figure 3 already achieves optimal performance while it still needs
boosting or modifications in the scenarios of Figure 2 and 4. By comparing the
CoT, we conclude that boosting of random forests in Figure 2 can achieve the
best overall performance.

Figure 4: CoT (95% level) for local model modification of random forests.
Three scenarios of model signals are considered and the number of trees $B$ is
set to be 1000 or 5000. In each subplot, the x axis with the number ”3,5,7,9”
denoting the different values of maximum depth of random forests model. And
”Full” represents the random forests with trees grown to full depth. The red
represents the model without modification while blue denotes the random
forests model with local model modification. For comparison purpose, we also
use green star symbols to present the average coverage of the expectation with
95% level over 100 query points.The rest is the same as Figure 2

##  5 Model Comparisons

In this section, we first demonstrate the procedure of constructing the test
statistics for model comparisons. Then we give an example by comparing random
forests with different maximum depths and GLM.

###  5.1 Test Statistics Construction

Suppose we have training data
$\big{(}Z_{i}\big{)}_{i=1}^{n}=\big{(}(Y_{i},X_{i})\big{)}_{i=1}^{n}$ and we
fit two estimators $\hat{f}_{1}$ and $\hat{f}_{2}$. We wish to test if the
predictions from these two estimators are significantly different. Note that
this difference (and its significance) may change with the query point. This
will be demonstrated over a set of query points $x_{1},\dots,x_{m}$;
evaluation of this difference at a single query point is a special case.

Suppose the directional derivatives for $\hat{f}_{1}$ and $\hat{f}_{2}$ are
given by $U_{i}^{(1)}(x)$ and $U_{i}^{(2)}(x)$ respectively for any query
point $x$ and $i=1,\dots,n$. Now consider the vector

$F=\begin{pmatrix}F_{1}\\\ F_{2}\end{pmatrix}\sim N(\mu,\Sigma);\text{
asymptotically, where }F_{p}=\begin{pmatrix}\hat{f}_{p}(x_{1})\\\ \vdots\\\
\hat{f}_{p}(x_{m})\end{pmatrix},\;p=1,2$

here $\Sigma=\begin{pmatrix}\Sigma^{(11)}&\Sigma^{(12)}\\\
\Sigma^{(21)}&\Sigma^{(22)}\end{pmatrix}$ is a $2m\times 2m$ with $m\times m$
blocks. Now $\Sigma^{(pq)}_{ij}$ is the covariance between
$\hat{f}_{p}(x_{i})$ and $\hat{f}_{q}(x_{j})$, $p,q=1,2$, $i,j=1,\dots,m$ and
can be estimated by

$\hat{\Sigma}^{(pq)}_{ij}=\frac{1}{n^{2}}\sum_{k=1}^{n}U_{k}^{(p)}(x_{i})U_{k}^{(q)}(x_{j})$

The consistency of the IJ covariance estimate between two random forests was
shown in ghosal2020boosting. Consistency of the covariance estimate between an
M-estimator and a random forest is shown in appendix C; a special case of this
consistency was shown previously in ghosal2021generalised where the
M-estimator was an MLE.

Now our null hypothesis is that the estimators $\hat{f}_{1}$ and $\hat{f}_{2}$
are not significantly different and thus $\mathds{E}[F_{1}-F_{2}]=0$ under the
null. Further the covariance matrix of $F_{1}-F_{2}$ will be consistently
estimated by
$\hat{\Sigma}=\hat{\Sigma}^{(22)}+\hat{\Sigma}^{(11)}-\hat{\Sigma}^{(12)}-\hat{\Sigma}^{(21)}$.
Thus our test-statistic will be

$(F_{1}-F_{2})^{\top}\hat{\Sigma}^{-1}(F_{1}-F_{2})\sim\chi^{2}_{m}$

where $\chi^{2}_{m}$ represents the chi-squared distribution with $m$ degrees
of freedom. Note that the estimation of
$\tilde{\Sigma}=\Sigma^{(22)}+\Sigma^{(11)}-\Sigma^{(12)}-\Sigma^{(21)}$ with
$\hat{\Sigma}$ involves uncertainty which will depend on $n$ and other details
of fitting the estimators $\hat{f}_{1}$ and $\hat{f}_{2}$ such as number of
trees for fitting random forests. This uncertainty is unaccounted for in our
procedures; we assume that it is small enough so that the null distribution of
the test-statistic above holds asymptotically.

In the case of random forests, note that we can directly use V-statistics
based correction method to obtain unbiased estimates of $\Sigma^{(11)}$ and
$\Sigma^{(22)}$. Then $\Sigma^{(12)}$ and $\Sigma^{(21)}$ are estimated by
using the uncorrected directional derivatives $U_{i}(x)=n\cdot
cov_{b}(N_{i,b},T_{b}(x))$ defined in section 3.2. This covariance estimate is
unbiassed when $f_{1}$ and $f_{2}$ are independent, conditional on the
training set since the Monte Carlo error in the directional derivatives of
each has expectation 0. When $f_{1}$ and $f_{2}$ are both ensembles, they can
alternatively be obtained using the same subsamples, in which case the
combination of corresponding ensemble members can be thought of as a single
ensemble and IJ applied accordingly. ghosal2020boosting found that choosing
subsamples independently between the two models improved predictive
performance we we take this approach throughout.

When a large number of query points are used, we must estimate a large
covariance matrix. Doing so accurately becomes increasingly challenging as the
dimension increases. These errors are compounded when taking the difference
between models and readily result in negative eigenvalues. For this paper we
only use a small $(\leq 5)$ number of query points to ensure a stable and
well-conditioned $\Sigma$, but discuss potential alternatives in section 9.

Note that for boosting models we can check if the boosting provides a
significant difference by comparing the base model and the final estimate (in
fact we can generalise this to comparing any two stages of a boosting model if
we boost for more than two steps). Using the notations in section 4, we wish
to test if the use of $\hat{f}_{2}$ was a significant addition, i.e, if the
predictions from $\hat{f}_{1}$ and $\hat{f}_{1}+\hat{f}_{2}$ are significantly
different.

In this case

$F_{p}=\begin{pmatrix}\hat{g}_{p}(x_{1})\\\ \vdots\\\
\hat{g}_{p}(x_{m})\end{pmatrix},p=1,2\text{ and
}\hat{\Sigma}^{(pq)}_{ij}=\frac{1}{n^{2}}\sum_{k=1}^{n}V_{k}^{(p)}(x_{i})V_{k}^{(q)}(x_{j}),$

where $\hat{g}_{1}=\hat{f}_{1}$, $\hat{g}_{2}=\hat{f}_{1}+\hat{f}_{2}$,
$V_{k}^{(1)}=U_{k}^{(1)}$, $V_{k}^{(2)}=U_{k}^{(1)}+U_{k}^{(2)}$. Thus our
test-statistic can be constructed as

$G^{\top}\left(\hat{\sigma}^{(22)}\right)^{-1}G\sim\chi^{2}_{m}$

where

$\displaystyle G$
$\displaystyle:=F_{2}-F_{1}=\begin{pmatrix}\hat{f}_{2}(x_{1})\\\ \vdots\\\
\hat{f}_{2}(x_{m})\end{pmatrix}$ $\displaystyle\hat{\sigma}^{(22)}_{ij}$
$\displaystyle:=\hat{\Sigma}^{(22)}_{ij}+\hat{\Sigma}^{(11)}_{ij}-\hat{\Sigma}^{(12)}_{ij}-\hat{\Sigma}^{(21)}_{ij}=\frac{1}{n^{2}}\sum_{k=1}^{n}U_{k}^{(2)}(x_{i})U_{k}^{(2)}(x_{j}),$

Similarly as above the estimate $\hat{\sigma}^{(22)}$ has degrees of freedom
which will depend on $n$ and other details of fitting the estimators
$\hat{f}_{1}$ and $\hat{f}_{2}$. These degrees of freedom are difficult to
calculate but we can assume that it is high enough so that the null
distribution of the test-statistic above holds asymptotically.

###  5.2 Comparisons between random forests and a Linear Model

We give an example of comparing random forests and linear model. We consider
the settings described in section 2 and set the number of query points to be
5. Following the procedure in section 5.1, we compare the select pairs of
models and calculate the power, i.e., the proportion of times the null
hypothesis was rejected out of 200 replicates (at a 5% level). These power
values are graphically presented in Figure 5. In each subplot, the numbers
$\\{3,5,7,9\\}$ on the x axis represent random forests with different maximum
depth. And ”Full” represents the random forests with trees grown to full
depth. ”LM” on the x axis represents the simple linear regression. The colors
in each heatmap plot denote the strength of the power. A darker rectangle
between two models means a higher power of rejecting the null hypothesis.

Figure 5: Model Comparisons of multiple random forests and linear model. Three
data generating processes are considered and the number of query points is set
to be 5. The number of trees $B$ is set to be 1000 or 5000. In each subplot,
we draw the heatmap for the power, i.e., the proportion of times the null
hypothesis was rejected out of 200 replicates (at a 5% level) between each
model pair. The numbers $\\{3,5,7,9\\}$ on the x or y axis represent the
random forests with different maximum depth. And ”Full” represents the random
forests with trees grown to full depth. ”LM” on the x axis represents the
linear model and we consider simple linear regression in this example.

We see that for the Friedman signal the linear model is found to be
significantly different from all the random forest models, whereas the effect
is much less pronounced for the Linear and Constant signals. The comparison
between the different forests shows nearly no difference for Constant signal.
This is reasonable since the Constant signal can be easily fitted with any
maximum depth. As for Linear and Friedman signals, we observe a relatively
high power between the value ”3” and ”5,7,9, Full”. However, the difference
seems to be ignorable between the trees forests with higher maximum depth
lower number of trees. The choice of larger number of trees gives us more
significant differences. This is because we can get a more precise estimator
of the covariance with V-statistics based correction. A relatively small
number of trees could possibly lead to a mildly conservative test.

We also give an example to demonstrate the comparisons of boosting models in
Figure 6. We boost linear model with random forests with different maximum
depth and compare the base model and the final boosted estimate to check if
the boosting provides a significant difference. We observe that only under
Friedman Signal, the difference is very significant. However, it shows nearly
no difference for both Linear and Constant Signals.

Figure 6: Model Comparisons for boosting LM with random forests.Three data
generating processes are considered and the number of query points is set to
be 5. The number of trees $B$ is set to be 1000 or 5000. In each subplot, we
draw the heatmap for the power, i.e., the proportion of times the null
hypothesis was rejected out of 200 replicates (at a 5% level) between each
model pair. The numbers $\\{3,5,7,9\\}$ on the x axis represent the random
forests with different maximum depth. And ”Full” represents the random forests
with trees grown to full depth. ”LM” on the y axis represents using linear
model as the base model.The rest is the same as Figure 5.

##  6 IJ and Neural Network

Deep learning methods have recently achieved state-of-the-art performance on a
variety of prediction and learning tasks naturally leading to the need for
uncertainty quantification. Mathematically, we can formalize a neural network
with $L$ hidden layers and an activation function $\sigma$ as:

$\displaystyle\begin{split}\textrm{NN}\left(x;A^{(1)},b^{(1)},\ldots,A^{(L)},b^{(L)}\right)=\;A^{(L)}\sigma\left\\{\ldots
A^{(2)}\sigma\left(A^{(1)}x+b^{(1)}\right)\ldots+b^{(L-1)}\right\\}+b^{(L)},\end{split}$

where the dimension $m_{\ell}$ is the number of nodes at layer $\ell$,
$\ell=0,\ldots,L$, $x\in\mathbb{R}^{m_{0}}$ is the input signal,
$A^{(s)}\in\mathbb{R}^{m_{\ell}\times
m_{\ell-1}},b^{(s)}\in\mathbb{R}^{m_{\ell}}$ are the parameters that produce
the linear transformation of the $(\ell-1)$th layer, and the output is a
scalar with $m_{L}=1$.

The most immediate way to approach uncertainty quantification for neural
networks is to treat them as parametric models falling within the framework of
M-estimators which we explore above. However, as we demonstrate, a number of
technical issues lead IJ-based uncertainty quantification to have poor
statistical performance. Firstly, the structure of neural network is much more
complex than traditional models because of the nonlinear transformation with
activation functions and the high dimensions of the parameter space. In
addition, it is known that the optimization of neural network is quite
challenging [sun2019optimization] and often approached by using stochastic
gradient descent and early stopping [bottou1991stochastic]. Even when they
provide good predictive performance, these methods normally fail to guarantee
a global optimal solution. Thus the complex procedure of training neural
network also introduces another source of procedural variability because of
the random initialization, mini-batch gradient descent and early stopping
[huang2021quantifying]. This makes the uncertainty quantification even more
difficult.

###  6.1 Neural Network as an M-Estimator

(The paragraphs in this section are quite large. Feels like a wall of text.
Maybe we can break them up a bit?)

In this paper, we will focus on a small neural network model with only one
hidden layer and small number of hidden units. Instead of mini-batch gradient
descent, we use the whole training dataset to calculate gradient for
optimization and run enough epochs to ensure convergence. In addition, we also
consider training the network weights from fixed initialization parameters to
remove an important source of variability. Under these conditions, we can
assume that the global optimal point can be attained and there is no source of
procedural variability. Under these restrictions, we can regard a small neural
network as an M-estimator if a smooth (sigmoid) activation function is used.
By using the calculation in section 3.1, we can calculate obtain the IJ
directional derivatives for neural networks. We explore the empirical
performance of this application of IJ below.

####  6.1.1 Experiment Settings

We conduct a simulation experiment with Tensorflow2.0 in python. In our
experiments, we use both ReLU and sigmoid (logistic) activation functions.
Note that ReLU is not a smooth function so it will violate IJ assumptions.
However, as the most commonly used activation function, we will still examine
the application of IJ. We also vary the number of hidden units in the range
$\\{1,3,5,10,20\\}$ to see how the performance changes with the complexity of
the neural network structure. During the process of optimization, we use the
whole training dataset to calculate the gradient and the ”Adam” algorithm for
optimization [kingma2014adam]. The number of epochs is set to be 1000 and the
learning rate is 0.01. In our experimental settings, these were found to
reliably result in solutions that had converged to a local optimum. In these
experiments we observe that the Hessian matrix
$\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta},Z)]$ in
equation (3.1) is generally ill-conditioned because of the numerical issues.
To resolve this, we add a small identity matrix $0.001*I$ to it before taking
the inverse if its condition number is larger than $10^{5}$. We use the same
settings as in section 2. We record both Monto Carlo variance (the sample
variance across 200 repetitions) and the CoE for the predictors across the
query points and generate its boxplots, in order to investigate the
performance of IJ in representing variance and the distribution of neural
networks. The result is presents in Figure 7. Sub-figure (a) is for Monto
Carlo variance while sub-figure (b) records the CoE. At each sub-figure, we
use red to represent the training process with random initialization of neural
network parameters while blue denotes the case when the initialization is
fixed with a specific unchanged random seed over the replications. The numbers
at the x axis represents different number of hidden units of the neural
network structure. Both ReLU and sigmoid activation function are considered.
(The last sentence is a repetition)

(a) Boxplot for Monte Carlo Variance.

(b) Boxplot for CoE.

Figure 7: IJ of neural network with different number of hidden units. Three
data generating processes and both ReLU and sigmoid activation functionare
considered. (a) is for Monto Carlo Variance while (b) records the CoE. At each
sub-figure, we use red to represent the training process with random
initialization while blue is for the case when the initialization is fixed
with a specific unchanged random seed over the replications. The numbers at
the x axis represents different number of hidden units of the neural network
structure.

####  6.1.2 Discussions of Results

We first observe that for both ReLU and sigmoid activation functions, the
Monto Carlo variance increases as the number of hidden units becomes larger
under the Linear and Constant signals. This can be explained as the result of
overfitting. In contrast, under the Friedman signal, if we use random
initialization and sigmoid function, the Monto Carlo variance is extremely
large when the number of hidden units is 1 and actually decreases with the
number of hidden units. However, with fixed initialization, the Monto Carlo
variance is always quite small. This behavior results from the lack of
smoothness of the ReLU function, which results in many local optima. As a
result, even with very careful optimization procedures, it can still fall into
these local optimal points, which introduces an extra source of procedure
variability. When the number of hidden units is increased, the difference of
the predictors corresponding to each local optimal point becomes smaller and
the procedure variability can be thus reduced. Starting the optimization from
fixed initialization weights removes the effect, revealing initialization to
be an important – and difficult to quantify – source of variability.
Unfortunately, there do not appear to be good non-random choices of initial
weights. This can also happen when we consider sigmoid activation functions
under Friedman signal. However, the effect is much smaller in this case, since
the activation function is smoother, and the optimization process
correspondingly less complex.

Examining the boxplots for different signals, it is not surprising to find
that the ReLU activation function can achieve very good coverage for both
Linear and Constant signals, because this activation function preserves the
linear or constant forms of the function resulting in minimal bias. However,
for the Friedman signal, it suffers from the problem of local optimal points
as we have mentioned above when the number of hidden units is small. As a
result, it leads to a serious undercoverage when we consider random
initialization. Then by fixing the initialization or increasing the number of
hidden units, the coverage can be much improved. As for the sigmoid activation
function, it can successfully achieve reasonable coverage under the Friedman
signal but it suffers from undercoverage when the number of hidden units is
increased under Linear and Constant signals. This can be explained by
observing that when the number of hidden units is relatively large, many
hidden nodes will produce outputs very close to zero or one after the
transformation of sigmoid function. This can lead to an ill-conditioned
Hessian matrix, which results in the high inaccuracy of the IJ estimator.

Overall, we only achieve to get a valid IJ estimator for sigmoid function when
the number of hidden units is quite small and thus the model is close to being
parametric. But larger networks result in numerical challenges. For ReLU
activation functions, even though it can achieve promising performance for
Linear and Constant signals, it fails to produce satisfactory CoE for Friedman
signal. We speculate this is a result of the discontinuties in the derivative
of the ReLU activation function violate IJ assumptions. For a small network
the resulting discontinuities affect only a small enough fraction of the data
to be ignored; as the size of the network grows the set of resulting
discontinuities becomes dense in the data, rendering variation in the IJ
tangent plane a poor approximation to sample variance. More work is needed to
quantify the procedure variability introduced by complex neural network
structure and the numerical issue of the ill-conditioned Hessian matrix should
also be further carefully studied.

###  6.2 Model Comparisons with Neural Network

As a final exercise, we compare neural network with multiple random forests. A
small neural network with one hidden layer of 5 units and sigmoid activation
function (for which IJ has reasonable performance) is considered. The result
is presented in Figure 8. Here we see an obvious difference between neural
network and random forests under Friedman signal while both of them are quite
similar under the other two. Here constraining the maximum depth of trees
likely leads to additional bias in the random forest and high power to detect
differences between the two model classes.

Figure 8: Model Comparisons of multiple random forests and neural network. A
small neural network with one hidden layer of 5 units and sigmoid activation
function is used. Three data generating processes are considered and the
number of query points is set to be 5. The number of trees $B$ is set to be
1000 or 5000. In each subplot, we draw the heatmap for the power, i.e., the
proportion of times the null hypothesis was rejected out of 200 replicates (at
a 5% level) between each model pair. The numbers $\\{3,5,7,9\\}$ on the x axis
represent the random forests with different maximum depth. And ”Full”
represents the random forests with trees grown to full depth.The rest is the
same as Figure 5.

##  7 Ensemble Models as V-statistics

Section section 3.2 introduced IJ calculations for random forests. In this
section we observe that the IJ derivation applies to any ensemble obtained
from subsamples taken with replacement. In particular, this represents an
alternative to bootstrap methods for uncertainty quantification; we can employ
an ensemble of models which generally exhibit lower variance, and provide
uncertainty quantification and with little additional cost. Recall that for
random forests model, we use $T_{b}(x)$ to represent the prediction of the
tree at the data point $x$, where this tree is trained on the $b$th subsample
through bagging process. Generally, we can apply this framework to any other
models by denoting $T_{b}(x)$ as the prediction of the model at the data point
$x$, where the model is trained on the $b$th subsample through bagging
process. Following exactly the procedure in section 3.2, we can calculate the
covariance matrix and directional derivatives for uncertainty quantification
or model comparisons. We advocate for this procedure over bootstrapping as it
allows us to take advantage of the ensemble of models to both improve
prediction and provide uncertainty quantification.

Do we use the exact same formula for random forests and XGBoost? For the
bagged random forests $T_{b}$ does not depend on the other trees but for
XGBoost $T_{b}$ depends on all the trees $T_{1},\dots,T_{b-1}$. I feel like
there needs to be a correction for that. $T_{b}$ is probably a good substitute
for $T_{*}$ when independent but not when the dependence is based on the index
$b$. Or is the fact that the bootstrap sample/subsample is independent for
these enough? Sorry if the argument is a bit scatter-brained! :)

(jeff: To clarify it, each $T_{b}$ is a whole XGBoost model. So we borrow the
ideas of random forests to use any machine learning model as $T_{b}$.)

Figure 9: CoE (95% level) for ensemble XGBoost. We consider three data
generating processes. The number of subsamples is set to be 1000 with size
200. We draw the boxplots for the CoE over 100 query points. Figure 10: Model
Comparisons of multiple random forests and ensemble XGBoost. Three data
generating processes are considered and the number of query points is set to
be 5. The number of trees $B$ is set to be 1000 or 5000. In each subplot, we
draw the heatmap for the power, i.e., the proportion of times the null
hypothesis was rejected out of 200 replicates (at a 5% level) between each
model pair. The numbers $\\{3,5,7,9\\}$ on the x axis represent the random
forests with different maximum depth. And ”Full” represents the random forests
with trees grown to full depth.The rest is the same as Figure 5.

In the following, we use an example of XGBoost for demonstrations. XGBoost
(Extreme Gradient Boosting) is a scalable, distributed gradient-boosted
decision tree machine learning method that has gained great popularity and
attention recently in the field of machine learning. We obtain an ”ensemble
XGBoost” model by training XGBoost multiple times on subsamples taken with
replacement, falling within the framework of V-statistics. We use the xgboost
in python for implementation with its default settings for hyperparameters
(please mention these defaults, like number of tress, depth of trees,
subsample fraction, etc). The number of trees is 100 and the maximum depth is
6. The ensemble is generated from 1000 subsamples, each of size 200. The
simulation settings in section 2 (Incorrect reference. I’ve corected some but
there are many more in throughout) were employed to calculate CoE of ensemble
XGBoost and compare it with multiple random forests models. We show the result
of CoE in Figure 9, which demonstrates a very good performance for the
coverage. In Figure 10, we observe that for the Linear and Constant signals,
the difference between ensemble XGBoost and random forests is much less
pronounced. For the Friedman signal, the difference is quite obvious,
especially when the maximum depth is small. We can expect that as the maximum
depth increases, random forest will tend to behave similarily to ensemble
XGBoost.

##  8 Real Data Analysis

We study the practical application of our methods to Beijing Housing Data
[house]. These data contain the price of Beijing houses from 2011 to 2017,
(scraped from www.lianjia.com), along with covariates including location,
size, building features, age and configuration. We prepossess pre-process the
data by removing the rows with missing values and dropping the columns of url,
transaction and community id, and the average price by square (since this can
be calculated by the other two columns: total price and the square of house).
The resulting dataset has dimension $15360\times 51$ after converting all
categorical features to dummy variables. We take logarithm of the total price
and use it as the outcome to be predicted from the remaining covariates. In
order to explore the performance of the coverage of model predictions, we
randomly split the whole data into 50 disjoint training sets
$\\{\mathcal{D}_{i}\\}_{i=1}^{50}$ each of size 3000 and and use the remaining
3260 as a test set (Just to clarify: Do we do this training just once? Or are
we doing cross-validation? Saying ”50 different training sets” feels like
we’re doing CV. But in CV there is usually some hyperparameter optimisation so
we should mention what we’re optimising.). (jeff: We didn’t do any cross-
validation here. We create 50 training sets only for evaluating the
reproduction confidence interval.) We randomly draw 100 samples from the test
set as the query points at which to evaluate the coverage of model
predictions.

For this section, we consider ”reproduction intervals” as produced in
zhou2018boulevard. Rather than being defined to cover a population quantity,
these are defined to cover the prediction that would be made by applying the
same learning procedure on an independent sample. That is, for
$\hat{f}_{1}(x),\hat{f}_{2}(x)$ independently
$N(\mathds{E}[\hat{f}_{i}(x)],V(x))$ (What are $\hat{f}_{1}(x),\hat{f}_{2}$
here? I don’t see any model comparisons in the pictures - the models are
evaluated separately. And why do they have the same variance?) then
$\hat{f}_{1}(x)-\hat{f}_{2}(x)\sim N(0,2V(x))$. Thus

$I(x)=\left[\hat{f}(x)-1.96\sqrt{2\hat{V}(x)},\hat{f}(x)+1.96\sqrt{2\hat{V}(x)}\right]$

For our assessment, at each query point we construct a reproduction interval
$I_{i}(x)$ for each $D_{i}$. We then measure the proportion of replicates
falling within $I_{i}(x)$ and report the Coverage of Reproduction (CoR) by
averaging these proportions:

$\mbox{CoR}(x)=\frac{1}{50*49}\sum_{i=1}^{50}\sum_{j\neq i}(\hat{f}_{j}(x)\in
I_{i}(x)).$

We also evaluate the mean square error (MSE) for each estimated model
$\hat{f}_{i}$ on the testing set, denoted by $e_{i}$. Then we report the
average MSE $\frac{1}{50}\sum_{i=1}^{50}e_{i}$ of all the models fitted on
each training segment.

We consider five different models: (1) LM (linear model), (2) LM+RF (boost
linear model with a random forest), (3) RF (random forest) , (4) RF+RF (boost
a random forest with another random forest) (5) XGB (XGBoost). We use random
forests with full depth and 1000 trees here and report the boxplot of CoR and
MSE of each model. The result is presented in Figure 11(a). We can observe
that all these models achieve sufficient coverage and the models that involve
random forests tend to be more conservative for CoR. As a byproduct, for a
fixed query point, we plot the confidence intervals of predictions by the
models fitted on one specific segment of the training data shown in Figure
11(b). In addition, we also implement the model comparisons test for each pair
of model predictions and use the horizontal lines to group the pairs of model
predictions which are not statistically distinguishable (if possible make this
to look like a dendogram. Right now it’s not really clear. Or you could make a
heatmap like the previous model comparisons.). We observe that the prediction
of random forests is relatively similar with XGB, boosted random forests and
modified random forests but that linear models do appear to give a
statistically distinct prediction, even when later modified by a random
forest.

(a) CoR and MSE

(b) Model Comparisons

Figure 11: Analysis of Beijing Housing Data. (a) is the boxplot of CoR and MSE
of each model. x-axis repsents five different model combinations: LM (linear
model), LM+RF (boost linear model with a random forest), RF (random forest) ,
RF+RF (boost a random forest with another random forest), and XGB (XGBoost).
(The models are described in the text so no need to repeat here. Also I prefer
if the MSE is in a separate plot - or atleast make it’s corresponding y-axis
tick marks the same blue color.) (b) is the confidence intervals of
predictions for these different model combinations. The bottom horizontal
lines group the pairs of model predictions which are not statistically
distinguishable.

##  9 Discussions

In this paper we revisited the classical but underused IJ method, first
described in efron1982jackknife to derive variance estimates. The underlying
IJ directional derivatives are highly versatile and can be used to construct
estimates of the covariance between two different models trained on the same
data. With this we devised a statistical test to quantify how exactly a model
differs from another. In the case of ensemble methods, IJ is biased upwards.
However bias corrections can be successfully derived for ensembles obtained by
subsampling with replacement.

We demonstrated that IJ is an efficient method for model comparisons or
providing uncertainty estimates for combinations of models in the case of
random forests or M estimators, where the specific expressions for directional
derivatives are derived. However, as for more general machine learning models,
we propose to emsemble them as a V-statistic, which provides a generic means
of producing uncertainty quantification.

There are also some ideas related to the content of this paper that could be
explored in the future. Even though we only worked with Gaussian responses in
this paper, the methods can be readily applied to non-Gaussian responses from
an exponential family. The software corresponding to this paper can work with
such data, where the random forest “boosts” are fitted on pseudo-residuals as
defined in ghosal2021generalised.

Closed forms for the IJ could also be explored for any general V-statistic or
more general ensemble methods. For any two given models proving consistency of
the IJ covariance estimate (or disproving it or constructing suitable
regularity conditions) is also a theoretical challenge.

Note that we have ignored the uncertainty in the estimate $\hat{\Sigma}$
throughout. If there are a large number of query points, IJ estimates of the
covariance matrix $\Sigma$ could be very unstable. A potential future
direction would be to characterize the distribution of this estimate,
especially for random forests covariances, and modify our uncertainty
quantification correspondingly.

We could also extend the boosting framework discussed in section 4 to work
with ensembles where the contribution from each boosting stage may be random
or vary during the training itself - for example the Boulevard boosting model
where the weight for each tree depends on the eventual total (random) number
of trees [zhou2018boulevard] and it keeps changing during training as more
trees are added to the boosting ensemble.

Finally, we also need to resolve the problem of calculating IJ directional
derivatives of neural networks. Our simulations show many challenges to a
naïve framework for this, including the technical issues of optimization,
numerical issues of sigmoid activation functions in larger number of hidden
units. Overcoming these requires either a more direct calculation of
directional derivatives, possibly through automatic differentiation, or a new
representation of neural network training.

##  Appendix A Asymptotics for M-estimators

##### The Jacobian

Suppose $f:\mathds{R}^{m+n}\to\mathds{R}^{\ell}$ is a differentiable function.
Then for any $(x_{0},y_{0})\in\mathds{R}^{m+n}$ (where
$x_{0}\in\mathds{R}^{m},y_{0}\in\mathds{R}^{n}$) the (partial) Jacobian is a
matrix [in $\mathds{R}^{\ell\times m}$] defined by

$\displaystyle J_{f,x}(x_{0},y_{0})$
$\displaystyle=\begin{pmatrix}\frac{\partial f}{\partial
x_{1}}(x_{0},y_{0})\cdots\frac{\partial f}{\partial
x_{m}}(x_{0},y_{0})\end{pmatrix}$
$\displaystyle=\begin{pmatrix}(\nabla_{x}^{\top}f_{1})(x_{0},y_{0})\\\
\vdots\\\
(\nabla_{x}^{\top}f_{\ell})(x_{0},y_{0})\end{pmatrix}=\left(\left(\frac{\partial
f_{i}}{\partial
x_{j}}(x_{0},y_{0})\right)\right)_{\begin{subarray}{c}i=1,\dots,\ell\\\
j=1,\dots,m\end{subarray}}$

##### The Delta Method

Suppose $X,\\{X_{n}\\}_{n=1}^{\infty}$ are random variables in
$\mathds{R}^{m}$, $c\in\mathds{R}^{m}$ is constant and $r_{n}\to\infty$ is a
sequence such that $r_{n}(X_{n}-c)\xrightarrow{d}X$. Also
$f:\mathds{R}^{m}\to\mathds{R}^{\ell}$ is a differentiable function with
Jacobian $J_{g}$. Then by the Taylor expansion

$\displaystyle f(X_{n})$ $\displaystyle=f(c)+J_{f}(\tilde{c})(X_{n}-c),$ where
$\displaystyle\|\tilde{c}-c\|\leq\|X_{n}-c\|$ $\displaystyle\implies
r_{n}(f(X_{n})-f(c))$ $\displaystyle=J_{f}(\tilde{c})\cdot r_{n}(X_{n}-c)$

Now by Slutsky’s theorem
$X_{n}-c=r_{n}(X_{n}-c)\cdot\frac{1}{r_{n}}\xrightarrow{c}X\cdot 0=0$, i.e,
$X_{n}\xrightarrow{d}c\implies\tilde{c}\xrightarrow{d}c\implies
J_{f}(\tilde{c})\xrightarrow{d}J_{f}(c)$. Applying Slutsky’s theorem again we
see that

$\displaystyle J_{f}(\tilde{c})\cdot r_{n}(X_{n}-c)$
$\displaystyle\xrightarrow{d}J_{f}(c)X$ $\displaystyle r_{n}(f(X_{n})-f(c))$
$\displaystyle\xrightarrow{d}J_{f}(c)X$

Now suppose we define the M-estimator to be
$\hat{\theta}=\operatorname*{arg\,max}_{\theta\in\Theta}\mathds{E}_{\hat{\mathscr{D}}}[m(\theta,Z)]$
where $\Theta\subseteq\mathds{R}^{p}$ is the parameter space and the data
$\\{Z_{i}\\}_{i=1}^{n}$ are independently generated from the distribution
$\mathscr{D}$ over $\mathscr{Z}\subseteq\mathds{R}^{d}$. We assume $m$ has a
continuous second derivative with respect to $\theta$. Also
$\hat{\mathscr{D}}$ is the empirical distribution over the training data
$(Z_{i})_{i=1}^{n}$ from $\mathscr{Z}^{n}$, thereby making
$\mathds{E}_{\hat{\mathscr{D}}}[m(\theta,Z)]=\frac{1}{n}\sum_{i=1}^{n}m(\theta,Z_{i})$.
Finally let
$\theta_{0}=\operatorname*{arg\,max}_{\theta\in\Theta}\mathds{E}_{\mathscr{D}}[m(\theta,Z)]$
be the real value of the M-estimate.

##### Asymptotic Consistency

Assume the following conditions

$\displaystyle\sup_{\theta\in\Theta}\left|\mathds{E}_{\hat{\mathscr{D}}}[m(\theta,Z)]-\mathds{E}_{\mathscr{D}}[m(\theta,Z)]\right|\xrightarrow{p}0$
$\displaystyle\sup_{\theta:\|\theta-\theta_{0}\|>\epsilon}\mathds{E}_{\mathscr{D}}[m(\theta,Z)]<\mathds{E}_{\mathscr{D}}[m(\theta_{0},Z)],\;\forall\,\epsilon>0$
$\displaystyle\mathds{E}_{\hat{\mathscr{D}}}[m(\hat{\theta},Z)]\geq\mathds{E}_{\hat{\mathscr{D}}}[m(\theta_{0},Z)]-o_{P}(1)$

The first assumption states that the law of large numbers is uniform over
$\Theta$ and the second one states that the global extremum $\theta_{0}$ is
well-separated. Then from the second assumption $\forall\,\epsilon>0$,
$\exists\,\delta>0$ such that

$\displaystyle P(\|\hat{\theta}-\theta_{0}\|\geq\epsilon)$ $\displaystyle\leq
P(\mathds{E}_{\mathscr{D}}[m(\theta_{0},Z)]-\mathds{E}_{\mathscr{D}}[m(\hat{\theta},Z)]>\delta)$
$\displaystyle=P\left((\mathds{E}_{\mathscr{D}}[m(\theta_{0},Z)]-\mathds{E}_{\hat{\mathscr{D}}}[m(\theta_{0},Z)])+(\mathds{E}_{\hat{\mathscr{D}}}[m(\theta_{0},Z)]-\mathds{E}_{\hat{\mathscr{D}}}[m(\hat{\theta},Z)])\right.$
$\displaystyle\qquad\qquad\left.+(\mathds{E}_{\hat{\mathscr{D}}}[m(\hat{\theta},Z)]-\mathds{E}_{\mathscr{D}}[m(\hat{\theta},Z)])\geq\delta\right)$
$\displaystyle\leq
P(\mathds{E}_{\mathscr{D}}[m(\theta_{0},Z)]-\mathds{E}_{\hat{\mathscr{D}}}[m(\theta_{0},Z)]\geq\delta/3)+P(\mathds{E}_{\hat{\mathscr{D}}}[m(\theta_{0},Z)]-\mathds{E}_{\hat{\mathscr{D}}}[m(\hat{\theta},Z)]\geq\delta/3)$
$\displaystyle\qquad\qquad+P(\mathds{E}_{\hat{\mathscr{D}}}[m(\hat{\theta},Z)]-\mathds{E}_{\mathscr{D}}[m(\hat{\theta},Z)]\geq\delta/3)$

Now it is easily seen that by the first and third terms above has a limit of 0
(by the first assumption) and the second term also has a limit of 0 (by the
third assumption). Thus $\forall\,\epsilon>0$,
$P(\|\hat{\theta}-\theta_{0}\|\geq\epsilon)\to 0$, i.e.,
$\hat{\theta}\to\theta$ is probability and thus almost everywhere.

We can also do a similar proof of consistency for the corresponding
Z-estimators. Define $g(\theta,Z)=\nabla_{\theta}m(\theta,Z)$ and suppose
$\hat{\theta}$ is an (approximate) solution of
$\mathds{E}_{\hat{\mathscr{D}}}[g(\theta,Z)]=0$ and $\theta_{0}$ is an (exact)
solution of $\mathds{E}_{\mathscr{D}}[g(\theta,Z)]=0$ satisfying the following
conditions

$\displaystyle\sup_{\theta\in\Theta}\|\mathds{E}_{\hat{\mathscr{D}}}[g(\theta,Z)]-\mathds{E}_{\mathscr{D}}[g(\theta,Z)]\|\xrightarrow{p}0$
$\displaystyle\sup_{\theta:\|\theta-\theta_{0}\|>\epsilon}\|\mathds{E}_{\mathscr{D}}[g(\theta,Z)]\|>0,\;\forall\,\epsilon>0$
$\displaystyle\|\mathds{E}_{\hat{\mathscr{D}}}[g(\hat{\theta},Z)]-\mathds{E}_{\hat{\mathscr{D}}}[g(\theta_{0},Z)]\|\xrightarrow{p}0$

Then by the second assumption $\forall\,\epsilon>0$, $\exists\,\delta>0$ such
that

$\displaystyle P(\|\hat{\theta}-\theta_{0}\|\geq\epsilon)$ $\displaystyle\leq
P(\|\mathds{E}_{\mathscr{D}}[g(\hat{\theta},Z)]\|>\delta)$
$\displaystyle=P\left((\|\mathds{E}_{\mathscr{D}}[g(\hat{\theta},Z)]\|-\|\mathds{E}_{\hat{\mathscr{D}}}[g(\hat{\theta},Z)]\|)+(\|\mathds{E}_{\hat{\mathscr{D}}}[g(\hat{\theta},Z)]\|-\|\mathds{E}_{\hat{\mathscr{D}}}[g(\theta_{0},Z)]\|)\right.$
$\displaystyle\qquad\qquad\left.+(\|\mathds{E}_{\hat{\mathscr{D}}}[g(\theta_{0},Z)]\|-\|\mathds{E}_{\mathscr{D}}[g(\theta_{0},Z)]\|)\geq\delta\right)$
$\displaystyle\leq
P(\|\mathds{E}_{\mathscr{D}}[g(\hat{\theta},Z)]-\mathds{E}_{\hat{\mathscr{D}}}[g(\hat{\theta},Z)]\|\geq\delta/3)+P(\|\mathds{E}_{\hat{\mathscr{D}}}[g(\hat{\theta},Z)]-\mathds{E}_{\hat{\mathscr{D}}}[g(\theta_{0},Z)]\|\geq\delta/3)$
$\displaystyle+P(\|\mathds{E}_{\hat{\mathscr{D}}}[g(\theta_{0},Z)]-\mathds{E}_{\mathscr{D}}[g(\theta_{0},Z)]\|\geq\delta/3)\qquad[\because\|a\|-\|b\|\geq
x\implies\|a-b\|\geq x]$

Now it is easily seen that by the first and third terms above has a limit of 0
(by the first assumption) and the second term also has a limit of 0 (by the
third assumption). Thus $\forall\,\epsilon>0$,
$P(\|\hat{\theta}-\theta_{0}\|\geq\epsilon)\to 0$, i.e.,
$\hat{\theta}\to\theta$ is probability and thus almost everywhere.

Note that if we assume $\hat{\theta}$ is always an exact solution of
$\mathds{E}_{\hat{\mathscr{D}}}[g(\theta,Z)]=0$ then we can drop the third
assumption and the proof is also simplified as follows - by the second
assumption $\forall\,\epsilon>0$, $\exists\,\delta>0$ such that

$\displaystyle P(\|\hat{\theta}-\theta_{0}\|\geq\epsilon)$ $\displaystyle\leq
P(\|\mathds{E}_{\mathscr{D}}[g(\hat{\theta},Z)]\|>\delta)$
$\displaystyle=P\left(\|\mathds{E}_{\mathscr{D}}[g(\hat{\theta},Z)]\|-\|\mathds{E}_{\mathscr{D}}[g(\theta_{0},Z)]\|)\geq\delta\right)$
$\displaystyle\leq
P(\|\mathds{E}_{\mathscr{D}}[g(\hat{\theta},Z)]-\mathds{E}_{\mathscr{D}}[g(\theta_{0},Z)]\|\geq\delta)\qquad[\because\|a\|-\|b\|\geq
x\implies\|a-b\|\geq x]$

and this term has a limit of 0 by the first condition.

##### Asymptotic Normality

We know that
$\hat{\theta}=\operatorname*{arg\,max}_{\theta\in\Theta}\mathds{E}_{\hat{\mathscr{D}}}[m(\theta,Z)]$
also satisfies that
$\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}m(\theta,Z)]=0$. Thus

$\displaystyle 0$
$\displaystyle=\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}m(\hat{\theta},Z)]$
$\displaystyle=\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}m(\theta_{0},Z)]+\mathds{E}_{\hat{\mathscr{D}}}\left[J_{\nabla_{\theta}m,\theta}(\tilde{\theta},Z)\right](\hat{\theta}-\theta_{0})$
$\displaystyle\text{where
}\|\tilde{\theta}-\theta_{0}\|\leq\|\hat{\theta}-\theta_{0}\|$
$\displaystyle\implies\hat{\theta}-\theta_{0}$
$\displaystyle=-\left(\mathds{E}_{\hat{\mathscr{D}}}\left[J_{\nabla_{\theta}m,\theta}(\tilde{\theta},Z)\right]\right)^{-1}\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}m(\theta_{0},Z)]$

Now $J_{\nabla_{\theta}m,\theta}=\nabla_{\theta}^{2}m$. Also

$\displaystyle\mathds{E}_{\hat{\mathscr{D}}}\left[\nabla_{\theta}^{2}m(\tilde{\theta},Z)\right]$
$\displaystyle\xrightarrow{\text{a.e.}}\mathds{E}_{\mathscr{D}}\left[\nabla_{\theta}^{2}m(\tilde{\theta},Z)\right]\qquad[\text{by
the uniform strong law of large numbers}]$
$\displaystyle\xrightarrow{\text{a.e.}}\mathds{E}_{\mathscr{D}}\left[\nabla_{\theta}^{2}m(\theta_{0},Z)\right]\qquad[\text{by
continuity of }\nabla_{\theta}^{2}m]$

And by the central limit theorem

$\displaystyle\sqrt{n}\left(\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}m(\theta_{0},Z)]-\mathds{E}_{\mathscr{D}}[\nabla_{\theta}m(\theta_{0},Z)]\right)\sim
N_{p}(0,Var_{\mathscr{D}}[\nabla_{\theta}m(\theta_{0},Z)])$

where $\mathds{E}_{\mathscr{D}}[\nabla_{\theta}m(\theta_{0},Z)]=0$ and
$Var_{\mathscr{D}}[\nabla_{\theta}m(\theta_{0},Z)]$ is a $p\times p$ matrix
whose elements are

$\displaystyle\left[Var_{\mathscr{D}}[\nabla_{\theta}m(\theta_{0},Z)]\right]_{ij}$
$\displaystyle=Cov_{\mathscr{D}}\left[\frac{\partial
m}{\partial\theta_{i}}(\theta_{0},Z),\frac{\partial
m}{\partial\theta_{j}}(\theta_{0},Z)\right]$
$\displaystyle=\mathds{E}_{\mathscr{D}}\left[\left(\frac{\partial
m}{\partial\theta_{i}}\cdot\frac{\partial
m}{\partial\theta_{j}}\right)(\theta_{0},Z)\right],\;\;i,j=1,\dots,p$
$\displaystyle\implies Var_{\mathscr{D}}[\nabla_{\theta}m(\theta_{0},Z)]$
$\displaystyle=\mathds{E}_{\mathscr{D}}[(\nabla_{\theta}m\nabla_{\theta}^{\top}m)(\theta_{0},Z)]$

Finally we can apply Slutsky’s theorem to obtain

$\displaystyle\sqrt{n}(\hat{\theta}-\theta_{0})$
$\displaystyle\xrightarrow{d}N_{p}(0,\Sigma(\theta_{0})),$ where
$\displaystyle\Sigma(\theta)=\left(\mathds{E}_{\mathscr{D}}\left[\nabla_{\theta}^{2}m(\theta,Z)\right]\right)^{-1}\mathds{E}_{\mathscr{D}}[(\nabla_{\theta}m\nabla_{\theta}^{\top}m)(\theta,Z)]\left(\mathds{E}_{\mathscr{D}}\left[\nabla_{\theta}^{2}m(\theta,Z)\right]\right)^{-1}$

If the prediction function corresponding to any parameter $\theta\in\Theta$
and query point $x\in\mathscr{Z}$ is $\eta(\theta,x)$ then by the delta method

$\sqrt{n}(\eta(\hat{\theta},x)-\eta(\theta_{0},x))\xrightarrow{d}N_{p}(0,V(x)),\text{
where
}V(x)=\nabla_{\theta}^{\top}\eta(\theta_{0},x)\Sigma(\theta_{0})\nabla_{\theta}\eta(\theta_{0},x)$

##  Appendix B Directional derivatives for M-estimators

Given the training dataset $(Z_{i})_{i=1}^{n}$ and any probability vector $P$
in $\mathds{R}^{n}$ define

$\mathds{E}_{\hat{\mathscr{D}}(P)}[m(\theta,Z)]=\sum_{k=1}^{n}P_{k}m(\theta,Z_{k})$

Now for some $i\in\\{1,\dots,n\\}$ and $\epsilon>0$ if $P=P(i,\epsilon)$ is
given by $P_{k}=(1-\epsilon)\frac{1}{n}+\epsilon\mathds{1}_{k=i}$ then define

$\hat{\theta}(i,\epsilon)=\operatorname*{arg\,max}_{\theta\in\Theta}\mathds{E}_{\hat{\mathscr{D}}(P(i,\epsilon))}[m(\theta,Z)]$

Thus $\hat{\theta}(i,\epsilon)$ satisfies
$\mathds{E}_{\hat{\mathscr{D}}(P(i,\epsilon))}[\nabla_{\theta}m(\theta,Z)]=0$.
Fixing $i$, if we can show that $\hat{\theta}(i,\epsilon)$ has a derivative at
$\theta=0$ then
$U_{\hat{\theta},i}=\frac{\partial\hat{\theta}}{\partial\epsilon}(i,0)$ will
be the directional derivative for the M-estimator. To show this suppose
$f:\Theta\times\mathds{R}_{+}\to\mathds{R}^{p}$ is given by

$\displaystyle f(\theta,\epsilon)$
$\displaystyle=\mathds{E}_{\hat{\mathscr{D}}(P(i,\epsilon))}[\nabla_{\theta}m(\theta,Z)]=\sum_{k=1}^{n}P_{k}\nabla_{\theta}m(\theta,Z_{k})$
$\displaystyle=\sum_{k=1}^{n}\left((1-\epsilon)\frac{1}{n}+\epsilon\mathds{1}_{k=i}\right)\nabla_{\theta}m(\theta,Z_{k})$

Now we know that $f(\hat{\theta},0)=0$ and we will assume that the Jacobian at
$(\hat{\theta},0)$, given by
$J_{f,\theta}(\hat{\theta},0)=\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta})]$,
is invertible. Then by the implicit function theorem there exists an open set
$S\ni 0$ and a continuously differentiable function
$h:\mathds{R}_{+}\to\mathds{R}^{p}$ such that $f(h(\epsilon),\epsilon)=0$ for
all $\epsilon\in S$. So we can express $\hat{\theta}(i,\epsilon)=h(\epsilon)$
making $h(0)=\hat{\theta}$. Further the derivative of $h$ is given by
$\frac{\partial
h}{\partial\epsilon}(\epsilon)=-\left(J_{f,\theta}(h(\epsilon),\epsilon)\right)^{-1}\frac{\partial
f}{\partial\epsilon}(h(\epsilon),\epsilon)$ over $S$. Hence when $\epsilon=0$,
we get the directional derivative to be

$\displaystyle U_{\hat{\theta},i}$
$\displaystyle=\frac{\partial\hat{\theta}}{\partial\epsilon}(i,0)=\frac{\partial
h}{\partial\epsilon}(0)=-\left(J_{f,\theta}(h(0),0)\right)^{-1}\frac{\partial
f}{\partial\epsilon}(h(0),0)$
$\displaystyle=-\left(J_{f,\theta}(\hat{\theta},0)\right)^{-1}\frac{\partial
f}{\partial\epsilon}(\hat{\theta},0)=-\left(\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta})]\right)^{-1}\frac{\partial
f}{\partial\epsilon}(\hat{\theta},0)$
$\displaystyle=-\left(\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta})]\right)^{-1}\left(\sum_{k=1}^{n}\left(\mathds{1}_{k=i}-\frac{1}{n}\right)\nabla_{\theta}m(\hat{\theta},Z_{k})\right)$
$\displaystyle=-\left(\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta})]\right)^{-1}\left(\nabla_{\theta}m(\hat{\theta},Z_{i})-\mathds{E}_{\hat{\mathscr{D}}}\left[\nabla_{\theta}m(\hat{\theta},Z)\right]\right)$
$\displaystyle=-\left(\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta})]\right)^{-1}\nabla_{\theta}m(\hat{\theta},Z_{i})$

Finally if the prediction function corresponding to any parameter
$\theta\in\Theta$ and query point $x\in\mathscr{Z}$ is $\eta(\theta,x)$ then
by the delta method we define its directional derivative to be

$U_{i}(x)=\nabla_{\theta}^{\top}\eta(\hat{\theta},x)U_{\hat{\theta},i}=-\nabla_{\theta}^{\top}\eta(\hat{\theta},x)\left(\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta})]\right)^{-1}\nabla_{\theta}m(\hat{\theta},Z_{i})$

##  Appendix C Consistency of the IJ

##### Consistency of the IJ variance estimate

We see that the IJ variance estimate for $\sqrt{n}\cdot\eta(\hat{\theta},x)$
is

$\displaystyle n\cdot\frac{1}{n^{2}}\sum_{i=1}^{n}U_{i}(x)^{2}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\nabla_{\theta}^{\top}\eta(\hat{\theta},x)U_{\hat{\theta},i}U_{\hat{\theta},i}^{\top}\nabla_{\theta}\eta(\hat{\theta},x)$
$\displaystyle=\frac{1}{n}\nabla_{\theta}^{\top}\eta(\hat{\theta},x)\left(\sum_{i=1}^{n}U_{\hat{\theta},i}U_{\hat{\theta},i}^{\top}\right)\nabla_{\theta}\eta(\hat{\theta},x)$
$\displaystyle=\frac{1}{n}\nabla_{\theta}^{\top}\eta(\hat{\theta},x)\left(\sum_{i=1}^{n}\left(\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta})]\right)^{-1}\nabla_{\theta}m(\hat{\theta},Z_{i})\nabla_{\theta}^{\top}m(\hat{\theta},Z_{i})\left(\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta})]\right)^{-1}\right)\nabla_{\theta}\eta(\hat{\theta},x)$
$\displaystyle=\nabla_{\theta}^{\top}\eta(\hat{\theta},x)\left(\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta})]\right)^{-1}\left(\frac{1}{n}\sum_{i=1}^{n}\nabla_{\theta}m(\hat{\theta},Z_{i})\nabla_{\theta}^{\top}m(\hat{\theta},Z_{i})\right)\left(\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta})]\right)^{-1}\nabla_{\theta}\eta(\hat{\theta},x)$
$\displaystyle=\nabla_{\theta}^{\top}\eta(\hat{\theta},x)\left(\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta})]\right)^{-1}\mathds{E}_{\hat{\mathscr{D}}}[(\nabla_{\theta}m\nabla_{\theta}^{\top}m)(\hat{\theta},Z)]\left(\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}^{2}m(\hat{\theta})]\right)^{-1}\nabla_{\theta}\eta(\hat{\theta},x)$

Hence by the strong law of large numbers and that
$\hat{\theta}\xrightarrow{\text{a.e.}}\theta_{0}$ it immediately follows that
the IJ variance estimate for M-estimator predictions is consistent.

##### Consistency of the IJ covariance estimate

Suppose we train an M-estimator $\hat{\theta}$ and a random forest $\hat{F}$
on the same training data - then we look at its predictions at two different
points. Suppose these predictions are given by $\eta(\hat{\theta},x_{1})$ and
$\hat{F}(x_{2})$ respectively. Below we show that the IJ covariance estimate
for these predictions are consistent.

Simplifying the M-estimator prediction: Note that by Taylor expansion

$\eta(\hat{\theta},x_{1})=\eta(\theta_{0},x_{!})+\nabla_{\theta}^{\top}\eta(\tilde{\theta},x_{1})(\hat{\theta}-\theta_{0}),\text{
where }\|\tilde{\theta}-\theta_{0}\|\leq\|\hat{\theta}-\theta_{0}\|$

But
$\hat{\theta}\xrightarrow{\text{a.e.}}\theta_{0}\implies\tilde{\theta}\xrightarrow{\text{a.e.}}\theta_{0}$.
Also from the proof of the asymptotic normality for $\hat{\theta}$ we saw that

$\displaystyle\hat{\theta}-\theta_{0}$
$\displaystyle-\left(-\left(\mathds{E}_{\mathscr{D}}\left[\nabla_{\theta}^{2}m(\theta_{0},Z)\right]\right)^{-1}\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}m(\theta_{0},Z)]\right)\xrightarrow{\text{a.e.}}0$
$\displaystyle\implies\eta(\hat{\theta},x_{1})$
$\displaystyle-\left(\eta(\theta_{0},x_{1})+v_{\mathscr{D}}(x_{1})^{\top}\mathds{E}_{\hat{\mathscr{D}}}[\nabla_{\theta}m(\theta_{0},Z)]\right)\xrightarrow{\text{a.e.}}0,$
where $\displaystyle
v_{\mathscr{D}}(x)=-\nabla_{\theta}^{\top}\eta(\theta_{0},x)\left(\mathds{E}_{\mathscr{D}}\left[\nabla_{\theta}^{2}m(\theta_{0},Z)\right]\right)^{-1}$

Simplifying the random forest prediction: Now consider the tree kernel $T$ of
the random forest and its first Hajek projection given by
$\mathring{T}(x;Z_{1},\dots,Z_{k})=\mathds{E}[T]+\sum_{j=1}^{k}T_{1}(x;Z_{j})$.
Then the first Hajek projection of the random forest is given by

$\mathring{F}(x_{2})=\mathds{E}[T]+\frac{k}{n}\sum_{j=1}^{n}T_{1}(x_{2};Z_{j}).$

Now suppose the subsample size of the training data for each tree $k=k_{n}$
varies with $n$ such that

$\displaystyle\lim_{n\to\infty}$ $\displaystyle\frac{k_{n}}{n}\to 0\text{ and
for all query points
}x,\lim_{n\to\infty}\frac{k_{n}\zeta_{1,k_{n}}(x)}{\zeta_{k_{n},k_{n}}(x)}\neq
0,$ $\displaystyle\text{ where }\zeta_{c,k}(x)$
$\displaystyle=cov(T(x;Z_{1},\dots,Z_{c},Z_{c+1},\dots,Z_{k}),T(x;Z_{1},\dots,Z_{c},Z^{\prime}_{c+1},\dots,Z^{\prime}_{k}))$
$\displaystyle=var(\mathds{E}[T(x;Z_{1},\dots,Z_{k})\mid
Z_{1}=z_{1},\dots,Z_{c}=z_{c}])$

Then this condition, alongwith Lemma 13 in wager2017estimation, shows that we
can replace $\hat{F}$ with $\mathring{F}$ without affecting the variance of
the random forest. It can also be easily extended to the covariance of the
random forest with any other well-behaved model. Note that this condition was
also used in Theorem 1 of ghosal2020boosting. We can also use a slightly
weaker condition

$\lim_{n\to\infty}\frac{k_{n}}{n}\cdot\frac{\zeta_{k_{n},k_{n}}(x)}{k_{n}\zeta_{1,k_{n}}(x)}=0,$

as seen in Theorem 1 of peng2019asymptotic.

Simplifying the theoretical covariance: Since $\hat{F}$ is replaceable with
$\mathring{F}$, it follows that

$\displaystyle cov_{\mathscr{D}}(\eta(\hat{\theta},x_{1}),\hat{F}(x_{2}))$
$\displaystyle=cov_{\mathscr{D}}(\eta(\hat{\theta},x_{1}),\mathring{F}(x_{2}))$
$\displaystyle=cov_{\mathscr{D}}\left(v_{\mathscr{D}}(x_{1})^{\top}\left[\frac{1}{n}\sum_{j=1}^{n}\nabla_{\theta}m(\theta_{0},Z_{j})\right],\frac{k}{n}\sum_{j=1}^{n}T_{1}(x_{2};Z_{j})\right)$
$\displaystyle=\frac{1}{n}\cdot\frac{k}{n}\cdot n\cdot
cov\left(v_{\mathscr{D}}(x_{1})^{\top}\nabla_{\theta}m(\theta_{0},Z_{n}),T_{1}(x_{2};Z_{n})\right)$
$\displaystyle=\frac{k}{n}\cdot\mathds{E}_{\mathscr{D}}[v_{\mathscr{D}}(x_{1})^{\top}\nabla_{\theta}m(\theta_{0},Z_{n})\cdot
T_{1}(x_{2};Z_{n})],\;[\text{since }\mathds{E}_{\mathscr{D}}[T_{1}]=0]$

Decomposing the M-estimator directional derivative: Define the empirical
version of $v_{\mathscr{D}}$ to be
$v_{\hat{\mathscr{D}}}(x)=-\nabla_{\theta}^{\top}\eta(\hat{\theta},x)\left(\mathds{E}_{\hat{\mathscr{D}}}\left[\nabla_{\theta}^{2}m(\hat{\theta},Z)\right]\right)^{-1}$.
Then we decompose the directional derivative for $\eta(\hat{\theta},x_{1})$,
given by
$U_{i}(x_{1})=v_{\hat{\mathscr{D}}}(x_{1})^{\top}\nabla_{\theta}m(\hat{\theta},Z_{i})$
as $U_{i}(x_{1})=B_{i}(x_{1})+S_{i}(x_{1})$, where
$B_{i}(x)=v_{\mathscr{D}}(x)^{\top}\nabla_{\theta}m(\theta_{0},Z_{i})$.

Decomposing the random forest directional derivative: We know that for a
random forest the directional derivatives are given by
$U_{i}^{\prime}(x_{2})=n\cdot cov_{b}(N_{i,b},T_{b}(x_{2}))$ where the
covariance is over $b=1,\dots,B$ trees, $N_{i,b}$ is the number of times the
$i$th training datapoint is in the $b$th tree and $T_{b}$ is the $b$th tree
kernel. Now

$\displaystyle U_{i}^{\prime}(x_{2})$ $\displaystyle=n\cdot
cov_{b}(N_{i,b},T_{b}(x_{2}))$
$\displaystyle=n\left(\mathds{E}_{V\sim\hat{\mathscr{D}}}[\hat{F}(x_{2})\mid
V_{1}=Z_{i}]-\mathds{E}_{V\sim\hat{\mathscr{D}}}[\hat{F}(x_{2})]\right)$
$\displaystyle=k\left(\mathds{E}_{V\underset{\sim}{\subset}\hat{\mathscr{D}}}[T(x_{2})\mid
V_{1}=Z_{i}]-\mathds{E}_{V\underset{\sim}{\subset}\hat{\mathscr{D}}}[T(x_{2})]\right)$
$\displaystyle=k(A_{i}(x_{2})+R_{i}(x_{2})),\text{ where
}A_{i}(x)=\mathds{E}_{V\underset{\sim}{\subset}\hat{\mathscr{D}}}[\mathring{T}(x)\mid
V_{1}=Z_{i}]-\mathds{E}_{V\underset{\sim}{\subset}\hat{\mathscr{D}}}[\mathring{T}(x)]$

Here $V\underset{\sim}{\subset}\hat{\mathscr{D}}$ denotes subsampling (without
replacement) from the empirical distribution.

Finally we know the following:

$\displaystyle\mathds{E}\left[\frac{k}{n^{2}}\sum_{i=1}^{n}B_{i}(x_{1})A_{i}(x_{2})\right]-cov_{\mathscr{D}}(\eta(\hat{\theta},x_{1}),\mathring{F}(x_{2}))$
$\displaystyle\xrightarrow{\text{a.e.}}0\quad(\text{Lemma }\ref{lem:covcons})$
(C.1) $\displaystyle\frac{1}{n^{2}}\sum_{i=1}^{n}S_{i}(x_{1})^{2}$
$\displaystyle\xrightarrow{p}0\quad(\text{Lemma }\ref{lem:m_est_extra})$ (C.2)
$\displaystyle\frac{k^{2}}{n^{2}}\sum_{i=1}^{n}R_{i}(x_{2})^{2}$
$\displaystyle\xrightarrow{p}0\quad(\text{Lemma 13,
\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{wager2017estimation}{\@@citephrase{(}}{\@@citephrase{)}}}})$

Thus using the Cauchy-Schwartz inequality we see that

$\mathds{E}\left[\frac{k}{n^{2}}\sum_{i=1}^{n}U_{i}(x_{1})U^{\prime}_{i}(x_{2})\right]-cov_{\mathscr{D}}(\eta(\hat{\theta},x_{1}),\mathring{F}(x_{2}))\xrightarrow{\text{a.e.}}0$

thus proving that the IJ covariance estimate is consistent

######  Lemma 1.

Equation (C.1) holds

###### Proof.

We know from wager2017estimation that

$A_{i}(x)=\frac{n-k}{n}\left[T_{1}(x;Z_{i})-\frac{1}{n-1}\sum_{j\neq
i}T_{1}(x;Z_{j})\right]$

Thus

$\displaystyle\mathds{E}[B_{i}(x_{1})A_{i}(x_{2})]$
$\displaystyle=\frac{n-k}{n}\cdot\mathds{E}[v_{\mathscr{D}}(x_{1})^{\top}\nabla_{\theta}m(\theta_{0},Z_{i})\cdot
T_{1}(x_{2};Z_{i})]$
$\displaystyle\mathds{E}\left[\frac{k}{n^{2}}\sum_{i=1}^{n}B_{i}(x_{1})A_{i}(x_{2})\right]$
$\displaystyle=\frac{k}{n^{2}}\cdot
n\cdot\frac{n-k}{n}\cdot\mathds{E}[v_{\mathscr{D}}(x_{1})^{\top}\nabla_{\theta}m(\theta_{0},Z_{n})\cdot
T_{1}(x_{2};Z_{n})]$
$\displaystyle\implies\mathds{E}\left[\frac{k}{n^{2}}\sum_{i=1}^{n}B_{i}(x_{1})A_{i}(x_{2})\right]$
$\displaystyle-
cov_{\mathscr{D}}(\eta(\hat{\theta},x_{1}),\mathring{F}(x_{2}))\xrightarrow{\text{a.e.}}0\qquad\left[\because
1-\frac{k}{n}\to 1\right]$

∎

######  Lemma 2.

Equation (C.2) holds

###### Proof.

We use $v_{\hat{\mathscr{D}}}$ and $v_{\mathscr{D}}$ without its argument
$x_{1}$ to reduce notational clutter. Note that

$\displaystyle
S_{i}(x_{1})=v_{\hat{\mathscr{D}}}^{\top}\nabla_{\theta}m(\hat{\theta},Z_{i})-v_{\mathscr{D}}^{\top}\nabla_{\theta}m(\theta_{0},Z_{i})$
$\displaystyle\implies\frac{1}{n^{2}}\sum_{i=1}^{n}S_{i}(x_{1})^{2}$
$\displaystyle=\frac{1}{n^{2}}\sum_{i=1}^{n}\left(v_{\hat{\mathscr{D}}}^{\top}\nabla_{\theta}m(\hat{\theta},Z_{i})-v_{\mathscr{D}}^{\top}\nabla_{\theta}m(\theta_{0},Z_{i})\right)^{2}$
$\displaystyle=\frac{1}{n^{2}}\sum_{i=1}^{n}\left((v_{\hat{\mathscr{D}}}-v_{\mathscr{D}})^{\top}\nabla_{\theta}m(\hat{\theta},Z_{i})-v_{\mathscr{D}}^{\top}\left(\nabla_{\theta}m(\hat{\theta},Z_{i})-\nabla_{\theta}m(\theta_{0},Z_{i})\right)\right)^{2}$
$\displaystyle=\frac{1}{n^{2}}\sum_{i=1}^{n}2\left(\left((v_{\hat{\mathscr{D}}}-v_{\mathscr{D}})^{\top}\nabla_{\theta}m(\hat{\theta},Z_{i})\right)^{2}+\left(v_{\mathscr{D}}^{\top}(\nabla_{\theta}m(\hat{\theta},Z_{i})-\nabla_{\theta}m(\theta_{0},Z_{i}))\right)^{2}\right)$
$\displaystyle\leq\frac{2}{n^{2}}\sum_{i=1}^{n}\|v_{\hat{\mathscr{D}}}-v_{\mathscr{D}}\|^{2}\|\nabla_{\theta}m(\hat{\theta},Z_{i})\|^{2}+\frac{2}{n^{2}}\sum_{i=1}^{n}\|v_{\mathscr{D}}\|^{2}\|\nabla_{\theta}m(\hat{\theta},Z_{i})-\nabla_{\theta}m(\theta_{0},Z_{i})\|^{2}$
$\displaystyle=\frac{2\|v_{\hat{\mathscr{D}}}-v_{\mathscr{D}}\|^{2}}{n}\cdot\mathds{E}_{\hat{\mathscr{D}}}\left[\|\nabla_{\theta}m(\hat{\theta},Z)\|^{2}\right]+\frac{2\|v_{\mathscr{D}}\|^{2}}{n^{2}}\sum_{i=1}^{n}\|\nabla_{\theta}^{2}m(\tilde{\theta}_{i},Z_{i})(\hat{\theta}-\theta_{0})\|^{2}$
$\displaystyle\qquad[\text{ where
}\|\tilde{\theta}_{i}-\theta_{0}\|\leq\|\hat{\theta}-\theta_{0}\|]$
$\displaystyle\leq\frac{2\|v_{\hat{\mathscr{D}}}-v_{\mathscr{D}}\|^{2}}{n}\cdot\mathds{E}_{\hat{\mathscr{D}}}\left[\|\nabla_{\theta}m(\hat{\theta},Z)\|^{2}\right]+\frac{2\|v_{\mathscr{D}}\|^{2}}{n^{2}}\sum_{i=1}^{n}\|\nabla_{\theta}^{2}m(\tilde{\theta}_{i},Z_{i})\|_{1}^{2}\cdot\|\hat{\theta}-\theta_{0}\|^{2}$
$\displaystyle=\frac{2\|v_{\hat{\mathscr{D}}}-v_{\mathscr{D}}\|^{2}}{n}\cdot\mathds{E}_{\hat{\mathscr{D}}}\left[\|\nabla_{\theta}m(\hat{\theta},Z)\|^{2}\right]+\frac{2\|v_{\mathscr{D}}\|^{2}\cdot\|\hat{\theta}-\theta_{0}\|^{2}}{n}\cdot\mathds{E}_{\hat{\mathscr{D}}}\left[\|\nabla_{\theta}^{2}m(\tilde{\theta},Z)\|_{1}^{2}\right]$

Now it is easily seen that

$\displaystyle\|v_{\hat{\mathscr{D}}}-v_{\mathscr{D}}\|$
$\displaystyle\xrightarrow{p}0$ $\displaystyle\|\hat{\theta}-\theta_{0}\|^{2}$
$\displaystyle\xrightarrow{p}0$
$\displaystyle\mathds{E}_{\hat{\mathscr{D}}}\left[\|\nabla_{\theta}m(\hat{\theta},Z)\|^{2}\right]$
$\displaystyle\xrightarrow{a.e.}\mathds{E}_{\mathscr{D}}\left[\|\nabla_{\theta}m(\theta_{0},Z)\|^{2}\right]$
$\displaystyle\mathds{E}_{\hat{\mathscr{D}}}\left[\|\nabla_{\theta}^{2}m(\tilde{\theta},Z)\|_{1}^{2}\right]$
$\displaystyle\xrightarrow{a.e.}\mathds{E}_{\mathscr{D}}\left[\|\nabla_{\theta}^{2}m(\theta_{0},Z)\|_{1}^{2}\right]$

where the last two limits are due to continuity and the uniform strong law of
large numbers. Putting them together we see that
$\frac{1}{n^{2}}\sum_{i=1}^{n}S_{i}(x_{1})^{2}\xrightarrow{p}0$. ∎

##  Appendix D Local linear bias corrections for non-Gaussian responses

The local linear bias correction proposed in lu2021unified and explored in
section 4.3 will not be suitable to use for responses which are not
continuous. But if the response is from a general exponential family then we
could use techniques from ghosal2021generalised to define a local modification
- specifically the idea of looking at squared error from the perspective of
log-likelihoods. First note that the $\widehat{Bias}(x)$ defined in section
4.3 can be rewritten as

$\displaystyle\widehat{Bias}(x)$
$\displaystyle=\operatorname*{arg\,max}_{b}\sum_{k=1}^{n}w_{k}(x)\left[-\left(Y_{k}-(\hat{f}(X_{k})+b)\right)^{2}\right],$
$\displaystyle w_{k}(x)$
$\displaystyle=\sum_{b=1}^{B}\mathds{1}\left\\{Z_{k}\notin I_{b},X_{k}\in
L_{b}(x)\right\\}$

Here $w_{k}(x)$, the out-of-bag weight, is the numerator of $v_{k}(x)$ as
defined in section 4.3, i.e.,
$v_{k}(x)=\frac{w_{k}(x)}{\sum_{\ell=1}^{n}w_{\ell}(x)}$. For general
responses the goal would be to maximise a similar quantity involving the log-
likelihood - suppose we’ve fitted an initial estimator (an MLE-type estimate
or a GLM) $\hat{f}_{1}$ and then we’ve fitted a random forest $\hat{f}_{2}$
based on generalised residuals as defined in ghosal2021generalised. Then for a
query point $x$ we define the local bias correction to be

$\widehat{Bias}(x)=\operatorname*{arg\,max}_{b}\sum_{k=1}^{n}w_{k}(x)\ell\left(\hat{f}_{1}(X_{k})+\hat{f}_{2}(X_{k})+b;Y_{k},X_{k}\right),$

where $(Y_{k},X_{k})$ is the $k$th training data and $\ell$ is the log-
likelihood function with arguments in the link-space. Then the final link-
space prediction for the query point $x$ will be
$\hat{f}_{1}(x)+\hat{f}_{2}(x)+\widehat{Bias}(x)$. This optimisation step may
not always have a closed form solution and could be computationally expensive
to implement for each query point separately (compared to the case for
continuous responses). Also, theoretically the local bias will behave as an
M-estimator so we could use the directional derivatives derived above to
quantify uncertainty of this bias - hence enabling construction of confidence
intervals and comparison tests with the unmodified predictor
$\hat{f}_{1}+\hat{f}_{2}$.