# Leveraging Computer Vision in the Intensive Care Unit (ICU) for Examining
Visitation and Mobility

Scott Siegel  Jiaqing Zhang  Sabyasachi Bandyopadhyay  Subhash Nerella
Brandon Silva  Tezcan Baslanti  Azra Bihorac  Parisa Rashidi

###### Abstract

Despite the importance of closely monitoring patients in the Intensive Care
Unit (ICU), many aspects are still assessed in a limited manner due to the
time constraints imposed on healthcare providers. For example, although
excessive visitations during rest hours can potentially exacerbate the risk of
circadian rhythm disruption and delirium, it is not captured in the ICU.
Likewise, while mobility can be an important indicator of recovery or
deterioration in ICU patients, it is only captured sporadically or not
captured at all. In the past few years, the computer vision field has found
application in many domains by reducing the human burden. Using computer
vision systems in the ICU can also potentially enable non-existing assessments
or enhance the frequency and accuracy of existing assessments while reducing
the staff workload. In this study, we leverage a state-of-the-art noninvasive
computer vision system based on depth imaging to characterize ICU visitations
and patients’ mobility. We then examine the relationship between visitation
and several patient outcomes, such as pain, acuity, and delirium. We found an
association between deteriorating patient acuity and the incidence of delirium
with increased visitations. In contrast, self-reported pain, reported using
the Defense and Veteran Pain Rating Scale (DVPRS), was correlated with
decreased visitations. Our findings highlight the feasibility and potential of
using noninvasive autonomous systems to monitor ICU patients.

##  1 Introduction

The Intensive Care Unit (ICU) can be a traumatic experience for many patients
that adversely affects their long-term physical and cognitive outcomes [1, 2].
Pain [3, 4], immobility [5, 6], and sleep deprivation [7, 8] are common
sources of patient distress in the ICU [9]. Prolonged exposure to these causes
can induce the onset of delirium [7, 8], post-ICU post-traumatic stress
disorder (PTSD) [10, 11, 11], and increase the mortality risk [12, 13].

Monitoring and assessing many aspects of patient care in the ICU is generally
performed by nurses or physicians [3], which can be an error-prone practice
[14, 15, 16]. Nurses are often required to work long hours and provide
required shift coverage to satisfy patient demand [17]. These demands often
lead to increased stress, sleep deprivation [18, 19], and instances of
provider burnout [20, 21, 22]. Providing care to patients while fatigued can
lead to mistakes that negatively affect patient outcomes [23, 18, 19].
Checking on the patient can also adversely affect the patient’s sleep and lead
to sleep deprivation [24, 25, 26]. Patient-care interventions are frequently
mentioned as a cause of poor sleep in the ICU [27]. Some studies have found
that patients rarely go more than 90 minutes without receiving some form of
intervention [28, 27]. Introducing AI-augmented approaches to monitor patients
would reduce providers’ workload and incidents of burnout.

Recent advances in AI and sensing technologies have catalyzed the development
of AI-augmented solutions to monitor patient mobility and sources of sleep
disruption in the critical care setting [29]. Ma et al. employed three
Microsoft Kinect sensors and the Fast Region-based Convolutional Network (Fast
R-CNN) to measure patient activities [15, 30]. Yeung et al. used YOLOv2 to
detect and measure the duration of transitional activities, such as getting in
or out of bed [31, 32]. Davoudi et al. used a collection of sensors and deep
neural networks to monitor sources of sleep disruption as well as characterize
the physical activity and facial expression patterns exhibited by delirious
patients [33]. Using sound sensors, Davoudi et al. found that average sound
pressure levels from the delirious patients’ rooms were roughly equivalent to
street traffic noise. This finding is concerning because patients have
reported sound as being one of the major etiologies of poor sleep quality
[24].

In this study, we examined how patient mobility and personnel-related
disturbances correlate with measures of patients’ health, i.e., self-reported
pain using DVPRS [34, 35], patient acuity [36, 37], and delirium [13]. In an
IRB-approved study at the University of Florida (UF), we recruited 44 adult
patients admitted to the surgical, medical, trauma, vascular, and neuromedcine
ICUs. To quantify patient mobility and sleep disturbances, we employed an
Azure Microsoft Kinect sensor [38] in the ICU to record privacy-preserving
depth video [39] of patients during their stay. We developed and deployed a
full-stack annotation tool that enabled human annotators to label the presence
and posture of individuals in hospital rooms. We used the labeled data to
train the YOLOv8 [40, 41] neural network to detect individuals in each depth
frame and classify the activity as lying in bed or standing. After training
the model on selected annotations, we ran the best-performing model on depth
video frame sequences that were temporally aligned with timestamps of clinical
outcomes of interest. Within each sequence, we extracted 1) the proportion of
lying in bed frames, 2) the average number of people in the room, and 3) the
variance of people in the room. These three metrics were examined against our
clinical outcomes of interest, i.e., self-reported pain scores, patient acuity
[36, 37], and delirium. Our findings revealed significant correlations between
indicators of patient outcomes. Patients reported being in less pain when
hospital rooms had higher visitations. In contrast, objective measurements,
including patient acuity and instances of delirium, were significantly less
stable when more individuals were in the room.

##  2 Methods

###  2.1 Study participants

The data collected in this study were obtained from adult patients admitted to
the surgical, medical, trauma, vascular, and neuromedcine ICUs at the
University of Florida Shands Health Hospital. This study was approved by the
University of Florida Institutional Review Board by IRB201400546,
IRB202101013, and IRB201900354. Informed consent was obtained from all
participants, and all procedures complied with the Declaration of Helsinki and
university guidelines and regulations.

The cohort used in this study includes 44 patients. Demographic information
for the 44 patients is shown in Table 1.

Table 1: Patients Characteristics Characteristics | Mean
---|---
| N | %
Age | 56.9 | 15.6
Gender |  |
Male | 27 | 61%
Female | 17 | 39%
Race |  |
White | 36 | 82%
Black | 5 | 11%
Asian | 2 | 5%
Other | 1 | 2%
Hospital stay |  |
(hours) | 93.6 |
Outcomes |  |
Pain score | 2.4 |
Acuity | 0.69 |
Delirium | 0.24 |

###  2.2 Data Collection

We captured depth frames in the ICU using Microsoft Kinect cameras [38], which
captured 3D volumetric data of the ICU environment while safeguarding the
physical identities of each person. Six compact capture and compute systems
were used for data collection, each equipped with a depth camera and a compact
computer. The cameras had a field of view (FOV) of $70^{\circ}$ x $60^{\circ}$
and produced 640 x 576 px image frames. Frames were saved at a frequency of 1
FPS. Clinical coordinators oversaw the deployment and use of these systems to
ensure patient care was not impeded.

After the data was captured by the local capture and compute system, it was
encrypted, compressed, and transferred through a secure link to a high-
performance server for data annotation and model training.

###  2.3 Data Annotation

To build our mobility detection dataset, we decided to develop an in-house
bounding box annotation tool that would enable remote users to complete
annotations without installing any software. We designed our depth annotation
tool as a full-stack web application that enables annotators to draw bounding
boxes around every person portrayed in an image and label each box as one of
the following postures: lying in bed, standing, sitting on a chair, and
sitting on a bed. The front-end user interface was developed using the React
JavaScript framework (Figure 1). The front end was connected to a back-end
MongoDB database using the FastAPI Python framework. Before starting the
annotation, annotators completed the necessary steps to acquire Institutional
Review Board (IRB) approval. They were provided with a Standard Operating
Procedures (SOP) document that comprehensively outlined the correct annotation
procedures for each frame, and they were assigned an initial training dataset.
An additional annotation reviewer evaluated the quality of the annotations.

Figure 1: The front-end user interface to annotate depth frames. Posture
classes include lying in bed, standing, sitting on a bed, and sitting on a
chair. The tool also includes an assisted mobility class, which is intended to
capture instances when someone is physically assisting the patient with their
mobility.

###  2.4 YOLOv8 Architecture

We used the YOLOv8 [42] architecture for our patient activity detection model.
The YOLOv8’s backbone is the Darknet-53 convolutional neural network with an
added C2f module [42, 43]. The C2f module, which stands for ”cross-stage
partial bottleneck with two convolutions,” improves detection accuracy by
combining high-level features with contextual information. To improve YOLOv8’s
ability to detect smaller objects, bounding box loss is computed using the
complete intersection over union (CIoU) [44] and distribution focal loss (DFL)
[45] loss functions. The YOLOv8 architecture also integrates feature pyramid
networks, which are adept at recognizing objects at different scales. This is
beneficial for our depth frame detection task, given the varying distances of
patients within the frames.

###  2.5 Training Pipeline

The annotated dataset in this study included 30,840 annotated depth frames
from 44 patients. Model performance was evaluated using 5-fold cross-
validation. Patient IDs were used to split training and test sets within each
fold to prevent data leakage [46]. In each fold, 20% of the patients were
allocated for testing, while the remaining patients were used for training.
Experiments were conducted on a single A100 GPU. The default hyperparameters
provided by Ultralytics [41] were used for our experiments. The model was
trained for 100 epochs using stochastic gradient descent (SGD) and a learning
rate of 0.01. We report the mean average precision (mAP), precision, recall,
and F1-Score using an intersection over union threshold (IoU) of 0.5 as an
average across five folds.

###  2.6 Characterizing mobility and sleep disturbances

The total number of detected lying in bed and standing labels within each
frame was used to measure patient visitations average and variance. Walking in
and out of rooms and checking different equipment can produce noise and
disturb patients’ sleep. We evaluated sleep disturbances during both day and
night as it has been shown that approximately 50% of ICU patients’ sleep
occurs during the day [24, 27, 47].

###  2.7 Patient Outcomes

We assessed patient outcomes by evaluating three distinct scores, namely self-
reported pain, patient acuity [36, 37], and occurrences of delirium. The
patients reported pain using the Defense and Veterans Pain Rating Scale
(DVPRS) [35], which is a numerical scale ranging from 0 to 10, where 0
indicated the absence of pain, 1-4 signified mild pain, 5-6 indicated moderate
pain, and 7-10 represented severe pain. Patient acuity was evaluated using the
methodology proposed by a recent study [36, 37]. Patient acuity was classified
into two states: stable and unstable. Patients were classified as unstable if,
during the specified time interval, they received continuous renal replacement
therapy (CRRT), required ventilation, received a vasopressor infusion, or
required a blood transfusion of 10 units or more in the past 24 hours.
Delirium was evaluated using the Richmond Agitation Sedation Scale (RASS), the
Confusion Assessment Method for the Intensive Care Unit (CAM-ICU), and the
Glasgow Coma Scale (GCS) [48, 49].

###  2.8 Associating model metrics with patient Outcomes

The timestamps corresponding to clinical events of interest were used to
select temporally aligned depth video frames. Patient acuity was recorded in
four-hour intervals. For pain record timestamps, video frames within the 15-30
minute period prior to the pain timestamp were selected. We excluded the
15-minute window preceding each timestamp as it was likely to encompass the
frames captured when the caregiver was in the room. Including these frames has
the potential to introduce bias. For delirium record timestamps, video frames
were chosen if they occurred in the 15-60 minute window preceding the recorded
time.

Model metrics were recorded for each respective sequence of video frames, then
evaluated for normality using D’Agostino and Pearson’s test [50, 51]. Because
the null hypothesis that each sample came from a normal distribution could be
rejected, the Mann-Whitney U rank test [52] was used to evaluate independence
between samples. Samples were also split based on whether the record
timestamps occurred during the day (7 AM - 7 PM) or night (7 PM - 7 AM). (need
to mention type 1 error correction)

##  3 Results

###  3.1 YOLOv8 experiments

In this section, we report the performance of the YOLOv8 models that were
trained on our mobility detection dataset and evaluated using 5-fold cross-
validation. In Table 2, we present the experimental results of our posture
detection models across five independent test sets and report average mAP,
Precision, Recall, and F1 Score. The results are stratified by the predicted
posture. The model obtained a mAP of 0.88 for the standing class and a mAP of
0.79 for the lying-in-bed class.

Table 2: Experimental results on test dataset | Performance
---|---
Class | Mean (STD)
| mAP | Precision | Recall | F1 Score
Lying in bed | 0.79 (0.21) | 0.79 (0.20) | 0.74 (0.18) | 0.77 (0.19)
Standing | 0.88 (0.17) | 0.90 (0.03) | 0.83 (0.21) | 0.84 (0.14)
All classes | 0.83 (0.17) | 0.84 (0.09) | 0.79 (0.18) | 0.81 (0.14)

###  3.2 Associations with Patient Outcomes

For ICU room visitations, our study uncovered noteworthy positive correlations
with subjective measures of patient outcomes, such as self-reported pain,
alongside negative correlations with objective measures like delirium and
acuity. This relationship persisted even when we divided the samples into
night and day occurrences, revealing consistent patterns across all periods.
The following sections explore each patient metric in more detail.

####  3.2.1 Self-Reported Pain

Patients consistently reported lower pain scores when more people were in the
room (Table 3). In table 3, pain groups were combined based on whether the
pain scores were classified into the no pain and mild pain groups versus the
moderate and severe pain groups. The Kruskal-Wallis H-test [53] was also used
to compare the four pain groups separately, but no significant difference was
found ($p=0.9$). The proportion of frames with the patient lying in bed was
higher for the no and mild pain group compared to the moderate and severe pain
group, but the difference wasn’t statistically significant.

Splitting samples based on whether they occurred during the day versus night
yielded similar trends as those found with the combined datasets (Table 3).
The difference between night and day samples for the average number of people
detected is shown in Figure 2.

Table 3: Pain Associations with Model Metrics Model Metrics | No/Mild pain | Moderate/High Pain | p-value
---|---|---|---
N | Mean (STD) | Max (Min) | N | Mean (STD) | Max (Min)
Visitation Average |  |  |  |  |  |  |
Day | 973 | 1.77 (0.65) | 4.32 (1.0) | 307 | 1.62 (0.58) | 3.45 (1.0) | $3.06e^{-04}$
Night | 699 | 1.45 (0.52) | 4.31 (1.0) | 320 | 1.48 (0.51) | 3.22 (1.0) | 0.219
Combined | 1672 | 1.64 (0.62) | 4.32 (1.0) | 627 | 1.55 (0.55) | 3.45 (1.0) | 0.008
Visitation Variance |  |  |  |  |  |  |
Day | 973 | 0.41 (0.4) | 3.0 (0.0) | 307 | 0.34 (0.34) | 2.32 (0.0) | 0.002
Night | 699 | 0.24 (0.28) | 1.88 (0.0) | 320 | 0.26 (0.31) | 2.27 (0.0) | 0.402
Combined | 1672 | 0.34 (0.36) | 3.0 (0.0) | 627 | 0.3 (0.33) | 2.32 (0.0) | 0.01
Lying in Bed Proportion |  |  |  |  |  |  |
Day | 973 | 0.89 (0.26) | 1.0 (0.0) | 307 | 0.87 (0.29) | 1.0 (0.0) | 0.305
Night | 699 | 0.95 (0.19) | 1.0 (0.0) | 320 | 0.91 (0.26) | 1.0 (0.0) | 0.093
Combined | 1672 | 0.91 (0.23) | 1.0 (0.0) | 627 | 0.89 (0.27) | 1.0 (0.0) | 0.549
Figure 2: Relationship between pain levels and the average number of people in
the room during the night versus day.

####  3.2.2 Patient Acuity

For patient acuity, hospital room visitations average and variance were both
significantly higher for unstable frame sequences compared to stable sequences
(Table 4). The same pattern was found for the proportion of lying-in-bed
samples, but the difference between groups wasn’t statistically significant.
Similar to the patterns found for the pain associations, average day versus
night hospital room visitations and variance were greater for stable patients
than unstable patients (Table 4). The difference between night and day samples
for the average number of people detected is shown in Figure 3.

Table 4: Patient Acuity Associations with Model Metrics Model Metrics | Stable | Unstable | p-value
---|---|---|---
N | Mean (STD) | Max (Min) | N | Mean (STD) | Max (Min)
Visitation Average |  |  |  |  |  |  |
Day | 294 | 1.64 (0.38) | 3.33 (1.04) | 133 | 1.85 (0.5) | 3.48 (1.04) | $5.05e^{-05}$
Night | 285 | 1.41 (0.3) | 2.83 (1.01) | 127 | 1.5 (0.38) | 2.93 (1.05) | 0.034
Combined | 579 | 1.52 (0.36) | 3.33 (1.01) | 260 | 1.68 (0.48) | 3.48 (1.04) | $3.96e^{-05}$
Visitation Variance |  |  |  |  |  |  |
Day | 294 | 0.56 (0.31) | 1.95 (0.05) | 133 | 0.68 (0.38) | 2.48 (0.05) | 0.001
Night | 285 | 0.38 (0.27) | 1.98 (0.01) | 127 | 0.42 (0.27) | 1.77 (0.05) | 0.035
Combined | 579 | 0.47 (0.31) | 1.98 (0.01) | 260 | 0.56 (0.36) | 2.48 (0.05) | $2.51e^{-04}$
Lying in Bed Proportion |  |  |  |  |  |  |
Day | 294 | 0.85 (0.25) | 1.0 (0.0) | 133 | 0.93 (0.14) | 1.0 (0.11) | 0.082
Night | 285 | 0.93 (0.2) | 1.0 (0.0) | 127 | 0.96 (0.13) | 1.0 (0.01) | 0.64
Combined | 579 | 0.89 (0.23) | 1.0 (0.0) | 260 | 0.95 (0.13) | 1.0 (0.01) | 0.161
Figure 3: Relationship between patient acuity and the average number of people
in the room during the night versus day.

####  3.2.3 Delirium

Associations with delirium exhibited similar trends to those found for patient
acuity. Hospital room visitation average and variance were significantly
higher for delirium frame sequences than non-delirious sequences (Table 5).
The same trend was found for the proportion of lying-in-bed samples, but the
difference between groups wasn’t statistically significant. Differences
between day versus night hospital room visitation average and variance were
more significant for delirium samples than for non-delirious samples.
Comparing day versus night differences for delirium samples isn’t as reliable
because the sample sizes were smaller compared to those for pain and acuity.

Table 5: Delirium Associations with Model Metrics Model Metrics | Non-Delirious | Delirium | p-value
---|---|---|---
N | Mean (STD) | Max (Min) | N | Mean (STD) | Max (Min)
Visitation Average |  |  |  |  |  |  |
Day | 79 | 1.51 (0.42) | 2.51 (1.0) | 26 | 2.06 (0.53) | 2.98 (1.01) | $7.99e^{-06}$
Night | 87 | 1.6 (0.59) | 4.43 (1.0) | 28 | 1.48 (0.51) | 2.83 (1.0) | 0.245
Combined | 166 | 1.56 (0.52) | 4.43 (1.0) | 54 | 1.76 (0.59) | 2.98 (1.0) | 0.038
Visitation Variance |  |  |  |  |  |  |
Day | 79 | 0.35 (0.31) | 1.44 (0.0) | 26 | 0.58 (0.41) | 2.01 (0.01) | 0.001
Night | 87 | 0.34 (0.28) | 1.43 (0.0) | 28 | 0.36 (0.4) | 1.34 (0.0) | 0.669
Combined | 166 | 0.34 (0.29) | 1.44 (0.0) | 54 | 0.47 (0.42) | 2.01 (0.0) | 0.047
Lying in Bed Proportion |  |  |  |  |  |  |
Day | 79 | 0.96 (0.13) | 1.0 (0.0) | 26 | 0.93 (0.16) | 1.0 (0.32) | 0.389
Night | 87 | 0.93 (0.22) | 1.0 (0.0) | 28 | 0.99 (0.03) | 1.0 (0.89) | 0.573
Combined | 166 | 0.94 (0.18) | 1.0 (0.0) | 54 | 0.96 (0.11) | 1.0 (0.32) | 0.846

##  4 Discussion

In this study, we employed the state-of-the-art YOLOv8 deep learning model
[41] to detect patient postures. The results demonstrated that the model had a
harder time detecting the ”lying in bed class compared to the standing class.
This result aligns with our expectations because individuals lying in bed were
mostly covered by thick covers that made it difficult to distinguish which
features belonged to the person.

The metrics extracted from the model were the proportion of lying in bed
frames and the average and variance of hospital room visitations. The
proportion of lying in bed frames was evaluated because we could assume that
if someone was lying in bed, that person was the patient. If every frame in a
sequence contained a person lying in bed, that would indicate that they were
immobile.

Among the metrics extracted from the model, hospital room visitations average
and variance exhibited substantial correlations with measures of patient
outcomes. These findings are logical because noise is a crucial factor
associated with sleep deprivation and the onset of delirium [54, 55], and each
individual represents a unique source of noise. This observation highlights
the importance of using noninvasive autonomous systems to monitor patients.
Individuals who tend to wake easily during sleep may be disturbed when
healthcare providers open doors or inspect equipment [24, 25, 26]. Minimizing
the need for healthcare providers to enter or exit hospital rooms would allow
patients to experience uninterrupted sleep and facilitate their recovery.
Furthermore, implementing a noninvasive system for tracking patient mobility
would allow care providers to remotely assess whether a patient is asleep and
delay room entry until the patient is awake.

Another intriguing discovery was the difference in the relationships
identified when comparing subjective and objective measures of patient
outcomes. Higher levels of hospital room visitations were linked to adverse
objective measures of patient outcomes, including patient acuity and delirium.
Conversely, self-reported pain exhibited a contrasting pattern. According to
Boring et al., patients tend to downplay their pain when in the presence of
friends or family members, as they wish to avoid causing undue worry [56].
Since we couldn’t identify anyone in each frame, it wasn’t possible to know
whether the people in the room were care providers or family members.
Incorporating each person’s relationship to the patient in video sequences
would make it possible to explore the impact of family presence on self-
reported pain records further.

While a negative correlation was observed between patient outcomes and the
number of people in each room, the current analysis does not provide adequate
evidence to establish a causal relationship. The patient might have pressed
the nurse call button [57] to request aid or their condition may have worsened
to the extent that emergency medical staff were called in by the patient’s
monitoring equipment. This would suggest that the number of individuals in the
room resulted from the patients’ compromised outcomes. Incorporating the
timing of each event that could potentially disturb a patient is needed to
evaluate whether there is a causal relationship between hospital room
visitations and patients’ outcomes.

Using the average and variance of the number of people in an ICU room only
considers one factor that can adversely affect patients’ sleep. Other stimuli,
such as light or other sources of noise, have also been shown to impact
patients’ sleeping habits [24, 25, 26]. Incorporating real-time data from
noise and light sensors would improve our system’s capacity to identify the
critical factors contributing to sleep deprivation and notify care providers
of potentially serious events, such as the onset of delirium.

Monitoring the average amount of time a patient was lying in bed indicates how
much they were resting, but it doesn’t reveal how much time they spent
sleeping. Combining depth video with other sensors, such as a heart rate
monitor, an RGB video camera, or a polysomnography device, would improve our
system’s ability to monitor sleeping patterns.

There are advantages and disadvantages associated with the utilization of
depth cameras within an ICU as opposed to RGB cameras. Depth cameras provide
3D information about a scene and preserve the privacy of each person in a
scene as long as the resolution is low. However, identifying which person is
the patient can be a challenging endeavor when employing a depth camera. In
the future, we will incorporate patient tracking using one of the object
tracking architectures that have been integrated into the YOLOv8 model [58,
59]. Using tracking information, we can follow the person lying in the
hospital bed and evaluate their posture whenever they move around the hospital
room. Knowing which person the patient is will enable better mobility
characterization throughout the patient’s recovery period in the ICU. We can
establish a comprehensive mobility index by amalgamating this data with
information from other modalities, such as accelerometer data. This index
would be valuable for preventing adverse events and enhancing patient
outcomes.

##  5 Conclusions

We developed a full-stack application that enabled annotators to label 30,840
depth frames across 44 patients recovering in an ICU. We used this dataset to
train the state-of-the-art YOLOv8 deep learning model to characterize patient
mobility and evaluate hospital room visitations.

Our trained model revealed significant relationships between indicators of
patient outcomes and hospital room visitations. Subjective assessments, such
as self-reported pain, were lower when rooms had higher visitations. In
contrast, objective measurements, including patient acuity and instances of
delirium, were significantly higher when more individuals were in the room.
These findings illustrate the following: (1) the number of people in a
hospital room can be a potential source of bias for self-reported pain scores;
(2) greater quantities of people in a hospital room may have adverse effects
on the patient, such as disturbing the patients’ sleep.

The dataset utilized in this report is still in development and will be
expanded to incorporate patient tracking, enabling granular characterization
of their mobility and the identification of instances where patients’ mobility
is being assisted by other care providers. We aim to integrate this system
into an intelligent suite of sensors backed by state-of-the-art AI
technologies, enabling real-time monitoring and analyses of patient outcomes
in the ICU. By further refining and expanding our dataset, we aim to enhance
the capabilities of our system to detect and respond to changes in patient
mobility and overall health. This integrated approach holds the potential to
revolutionize patient care, offering timely interventions, reducing human
error, and ultimately improving the quality of care provided in intensive care
units.

## References

  * [1] Kristin Halvorsen et al.  “Patients’ experiences of well-being when being cared for in the intensive care unit—An integrative review” Publisher: Wiley Online Library  In _Journal of clinical nursing_ 31.1-2, 2022, pp. 3–19
  * [2] Cecilia Cuzco et al.  “Patients’ experience while transitioning from the intensive care unit to a ward” Publisher: Wiley Online Library  In _Nursing in critical care_ 27.3, 2022, pp. 419–428
  * [3] Juliana Barr et al.  “Clinical practice guidelines for the management of pain, agitation, and delirium in adult patients in the intensive care unit” Publisher: LWW  In _Critical care medicine_ 41.1, 2013, pp. 263–306
  * [4] Mina F Nordness, Christina J Hayhurst and Pratik Pandharipande  “Current perspectives on the assessment and management of pain in the intensive care unit” Publisher: Taylor & Francis  In _Journal of Pain Research_ , 2021, pp. 1733–1744
  * [5] Joseph Adler and Daniel Malone  “Early mobilization in the intensive care unit: a systematic review” Publisher: LWW  In _Cardiopulmonary physical therapy journal_ 23.1, 2012, pp. 5–13
  * [6] Claire J Tipping et al.  “The effects of active mobilisation and rehabilitation in ICU on mortality and function: a systematic review” Publisher: Springer  In _Intensive care medicine_ 43, 2017, pp. 171–183
  * [7] Gerald L Weinhouse et al.  “Bench-to-bedside review: delirium in ICU patients-importance of sleep deprivation” Publisher: Springer  In _Critical care_ 13, 2009, pp. 1–8
  * [8] Yuliya Boyko, H Ording and Poul Jennum  “Sleep disturbances in critically ill patients in ICU: how much do we know?” Publisher: Wiley Online Library  In _Acta anaesthesiologica scandinavica_ 56.8, 2012, pp. 950–958
  * [9] John W Devlin et al.  “Clinical practice guidelines for the prevention and management of pain, agitation/sedation, delirium, immobility, and sleep disruption in adult patients in the ICU” Publisher: LWW  In _Critical care medicine_ 46.9, 2018, pp. e825–e873
  * [10] Dimitry S Davydow et al.  “Posttraumatic stress disorder in general intensive care unit survivors: a systematic review” Publisher: Elsevier  In _General hospital psychiatry_ 30.5, 2008, pp. 421–434
  * [11] Cássia Righy et al.  “Prevalence of post-traumatic stress disorder symptoms in adult critical care survivors: a systematic review and meta-analysis” Publisher: Springer  In _Critical care_ 23, 2019, pp. 1–13
  * [12] Margaret A Pisani et al.  “Days of delirium are associated with 1-year mortality in an older intensive care unit population” Publisher: American Thoracic Society  In _American journal of respiratory and critical care medicine_ 180.11, 2009, pp. 1092–1097
  * [13] Sébastien Ouimet, Brian P Kavanagh, Stewart B Gottfried and Yoanna Skrobik  “Incidence, risk factors and consequences of ICU delirium” Publisher: Springer  In _Intensive care medicine_ 33, 2007, pp. 66–73
  * [14] Jaime M Beecroft et al.  “Sleep monitoring in the intensive care unit: comparison of nurse assessment, actigraphy and polysomnography” Publisher: Springer  In _Intensive care medicine_ 34, 2008, pp. 2076–2083
  * [15] Andy J Ma et al.  “Measuring patient mobility in the ICU using a novel noninvasive sensor” Publisher: LWW  In _Critical care medicine_ 45.4, 2017, pp. 630–636
  * [16] Hsiang-Ling Wang and Yun-Fang Tsai  “Nurses’ knowledge and barriers regarding pain management in intensive care units” Publisher: Wiley Online Library  In _Journal of clinical nursing_ 19.21-22, 2010, pp. 3188–3196
  * [17] S Betsiou et al.  “Nursing errors in intensive care unit and their association with burnout, anxiety, insomnia and working environment: a cross-sectional study” Publisher: Hippokratio General Hospital of Thessaloniki  In _Hippokratia_ 26.3, 2022, pp. 110
  * [18] Amy L Weaver, Sonja E Stutzman, Charlene Supnet and DaiWai M Olson  “Sleep quality, but not quantity, is associated with self-perceived minor error rates among emergency department nurses” Publisher: Elsevier  In _International emergency nursing_ 25, 2016, pp. 48–52
  * [19] Mohamed Zaki Ramadan and Khalid Saad Al-Saleh  “The association of sleep deprivation on the occurrence of errors by nurses who work the night shift” Publisher: Editura Medicală Universitară  In _Current health sciences journal_ 40.2, 2014, pp. 97
  * [20] Laurent Papazian et al.  “High-level burnout in physicians and nurses working in adult ICUs: a systematic review and meta-analysis” Publisher: Springer  In _Intensive Care Medicine_ 49.4, 2023, pp. 387–400
  * [21] Jin Jun et al.  “Relationship between nurse burnout, patient and organizational outcomes: Systematic review” Publisher: Elsevier  In _International journal of nursing studies_ 119, 2021, pp. 103933
  * [22] Thomas W Reader, Brian H Cuthbertson and Johan Decruyenaere  “Burnout in the ICU: potential consequences for staff and patient well-being” Pages: 4–6 Publication Title: Intensive care medicine Volume: 34  Springer, 2008
  * [23] Linsey M Barker and Maury A Nussbaum  “The effects of fatigue on performance in simulated nursing work” Publisher: Taylor & Francis  In _Ergonomics_ 54.9, 2011, pp. 815–829
  * [24] Abbas Al Mutair, Abbas Shamsan, Samer Salih and Awad Al-Omari  “Sleep deprivation etiologies among patients in the intensive care unit: literature review” Publisher: LWW  In _Dimensions of Critical Care Nursing_ 39.4, 2020, pp. 203–210
  * [25] Shailesh Bihari et al.  “Factors affecting sleep quality of patients in intensive care unit” Publisher: American Academy of Sleep Medicine  In _Journal of Clinical Sleep Medicine_ 8.3, 2012, pp. 301–307
  * [26] Kirsten M Eliassen and Laila A Hopstock  “Sleep promotion in the intensive care unit—a survey of nurses’ interventions” Publisher: Elsevier  In _Intensive and Critical Care Nursing_ 27.3, 2011, pp. 138–142
  * [27] Jorunn Beck Edvardsen and Fredrik Hetmann  “Promoting sleep in the intensive care unit” Publisher: SAGE Publications Sage CA: Los Angeles, CA  In _SAGE Open Nursing_ 6, 2020, pp. 2377960820930209
  * [28] Marita Ritmala-Castren et al.  “Sleep and nursing care activities in an intensive care unit” Publisher: Wiley Online Library  In _Nursing & health sciences_ 17.3, 2015, pp. 354–361
  * [29] Behnood Gholami, Wassim M Haddad and James M Bailey  “AI in the ICU: In the intensive care unit, artificial intelligence can keep watch” Publisher: IEEE  In _IEEE Spectrum_ 55.10, 2018, pp. 31–35
  * [30] Ross Girshick  “Fast r-cnn”  In _Proceedings of the IEEE international conference on computer vision_ , 2015, pp. 1440–1448
  * [31] Serena Yeung et al.  “A computer vision system for deep learning-based detection of patient mobilization activities in the ICU” Publisher: Nature Publishing Group UK London  In _NPJ digital medicine_ 2.1, 2019, pp. 11
  * [32] Joseph Redmon and Ali Farhadi  “YOLO9000: better, faster, stronger”  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 7263–7271
  * [33] Anis Davoudi et al.  “Intelligent ICU for autonomous patient monitoring using pervasive sensing and deep learning” Publisher: Nature Publishing Group UK London  In _Scientific reports_ 9.1, 2019, pp. 8020
  * [34] Rosemary C Polomano et al.  “Psychometric testing of the Defense and Veterans Pain Rating Scale (DVPRS): A new pain scale for military population” Publisher: Oxford University Press  In _Pain Medicine_ 17.8, 2016, pp. 1505–1519
  * [35] Chester C Buckenmaier III et al.  “Preliminary validation of the Defense and Veterans Pain Rating Scale (DVPRS) in a military population”  In _Pain Medicine_ 14.1, 2013, pp. 110–123
  * [36] Benjamin Shickel et al.  “Deep multi-modal transfer learning for augmented patient acuity assessment in the intelligent ICU” Publisher: Frontiers Media SA  In _Frontiers in digital health_ 3, 2021, pp. 640685
  * [37] Yuanfeng Ren et al.  “Development of Computable Phenotype to Identify and Characterize Transitions in Acuity Status in Intensive Care Unit”  In _arXiv preprint arXiv:2005.05163_ , 2020
  * [38] Zhengyou Zhang  “Microsoft kinect sensor and its effect” Publisher: IEEE  In _IEEE multimedia_ 19.2, 2012, pp. 4–10
  * [39] Edward Chou et al.  “Privacy-preserving action recognition for smart hospitals using low-resolution depth images”  In _arXiv preprint arXiv:1811.09950_ , 2018
  * [40] Joseph Redmon, Santosh Divvala, Ross Girshick and Ali Farhadi  “You only look once: Unified, real-time object detection”  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 779–788
  * [41] Glenn Jocher, Ayush Chaurasia and Jing Qiu  “YOLO by Ultralytics” original-date: 2022-09-11T16:39:45Z, 2023  URL: https://github.com/ultralytics/ultralytics
  * [42] Juan Terven and Diana Cordova-Esparza  “A comprehensive review of YOLO: From YOLOv1 to YOLOv8 and beyond”  In _arXiv preprint arXiv:2304.00501_ , 2023
  * [43] Joseph Redmon and Ali Farhadi  “Yolov3: An incremental improvement”  In _arXiv preprint arXiv:1804.02767_ , 2018
  * [44] Zhaohui Zheng et al.  “Distance-IoU loss: Faster and better learning for bounding box regression” Issue: 07  In _Proceedings of the AAAI conference on artificial intelligence_ 34, 2020, pp. 12993–13000
  * [45] Xiang Li et al.  “Generalized focal loss: Learning qualified and distributed bounding boxes for dense object detection”  In _Advances in Neural Information Processing Systems_ 33, 2020, pp. 21002–21012
  * [46] Karim Lekadir et al.  “FUTURE-AI: guiding principles and consensus recommendations for trustworthy artificial intelligence in medical imaging”  In _arXiv preprint arXiv:2109.09658_ , 2021
  * [47] S José Closs  “Patients’ night-time pain, analgesic provision and sleep after surgery” Publisher: Elsevier  In _International Journal of Nursing Studies_ 29.4, 1992, pp. 381–392
  * [48] Michael E Reznik et al.  “Fluctuations of consciousness after stroke: associations with the confusion assessment method for the intensive care unit (CAM-ICU) and potential undetected delirium” Publisher: Elsevier  In _Journal of critical care_ 56, 2020, pp. 58–62
  * [49] Shobhit Jain and Lindsay M Iverson  “Glasgow coma scale”, 2018
  * [50] RALPH D’agostino and Egon S Pearson  “Tests for departure from normality. Empirical results for the distributions of b^2 and \sqrt{b}” Publisher: Oxford University Press  In _Biometrika_ 60.3, 1973, pp. 613–622
  * [51] Ralph B d’Agostino  “An omnibus test of normality for moderate and large size samples” Publisher: Oxford University Press  In _Biometrika_ 58.2, 1971, pp. 341–348
  * [52] Henry B Mann and Donald R Whitney  “On a test of whether one of two random variables is stochastically larger than the other” Publisher: JSTOR  In _The annals of mathematical statistics_ , 1947, pp. 50–60
  * [53] William H Kruskal and W Allen Wallis  “Use of ranks in one-criterion variance analysis” Publisher: Taylor & Francis  In _Journal of the American statistical Association_ 47.260, 1952, pp. 583–621
  * [54] GCECM Mistraletti et al.  “Sleep and delirium in the intensive care unit” Publisher: Turino, Societa italiana di anestesiologia.  In _Minerva anestesiologica_ 74.6, 2008, pp. 329–334
  * [55] Julie S Cupka et al.  “The effect of non-pharmacologic strategies on prevention or management of intensive care unit delirium: a systematic review” Publisher: Faculty of 1000 Ltd  In _F1000Research_ 9, 2020
  * [56] Brandon L Boring et al.  “How and why patient concerns influence pain reporting: a qualitative analysis of personal accounts and perceptions of others’ use of numerical pain scales” Publisher: Frontiers  In _Frontiers in psychology_ 12, 2021, pp. 663890
  * [57] Jose Galinato, Mary Montie, Lance Patak and Marita Titler  “Perspectives of nurses and patients on call light technology” Publisher: LWW  In _CIN: Computers, Informatics, Nursing_ 33.8, 2015, pp. 359–367
  * [58] Yifu Zhang et al.  “Bytetrack: Multi-object tracking by associating every detection box”  In _European Conference on Computer Vision_ Springer, 2022, pp. 1–21
  * [59] Nir Aharon, Roy Orfaig and Ben-Zion Bobrovsky  “BoT-SORT: Robust associations multi-pedestrian tracking”  In _arXiv preprint arXiv:2206.14651_ , 2022
  * [60] OpenAI  “ChatGPT”, 2021  URL: https://openai.com/research/chatgpt


# Asymptotic Normality for Multivariate Random Forest Estimators

Kevin Li
Massachusetts Institute of Technology
<EMAIL_ADDRESS>I would like to thank my advisors Alberto Abadie and Victor
Chernozhukov for reviewing multiple drafts of this paper. In addition, Sophie
(Liyang) Sun, Ben Deaner, participants in the seminar class 14.386 (Spring
2020), and participants of MIT’s Econometrics Lunch seminar provided very
useful feedback. I would also like to thank Professor Stefan Wager for helping
me understand mechanics of covariance estimates in random forest models.

(December 7, 2020)

###### Abstract

Regression trees and random forests are popular and effective non-parametric
estimators in practical applications. A recent paper by Athey and Wager shows
that the random forest estimate at any point is asymptotically Gaussian; in
this paper, we extend this result to the multivariate case and show that the
vector of estimates at multiple points is jointly normal. Specifically, the
covariance matrix of the limiting normal distribution is diagonal, so that the
estimates at any two points are independent in sufficiently deep trees.
Moreover, the off-diagonal term is bounded by quantities capturing how likely
two points belong to the same partition of the resulting tree. Our results
relies on certain a certain stability property when constructing splits, and
we give examples of splitting rules for which this assumption is and is not
satisfied. We test our proposed covariance bound and the associated coverage
rates of confidence intervals in numerical simulations.

##  1 Introduction

Tree estimators and random forests are nonparmetric estimators that enjoy
widespread popularity in applied data science, powering models in e-commerice,
finance and macroneconomics, and medicine. Random forests, introduced by
Breiman [1], is a bagging algorithm that aggregates estimates from a
collection of tree estimators. Since each tree has alow bias but high
variance, aggregation helps performance by balancing the bias-variance trade-
off. Other than good empirical performance, random forests also enjoy the
following advantages over other nonparametric methods: they scale naturally
with high dimensional data, as the optimal cut at each may be computed in
parallel by quantizing each covariate; categorical variables and missing data
can be easily corporated; and they are more interpretable since variable
importance can be explicitly characterized.

One fruitful application of random forests in economics is in estimating
heterogeneneous treatment effects, i.e., a function of the form
$f(x)=\operatorname{\mathbf{E}}(Y_{1}-Y_{0}\mid X=x)$. In order to conduct
inference in econometric applications (e.g., to test the null $H_{0}:f(x)=0$),
knowledge about the rate of convergence or asymptotic distribution of the
estimator $\hat{f}(x)$ is required. Moreover, functions of $f$ are often of
practical interest: for example, when we wish to study the difference in
treatment effects for two different subpopulations, the quantity of interest
is

$f(x)-f(\bar{x}),\quad\text{where $x$ and $\bar{x}$ describe the two
subpopulations}.$ (1)

More generally, we might also be interested in a weigthed treatment effect,
where a subpopulation $x$ is given an importance weight $\mu$. In this case,
the functional of $f$ is

$\int_{x\in\mathcal{X}}f(x)d\mu,\quad\text{where $\mu$ is not necessarily the
density of $x$},$ (2)

and the integral is taken over the feature space $\mathcal{X}$. Inference on
functions of $f$ requires not only the asymptotic distribution of the point
estimate $f(x)$, but also the correlation between estimates at distinct points
$f(x)$ and $f(\bar{x})$. This paper studies the correlation structure of a
class of random forests models whose asymptotic distributions were worked out
in [2]. We find sufficient conditions under which the asymptotic covariance of
estimates at distinct poitns vanish relative to the variance; moreover, we
provide finite sample heuristics based on our calculations. To the best of our
knowledge, this is the first set of results on the correlation of random
forest estimates.

The present paper builds on and extends the results in [2], which in turn
builds on related work [3] on the concentration properties of tree estimators
and random estimators in general. See also [4], which extends the random
forest model considered here to estimate moment conditions. Stability results
established in this paper have appeared in [5], who study the notions of
algorithmic stability for random forests and logistic regression and derive
generalization error guarantees. Also closely related to our paper are [6] and
[7], concerning finite sample Gaussian approximations of sums and U-statistics
in high dimensions. In this context, our paper provides a stepping stone
towards applying the theory of finite sample U-statistics to random forests,
where bounds on the correlation structure plays a central role.

The paper is structured as follows. In Section 2, we introduce the random
forest model and and state the assumptions required for our results; Section 3
contains our maion theoretical contributions; Section 4 builds on Sections 3
and discusses heuristics useful in finite sample settings; Section 5 evalutes
our heuristics with numerical experiments, and Section 6 concludes. All proofs
are found in the appendix.

##  2 Setup

In this paper, we study Gaussian approximations of random forest estimators.
Throughout, we assume that a random sample
$(Z_{i})_{i=1}^{n}=(X_{i},Y_{i})_{i=1}^{n}$ is given, where
$X_{i}\subseteq\mathcal{X}\in\mathbf{R}^{p}$ is a vector of regressors and
$Y_{i}\in\mathbf{R}$ the response.

A tree estimator $T(\xi;Z_{1},\dots,Z_{s})$ works by recursively partitioning
the space of covariates $\mathcal{X}$ into a collection of non-overlapping
hyper-rectangles. The estimator starts with the entire space $\mathcal{X}$
(assumed to be bounded), and uses the points $Z_{1}$, …, $Z_{s}$ to produce a
coordinate $j\in\\{1,\dots,p\\}$ and an index $t\in\mathbf{R}$. This divides
$\mathcal{X}$ into two hyperrectangles,

$\\{x\in\mathcal{X}:x_{j}<t\\}\quad\text{and}\quad\\{x\in\mathcal{X}:x_{j}>t\\}.$
(3)

In the next step, each of these two sets is partitioned into two
hyperrectangles each, using the points that land in each. The process
continues until a particular hyperrectangle satisfies a stopping criterion.
This process corresponds to a tree in the natural way, so we will refer to a
generic hyperrectangle a “node”, and the leaves of the tree “terminal nodes”,
where splitting ceases. Given the partition of $\mathcal{X}$ into terminal
nodes $N_{1},\dots,N_{q}$, the prediction is the average of the responses

$T(x;\xi,Z_{i_{1}},\dots,Z_{i_{s}})=\sum_{i=1}^{q}\mathbf{1}(x\in
N_{i})\cdot\frac{1}{|N_{i}|}\sum_{i:X_{i}\in N_{i}}Y_{i},$ (4)

where $|N_{i}|$ is defined to be the number of samples $X_{i}$ inside $N_{i}$.
Here, $\xi$ plays the role of randomization, as the splitting decisions allows
for randomization.

Given a base learner $T$, a random forest estimator at $x$ is defined to be
the order $s$ $U$-statistic with kernel $T$, i.e.,

$\operatorname{RF}(x;Z_{1},\dots,Z_{n})=\frac{1}{\binom{n}{s}}\sum_{i_{1},\dots,i_{s}}\operatorname{\mathbf{E}}_{\xi}T(x;\xi,Z_{i_{1}},\dots,Z_{i_{s}}),$
(5)

where the randomization device $\xi$ is marginalized over.

In this paper, we study the asymptotic distribution of the vector

$\operatorname{RF}(x_{1},\dots,x_{q};Z_{1},\dots,Z_{n})\coloneqq\begin{bmatrix}\operatorname{RF}(x_{1};Z_{1},\dots,Z_{n})\\\
\vdots\\\
\operatorname{RF}(x_{q};Z_{1},\dots,Z_{n})\end{bmatrix}\in\mathbf{R}^{q},$ (6)

where the random forest is employed to produce an estimate at $q$ points
$x_{1},\dots,x_{q}$.

Previous work by (Athey and Wager) has shown that under regularity conditions

$\frac{\operatorname{RF}(x,Z_{1},\dots,Z_{n})-m(x)}{\sigma_{n}(x)}\xRightarrow{\;\rm
dist\;}N(0,1)\qquad\text{where }m(x)=\operatorname{\mathbf{E}}(Y\mid X=x)$ (7)

and $\sigma_{n}^{2}(x)$ the variance of $\operatorname{RF}(x)$. In this paper,
we extend their result to cover multivariate predictions and establish joint
normality

$\Sigma^{-1/2}\\{\operatorname{RF}(x_{1},\dots,x_{q})-m(x_{1},\dots,x_{q})\\}\xRightarrow{\;\rm
dist\;}N(0,I_{q\times q})$ (8)

where $m(x_{1},\dots,x_{q})$ is the vector with components
$\operatorname{\mathbf{E}}(Y\mid X=x_{k})$. This is the main technical result
of the paper. In addition, we provide numerical simulations to gauge the
finite sample behavior of the multivariate random forest estimator.

###  2.1 Assumptions

Our results builds on [2], so that we will work with the same model of random
forests to be maintained. Specifically, we will assume the following
properties of the underlying tree algorithm.

######  Assumption (Honesty).

Conditional on $\\{X_{i}\\}$, the tree structure (i.e., the splitting
coordinates and splitting indices) are independent of the tree estimates.

There are several ways to satisfy this assumption. The first is to calculate
splits based only on $X_{i}$: for example, a splitting coordinate and index is
chosen to minimize the average squared distance to the center of each of the
two nodes as in clustering algorithms. The second way is to perform data
splitting: partition the data into two sets $\mathcal{I}_{1}$ and
$\mathcal{I}_{2}$. Observations in $\mathcal{I}_{1}$ and
$X_{i}\in\mathcal{I}_{2}$ may be freely used during the splitting process,
while $Y_{i}\in\mathcal{I}_{2}$ are used to determine terminal node values.
Finally, a third method assumes the existence of side information
$\\{W_{i}\\}$: splits are made in the domain of $\\{X_{i}\\}$ using
$\\{W_{i}\\}$ are surrogate responses, with $\\{Y_{i}\\}$ being used for
prediction at the terminal nodes.

For simplicity, we will assume that the first scheme is used, namely that the
splitting algorithm uses only $X_{i}$ and not $Y_{i}$ in placing the splits.
Our results remain valid for the third scheme.

######  Assumption (Randomized Cyclic Splits).

At each non-terminal node, the splitting algorithm decides to do data-
independent split with probability $\delta>0$, by flipping a “$\delta$-coin”
with distribution is independent of everything else. The first time the coin
lands heads, the first coordinate is chosen as the splitting coordinate; the
second time, the second coordinate is chosen, and so on, such that that on the
$J$-th time the coin lands heads, the $(p\text{ mod }J)+1$-th coordinate is
chosen.

This is modification of the assumption in [2], in which each node has a
probability $\delta$ of being chosen at each split; the latter is most easily
implemented by flipping a $p\delta$ coin, and selecting one of the $p$
coordinates uniformly at random to split when the coin lands heads. Our
assumption does away with the randomization in the second step.

######  Assumption (The Splitting Algorithm is $(\alpha,k)$-Regular).

There exists some $\alpha\in(0,1/2)$ such that whenever a split occurs in a
node with $m$ sample points, the two hyper-rectangles contains at least
$\alpha m$ many points each. Moreover, splitting ceases at a node only if the
node contains less than $2k-1$ points for some $k$. In particular, this
implies that $\operatorname{\mathbf{Var}}T$ has bounded entries, since the
number of points in each terminal node is bounded above.

As shown in [3], this implies that with all but exponentially small
probability, the splitting axis also shrinks by a factor between $\alpha$ and
$1-\alpha$.

The following assumption is specific for our model. The first specifies that
the _potential_ splits do not depend on the data; in practice, this is
“essentially” without loss of generality as the domain of $X$ is fixed. For
example, the assumption is satisfied if all splits occur at indices that lie
on a specified grid, e.g., points which can be represented by a 64-bit
floating point number. In keeping with the $(\alpha,k)$ regularity assumption
above, we also assume that no potential split may split the splitting axis by
more than a proportion $\alpha$.

######  Assumption (Predetermined Splits).

The potential splits at each node does not depend on the realizations of the
data. Furthermore, the number of potential splits at each node is finite, and
each potential split shrinks the length of the corresponding axis by at most
$\alpha$.

We will also require other assumptions regarding the “stability” of splits;
these assumptions will be introduced later. Finally, we follow [2] and place
the following distributional assumptions on the data generating process.

######  Assumption (Distributional Assumptions on $(X,Y)$).

The covariate $X$ is distributed on the $[0,1]^{p}$ with a density that is
bounded away from zero and infinity; without loss of generality, we may assuem
that $X$ is uniformly distributed on the $p$-dimensional hypercube.

Furthermore, the regression function $x\mapsto\operatorname{\mathbf{E}}(Y\mid
X=x)$, $x\mapsto\operatorname{\mathbf{E}}(Y^{2}\mid X=x)$, and
$x\mapsto\operatorname{\mathbf{E}}(Y^{3}\mid X=x)$ are uniformly Lipschitz
continuous.

##  3 Gaussianity of Multivariate $U$-Statistics

###  3.1 Hajek Projections

For $m\geq 1$ and a statistic $f(Z_{1},\dots,Z_{m})\in\mathbf{R}^{q}$, recall
that the Hajek projection of $f$ is defined to be

$\mathring{f}(Z_{1},\dots,Z_{m})=\sum_{i=1}^{m}\operatorname{\mathbf{E}}[f(Z_{1},\dots,Z_{m})\mid
Z_{i}]-(m-1)\operatorname{\mathbf{E}}f(Z_{1},\dots,Z_{m}).$ (9)

In particular, when $Z_{1},\dots,Z_{m}$ is an IID sequence and $f$ is
symmetric in its arguments, we have

$\mathring{f}(Z_{1},\dots,Z_{m})=\sum_{i=1}^{m}f_{1}(Z_{i})-(m-1)\operatorname{\mathbf{E}}f,$
(10)

where $f_{1}(z)=\operatorname{\mathbf{E}}[f\mid Z_{1}=z]$.

The previous calculations applied to our setting where $f$ is the random
forest estimator yields

$\mathring{\operatorname{RF}}(Z_{1},\dots,Z_{n})-\mu=\sum_{i=1}^{n}\operatorname{\mathbf{E}}(\operatorname{RF}-\mu\mid
Z_{i})=\frac{1}{\binom{n}{s}}\sum_{i=1}^{n}\operatorname{\mathbf{E}}\biggl{[}\sum_{i_{1},\dots,i_{s}}\operatorname{\mathbf{E}}_{\xi}T(\xi;Z_{i_{1}},\dots,Z_{i_{s}})-\mu\mid
Z_{i}\biggr{]}.$ (11)

Since our samples are independent,
$\operatorname{\mathbf{E}}(\operatorname{\mathbf{E}}_{\xi}T(\xi;Z_{i_{1}},\dots,Z_{i_{s}})\mid
Z_{i})=\mu$ whenever $i\notin\\{i_{1},\dots,i_{s}\\}$. As
$\\{i_{1},\dots,i_{s}\\}$ runs over the size-$s$ subsets of $\\{1,\dots,n\\}$,
there are exactly $\binom{n-1}{s-1}$ subsets which contain $i$. For each of
these subsets

$\operatorname{\mathbf{E}}(\operatorname{\mathbf{E}}_{\xi}T(\xi;Z_{i_{1}},\dots,Z_{i_{s}})-\mu\mid
Z_{i})\eqqcolon T_{1}(Z_{i})-\mu,$ (12)

where $T_{1}(z)=\operatorname{\mathbf{E}}_{\xi}T(z,Z_{2},\dots,Z_{s})$.
Therefore,

$\mathring{\operatorname{RF}}-\mu=\frac{1}{\binom{n}{s}}\sum_{i=1}^{n}\binom{n-1}{s-1}(T_{1}(Z_{i})-\mu)=\frac{s}{n}\sum_{i=1}^{n}(T_{1}(Z_{i})-\mu).$
(13)

###  3.2 Asymptotic Gaussianity via Hajek Projections

The standard technique to derive the asymptotic distribution of a
$U$-statistic is to establish a lower bound on the variance of its Hajek
statistic; this is the approach taken by (Athey and Wager) and we follow the
approach here. Let $V$ be the variance of $\mathring{\operatorname{RF}}$;
using our previous result, we could write

$V=\operatorname{\mathbf{Var}}\biggl{[}\frac{s}{n}\sum_{i=1}^{n}(T_{1}(Z_{i})-\mu)\biggr{]}=\frac{s^{2}}{n}\operatorname{\mathbf{Var}}(T_{1}(Z_{1}))=\frac{s}{n}\operatorname{\mathbf{Var}}\biggl{[}\sum_{i=1}^{s}T_{1}(Z_{i})\biggr{]}=\frac{s}{n}\operatorname{\mathbf{Var}}\mathring{T},$
(14)

where $\mathring{T}$ is the Hajek projection of the statistic
$(Z_{1},\dots,Z_{s})\mapsto\operatorname{\mathbf{E}}_{\xi}T(\xi;Z_{1},\dots,Z_{s})$.

Under standard regularity conditions for a multivariate triangular array
Central Limit Theorem,

$V^{-1/2}(\mathring{\operatorname{RF}}-\mu)\xRightarrow{\;\rm dist\;}N(0,I).$
(15)

To establish the joint normality for $\operatorname{RF}$, write

$V^{-1/2}(\operatorname{RF}-\mu)=V^{-1/2}(\operatorname{RF}-\mathring{\operatorname{RF}})+V^{-1/2}(\mathring{\operatorname{RF}}-\mu).$
(16)

so all that is required is to prove that
$V^{-1/2}(\operatorname{RF}-\mathring{\operatorname{RF}})\xrightarrow{\;\operatorname{\mathbf{P}}\;}0$.
We will show that the quantity converges in squared mean. Setting
$e=V^{-1/2}(\operatorname{RF}-\mathring{\operatorname{RF}})$,

$\begin{split}\operatorname{\mathbf{E}}(e^{\intercal}e)&=\operatorname{\mathbf{E}}(\operatorname{RF}-\mathring{\operatorname{RF}})^{\intercal}V^{-1}(\operatorname{RF}-\mathring{\operatorname{RF}})=\operatorname{\mathbf{E}}\operatorname{tr}V^{-1}(\operatorname{RF}-\mathring{\operatorname{RF}})(\operatorname{RF}-\mathring{\operatorname{RF}})^{\intercal}\\\
&=\operatorname{tr}V^{-1}\operatorname{\mathbf{E}}(\operatorname{RF}-\mathring{\operatorname{RF}})(\operatorname{RF}-\mathring{\operatorname{RF}})^{\intercal}=\operatorname{tr}V^{-1/2}\operatorname{\mathbf{Var}}(\operatorname{RF}-\mathring{\operatorname{RF}})V^{-1/2}.\end{split}$
(17)

By Proposition 1, we have the following Hoeffding decomposition with respect
to the matrix $V^{-1}$

$\operatorname{RF}-\mathring{\operatorname{RF}}=\frac{1}{\binom{n}{s}}\biggl{[}\sum_{i<j}\binom{n-2}{s-2}(T_{2}(Z_{i},Z_{j})-\mu)+\sum_{i<j<k}\binom{n-3}{s-3}(T_{3}(Z_{i},Z_{j},Z_{k})-\mu)+\cdots\biggr{]},$
(18)

where $T_{2}$, $T_{3}$, etc. are the second and third order projections of $T$
obeying the orthogonality conditions

$\operatorname{\mathbf{E}}[(T_{k}-\mu)^{\intercal}V^{-1}(T_{k^{\prime}}-\mu)]=0,\qquad\text{for
$k\neq k^{\prime}$.}$ (19)

In addition, being projections of $T$, the higher order projections also
satisfy

$\operatorname{\mathbf{E}}[(T_{k}-\mu)^{\intercal}V^{-1}(T_{k}-\mu)]\leq\operatorname{\mathbf{E}}[(T-\mu)^{\intercal}V^{-1}(T-\mu)].$
(20)

Therefore, (14) and (17) imply that

$\operatorname{\mathbf{E}}(e^{\intercal}e)\leq\frac{s}{n}\operatorname{tr}(\operatorname{\mathbf{Var}}\mathring{T}^{-1}\operatorname{\mathbf{Var}}T).$
(21)

Athey and Wager obtains a bound on the diagonal elements of
$\operatorname{\mathbf{Var}}\mathring{T}^{-1}$ and
$\operatorname{\mathbf{Var}}T$, namely

$\frac{(\operatorname{\mathbf{Var}}T)_{kk}}{(\operatorname{\mathbf{Var}}\mathring{T})_{kk}}\leq
c(\log s)^{p},\qquad\text{for each $k=1,\dots,q$},$ (22)

for a constant $c$. In the next section, we shall bound the off-diagonal
elements of $\operatorname{\mathbf{Var}}\mathring{T}$ and
$\operatorname{\mathbf{Var}}T$ in order to establish
$\operatorname{\mathbf{E}}(e^{\intercal}e)\to 0$.

######  Proposition 1 (Hoeffding Decomposition for Multivariate
$U$-statistics).

Fix a positive definite matrix $M$. Let
$f(x_{1},\dots,x_{n})\in\mathbf{R}^{p}$ be a vector-valued function that is
symmetric in its arguments and let $X_{1},\dots,X_{n}$ be a random sample such
that $f(X_{1},\dots,X_{n})$ has finite variance. Then there exists functions
$f_{1},f_{2},\dots,f_{n}$ such that

$f(X_{1},\dots,X_{n})=\operatorname{\mathbf{E}}(f)+\sum_{i=1}^{n}f_{1}(X_{1})+\sum_{i<j}f_{2}(X_{i},X_{j})+\dots+f_{n}(X_{1},\dots,X_{n})$
(23)

where $f_{k}$ is a function of $k$ arguments, such that

$\operatorname{\mathbf{E}}f_{k}(X_{1},\dots,X_{k})=0\quad\text{and}\quad\operatorname{\mathbf{E}}[f_{k}(X_{1},\dots,X_{k})^{\intercal}Mf_{\ell}(X_{1},\dots,X_{l})]=0.$
(24)

###  3.3 Covariance Calculations

The aim of this subsection is to prove the following asymptotic behavior of
the covariance matrix $\operatorname{\mathbf{Var}}\mathring{T}$ under suitable
stability conditions of the splitting algorithm

$(\operatorname{\mathbf{Var}}\mathring{T})_{k,l}=o(s^{-\epsilon})\text{ for
all $k\neq l$ and some $\epsilon>0$}.$ (25)

The next proposition shows that this is enough to establish
$\operatorname{\mathbf{E}}(e^{\intercal}e)\to 0$.

######  Proposition 2.

Suppose $\operatorname{\mathbf{Var}}T$ and
$\operatorname{\mathbf{Var}}\mathring{T}$ satisfy the conditions in (25). Then

$\frac{s}{n}\operatorname{tr}\operatorname{\mathbf{Var}}\mathring{T}^{-1}\operatorname{\mathbf{Var}}T\to
0.$ (26)

######  Remark.

Proposition 2 does not require the fact
$(\operatorname{\mathbf{Var}}T)_{k,l}\to 0$ for the off-diagonal entries of
$\operatorname{\mathbf{Var}}T$, only that they are bounded. Moreover, we only
require the much weaker bound
$(\operatorname{\mathbf{Var}}\mathring{T})_{k,l}=o(\frac{1}{\log^{p}s})$.
However, tighter bounds are still useful for applications: for one they allow
us to control the “approximation” error
$\operatorname{\mathbf{E}}(e^{\intercal}e)$; second, they are useful for
bounding the off diagonal terms in $V^{-1}(\operatorname{RF}-\mu)$.

The rest of this section will establish the covariance bound for
$\operatorname{\mathbf{Var}}\mathring{T}$. Recall that

$\mathring{T}-\mu=\sum_{i=1}^{s}\operatorname{\mathbf{E}}(T\mid
Z_{i})\quad\text{so
that}\quad\operatorname{\mathbf{Var}}\mathring{T}=s\operatorname{\mathbf{Var}}(\operatorname{\mathbf{E}}(T\mid
Z_{1}))\text{ due to independence}.$ (27)

By the orthogonality condition,

$\operatorname{\mathbf{Var}}\operatorname{\mathbf{E}}(T\mid
Z_{1})=\operatorname{\mathbf{Var}}[\operatorname{\mathbf{E}}(T\mid
Z_{1})-\operatorname{\mathbf{E}}(T\mid
X_{1})]+\operatorname{\mathbf{Var}}[\operatorname{\mathbf{E}}(T\mid X_{1})],$
(28)

and furthermore, honesty implies

$\operatorname{\mathbf{E}}(T_{k}\mid
Z_{1})-\operatorname{\mathbf{E}}(T_{k}\mid
X_{1})=\operatorname{\mathbf{E}}(I_{k}\mid
X_{1})(Y_{1}-\operatorname{\mathbf{E}}(Y_{1}\mid I_{k}=1,X_{1}))$ (29)

where $T_{k}$ is the $k$-th coordinate of $T$ (i.e., the estimate of the tree
at $x_{k}$) and $I_{k}$ is the indicator variable for the event that $X_{1}$
belongs to the terminal node containing $x_{k}$. In particular, the off-
diagonal entry at $(j,k)$ of
$\operatorname{\mathbf{Var}}[\operatorname{\mathbf{E}}(T\mid
Z_{1})-\operatorname{\mathbf{E}}(T\mid X_{1})]$ is equal to

$\operatorname{\mathbf{E}}[\operatorname{\mathbf{E}}(I_{k}\mid
X_{1})\operatorname{\mathbf{E}}(I_{j}\mid
X_{1})(Y_{1}-\operatorname{\mathbf{E}}(Y_{1}\mid
X_{1},I_{k}=1)(Y_{1}-\operatorname{\mathbf{E}}(Y_{1}\mid X_{1},I_{j}=1)]$ (30)

Following the same truncation argument as in (Athey and Wager, Theorem 5),
this quantity is bounded by

$\operatorname{\mathbf{E}}[\operatorname{\mathbf{E}}(I_{k}\mid
X_{1})\operatorname{\mathbf{E}}(I_{j}\mid X_{1})].$ (31)

Alternatively, this bound also holds in the case when $Y$ is bounded. (Note
that a direct application of the Cauchy-Schwarz inequality yields a weaker
bound $\sqrt{\operatorname{\mathbf{E}}[\operatorname{\mathbf{E}}(I_{k}\mid
X_{1})^{2}\operatorname{\mathbf{E}}(I_{j}\mid X_{1})^{2}]}$.) Recall that
$\operatorname{\mathbf{E}}(I_{k}\mid X_{1}=x)$ is the probability that $x$ and
$x_{k}$ belong to the same terminal node; likewise
$\operatorname{\mathbf{E}}(I_{j}\mid X_{1}=x)$ is the probability that $x$ and
$x_{j}$ belong to the same terminal node. The following proposition shows the
intuitive result that for sufficiently large $s$, at least one of these
probabilities is small.

######  Proposition 3.

Fix two distinct points $x$ and $\bar{x}\in[0,1]^{p}$, and let $M(x,\bar{x})$
be the probability that $x$ and $\bar{x}$ belongs to the same terminal node.
If $\delta>1/2$, then

$M(x,\bar{x})=o(s^{-(1+\epsilon)})\quad\text{for some $\epsilon>0$}.$ (32)

The same holds for the quantity $M(x,\bar{x}\mid X_{1}=\bar{x})$, the
probability conditional on $X_{1}=\bar{x}$.

This proposition shows that the contribution of
$\operatorname{\mathbf{Var}}[\operatorname{\mathbf{E}}(T\mid
Z_{1})-\operatorname{\mathbf{E}}(T\mid X_{1})]$ to the cross covariances of
$\operatorname{\mathbf{Var}}\operatorname{\mathbf{E}}(T\mid Z_{1})$ is small
(i.e., smaller than the required $(\log s)/s$). The requirement that
$\delta>1/2$, while needed for our theoretical results, may not be needed in
practice. The reasons is that $\delta>1/2$ is required for a uniform bound on
the quantity

$\operatorname{\mathbf{E}}(I_{k}\mid X_{1})\operatorname{\mathbf{E}}(I_{j}\mid
X_{1}),$ (33)

while what the proposition demands is a bound on the expectation. Indeed, in
the extreme case $x=0$ and $\bar{x}=1$, it is easy to see that the expectation
satisfies the required bound even when $\delta\leq 1/2$. In light of this, we
could instead impose the high level condition that

$M(x,\bar{x})=o\biggl{(}\frac{1}{\log^{p}s\cdot s}\biggr{)},$ (34)

where
$M(x,\bar{x})=\operatorname{\mathbf{E}}[\operatorname{\mathbf{E}}(I_{k}\mid
X_{1})\operatorname{\mathbf{E}}(I_{j}\mid X_{1})]$ as above.

####  3.3.1 Bounding
$\operatorname{\mathbf{Var}}\operatorname{\mathbf{E}}(T\mid X_{1})$

We turn next to bound the cross covariances of
$\operatorname{\mathbf{Var}}[\operatorname{\mathbf{E}}(T\mid X_{1})]$. It will
be convenient to change notations, so that $x\mapsto x_{j}$, $\bar{x}\mapsto
x_{k}$, and let $x_{1}$ denote the value of $X_{1}$. Thus, we need to show
that

$\operatorname{\mathbf{E}}[(\operatorname{\mathbf{E}}(T\mid
X_{1}=x_{1})-\mu)(\operatorname{\mathbf{E}}(\bar{T}\mid
X_{1}=x_{1})-\bar{\mu}))]=o\biggl{(}\frac{\log s}{s}\biggr{)}.$ (35)

where $T$ and $\bar{T}$ are the tree estimates at $x$ and $\bar{x}$, and $\mu$
and $\bar{\mu}$ their (unconditional) expectations.

The plan is to show that if $\|x-x_{1}\|_{\infty}$ is small, then

$\operatorname{\mathbf{E}}(T\mid
X_{1}=x_{1})-\mu=\operatorname{\mathbf{E}}(T\mid
X_{1}=x_{1})-\operatorname{\mathbf{E}}(T)$ (36)

is also small. Intuitively, the knowledge of sample point $X_{1}=x_{1}$
changes the expectation of the tree estimate at $x$ only when $x_{1}$ affects
the position of the terminal node containing $x$. This is unlikely to happen
when $X_{1}=x_{1}$ is farther away; as soon as $X_{1}$ leaves the node
containing $x$ in the splitting process, its effect on the terminal node
diminishes.

Toward this end, fix $x$, and let $\Pi$ denote the terminal node that contains
$x$. Since the set of potential splits is unchanged when conditioning on
$X_{1}$, $\Pi$ is a discrete random variable and we may write

$\operatorname{\mathbf{E}}(T)=\sum_{\pi}\operatorname{\mathbf{P}}(\Pi=\pi)\mu_{\pi}\quad\text{and}\quad\operatorname{\mathbf{E}}(T\mid
X_{1}=x_{1})=\sum_{\pi}\operatorname{\mathbf{P}}(\Pi=\pi\mid
X_{1}=x)\mu_{\pi}^{\prime}$ (37)

where $\mu_{\pi}=\operatorname{\mathbf{E}}(T\mid\Pi=\pi)$ and
$\mu_{\pi}^{\prime}=\operatorname{\mathbf{E}}(T\mid\Pi=\pi,X_{1}=x)$.

Recall that $\Pi$ is generated by a recursive splitting procedure, and we can
make a natural correspondence between (37) and the expectation taken over a
directed acyclic graph (DAG) in the following way. Let $[0,1]^{p}$ be the root
of the DAG; for every potential split at $[0,1]^{p}$, there is a directed edge
to a new vertex, where the vertex is the hyperrectangle that contains $x$. If
a vertex is a leaf (i.e., splitting ceases at that node), then the vertex has
no outgoing edges; otherwise, a vertex corresponding to a non-terminal node
has an outgoing edge for each potential split, with each edge going to another
vertex that is the halfspace containing $x$.

To each terminal vertex $v$ in the DAG, which naturally correspond to a
terminal $\pi$ containing $x$, associate a value $f(v)\coloneqq\mu_{\pi}$.
Each edge $e=(v\to w)$ corresponding to a split $s$ from a node $v$ to a
halfspace $w$; associate with this edge the transition probability

$p(e)\coloneqq\operatorname{\mathbf{P}}(\text{$s$ is chosen at
$v$}\mid\text{current node is $v$})\eqqcolon\operatorname{\mathbf{P}}(w\mid
v).$ (38)

Given the transition probabilities, the value $f$ may be extended to each
vertex $v$ via the recursion

$f(v)\coloneqq\sum_{e:v\to w}\operatorname{\mathbf{P}}(w\mid
v)f(w).\qquad\text{Thus, $f(v)$ is the continuation value at $v$.}$ (39)

In this way, the expectation
$\operatorname{\mathbf{E}}(T)=f(\text{root})=f([0,1]^{p})$. If instead we
assign the terminal values $f^{\prime}(v)=\mu^{\prime}_{v}$ and the transition
probabilities

$p^{\prime}(e)=\operatorname{\mathbf{P}}(\text{$s$ is chosen at
$v$}\mid\text{current node is
$v$},X_{1}=x_{1})=\operatorname{\mathbf{P}}^{\prime}(w\mid v),$ (40)

then we recover $\operatorname{\mathbf{E}}(T\mid
X_{1}=x_{1})=f^{\prime}([0,1]^{p})$. Therefore, to bound
$\operatorname{\mathbf{E}}(T\mid X_{1}=x_{1})-\operatorname{\mathbf{E}}(T)$
requires bounding the difference in the continuation values.

We will need to assume that $p^{\prime}(e)\approx p(e)$; that is, We will
require the following assumption regarding differences in the splitting
probabilities.

######  Assumption (Splitting Stability).

For any node $v$, the total variation distance between the distributions
$\\{p(e)\\}_{e:v\to w}$ and $\\{p^{\prime}(e)\\}_{e:v\to w}$ is bounded by the
effective volume of $v$. Specifically, there exists some $\delta>0$ such that
for all $v$,

$\operatorname{TV}(p,p^{\prime})\leq\biggl{(}\frac{1}{s|v|}\biggr{)}^{1+\delta}$
(41)

up to some constant. Here, $|v|$ denotes the volume of the hyperrectangle at
$v$, i.e.,

$|v|=\biggl{|}\prod_{j=1}^{p}(a_{j},b_{j})\biggr{|}=\prod_{j=1}^{p}|b_{j}-a_{j}|.$
(42)

Recall that with probability exceeding $1/s$, the number of sample points
$X_{1},\dots,X_{s}$ in the hyperrectangle $v$ is approximately $s|v|$, from
our uniform distribution assumption. In general, since the density of $X$ is
bounded away from zero and infinity, the number of sample points will be
bounded (above and below) by a constant times $s|v|$. Therefore, the stability
assumption places a restriction on procedure used to select optimal splits:
namely, if the decision is made on the basis of $m$ points, then by moving any
one of the points changes the optimal split with probability bounded by
$m^{-(1+\delta)}$. In practice, most splitting procedures satisfy a stronger
bound. A set of sufficient conditions is given in the following proposition.

######  Proposition 4.

Assume that the optimal split is chosen based on the quantities

$f_{1}(\mu_{1},\dots,\mu_{Q}),\dots,f_{P}(\mu_{1},\dots,\mu_{Q})$ (43)

for some $Q\geq 1$, where $\mu_{1},\dots,\mu_{Q}$ are the sample averages of
the points being split

$\mu_{k}=\frac{1}{n_{v}}\sum_{i:X_{i}\in v}m_{k}(X_{i}).$ (44)

Specifically, the optimal split depends on the which $f_{i}$ achieves the
largest value, i.e., the value $\operatorname*{arg\,max}_{i}f_{i}(\mu)$. Here,
$f_{1},\dots,f_{P}$ are assumed to be Lipschitz, and the functions
$m_{1},\dots,m_{Q}$ are such that $m_{k}(X)$ is 1-subExponential. Then the
splitting stability assumption is satisfied.

######  Remark.

Recall that the $X_{i}$ are bounded—therefore, $m_{k}(X)$ being subExponential
allows the use of squared loss to compute the optimal split.

In general, the conditions in Proposition 4 is sufficient to guarantee the an
exponential bound instead of a polynomial one. Thus, the proposition should be
viewed as simply as providing a “plausibility argument” that stable splitting
rules are common in practice.

The next proposition shows that when the splitting probabilities are bounded
by $(s|v|)^{-(1+\delta)}$ as in the splitting stability assumption, so are the
continuation values.

######  Proposition 5.

Suppose the splitting probabilities satisfy a generic bound

$\operatorname{TV}(p,p^{\prime})\leq\log(s)\Delta(s|v|)\quad\text{at each node
$v$}.$ (45)

Then for any node $v$ containing $x$ but not $x_{1}$,

$|f(v)-f^{\prime}(v)|\leq C\Delta(s|v|)\quad\text{for some constant $C$}.$
(46)

(The constant $C$ does not depend on $v$.)

Finally, we can put the bounds on $\operatorname{TV}(p,p^{\prime})$ and
$|f-f^{\prime}|$ together and prove the required bound on the cross diagonals.

######  Proposition 6.

Suppose that the splitting rule is stable as in (41) and that
$\delta>1-\alpha$. For $x\neq x_{1}$,

$|\operatorname{\mathbf{E}}(T\mid
X_{1}=x_{1})-\operatorname{\mathbf{E}}(T)|=o\biggl{(}\frac{1}{s^{1+\epsilon}}\biggr{)}$
(47)

for some $\epsilon>0$. In particular, the off-diagonal entries of
$\operatorname{\mathbf{Var}}\operatorname{\mathbf{E}}(T\mid X_{1})$ are
$o(s^{-(1+\epsilon)})$.

Note that Proposition 3 required that $\delta>1/2$. Since $\alpha<1/2$ by
definition, this bound is necessarily more restrictive. As before, this
requirement is plausibly looser in applications: we use it to give the bound

$|v|\geq\alpha^{L},$ (48)

even though the RHS is only tight only at most $1/2^{L}$ proportion of the
nodes. Note that the “average” node has area approximately $(1/2)^{L}$, so
that $\delta>1/2$ may be more appropriate.

Putting everything together, we have show that required result that

$\operatorname{\mathbf{Var}}(\mathring{T})_{k,l}=o(s^{-\epsilon})\qquad\text{for
some $\epsilon>0$}.$ (49)

Therefore, by Proposition 2 and the trace calculations, the asymptotic joint
normality of the random forest estimator is established.

##  4 Heuristics and Practical Recommendations

The previous sections established the asymptotic normality result that

$V^{-1/2}(\mathring{\operatorname{RF}}-\mu)\xRightarrow{\;\rm
dist\;}N(0,I),\qquad\text{where
}V=\operatorname{\mathbf{Var}}\mathring{\operatorname{RF}}=\frac{s}{n}\operatorname{\mathbf{Var}}\mathring{T}.$
(50)

Recall $\mu$ is the expectation of $\operatorname{RF}$, while the target
function is $m(x)=\operatorname{\mathbf{E}}(Y\mid X=x)$. Athey and Wager shows
that $(\operatorname{RF}-\mu)/\sqrt{V}\xRightarrow{\;\rm dist\;}N(0,1)$ in the
univariate case; since we have shown that $V$ is diagonally dominant, their
result carries over, so that $V^{-1/2}(\operatorname{RF}-m)\xRightarrow{\;\rm
dist\;}N(0,I)$.

Moreover, (Athey and Wager) proposes a jackknife estimator that can
consistently estimate $\sqrt{V}$ in the scalar case. Our diagonally dominance
result suggests that the random forest estimates at $x$ and
$\bar{x}\in[0,1]^{p}$ are independent in the limit $n\to\infty$,

$\begin{split}\operatorname{\mathbf{Var}}(\operatorname{RF}(x)+\operatorname{RF}(\bar{x}))&=\operatorname{\mathbf{Var}}(\operatorname{RF}(x))+\operatorname{\mathbf{Var}}(\operatorname{RF}(\bar{x}))+2\operatorname{\mathbf{Cov}}(\operatorname{RF}(x)+\operatorname{RF}(\bar{x}))\\\
&\approx\operatorname{\mathbf{Var}}(\operatorname{RF}(x))+\operatorname{\mathbf{Var}}(\operatorname{RF}(\bar{x}))\end{split}$
(51)

so that the jackknife estimator for the scalar case can be used to obtain
confidence intervals. However, the approximation above depends on decay of the
off-diagonal terms. This section, we provide a back of the envelope bound for
the covariance term that may be useful for practioners. We stress that the
following calculations are mostly heuristics: as we have shown above, the
covariance term depends on quantities such as $M(x,\bar{x})$ which is heavily
dependent on the splitting algorithm. Since the aim is to produce a “usable”
result, we focus on heuristics.

To begin, recall that the covariance term is upper bounded by111This is a very
crude upper bound as we have dropped the quantity $\Delta(\alpha^{\ell}s)$
from the infinite series.

$M(x,\bar{x})+\log^{2}s\biggl{(}\sum_{\ell=0}^{\infty}p_{\ell}\biggr{)}\biggl{(}\sum_{\ell=0}^{\infty}\bar{p}_{\ell}\biggr{)}$
(52)

where $p_{\ell}=\operatorname{\mathbf{P}}(L\geq\ell)$ is the probability that
$x$ and $x_{1}$ are not separated after $\ell$ splits (c.f. proof of
Proposition 6) and likewise for $\bar{p}_{\ell}$. If we let $I$ (resp.
$\bar{I}$) denote the indicator that $X_{1}$ is in the terminal node of $x$
(resp. $\bar{x}$), then $I=1$ is equal to the event that $L=\log_{2}s$, so
that

$\operatorname{\mathbf{E}}(I\mid X_{1}=x_{1})=\operatorname{\mathbf{P}}(L=\log
s)\leq\frac{\operatorname{\mathbf{E}}L}{\log s}.$ (53)

Replacing the inequality with an approximation, we have
$\operatorname{\mathbf{E}}L=(\log s)\operatorname{\mathbf{E}}(I\mid
X_{1}=x_{1})$. All of this shows that the covariance term is bounded by

$(\log^{4}s)M(x,\bar{x}),\quad\text{up to appropriate constants}.$ (54)

Towards a useful heuristic, we will consider a bound on the correlation
instead of the covariance. In our notation, the result of (Athey and Wager)
lower bounds $M(x,x)$ (and $M(\bar{x},\bar{x})$) in their proof, while we
upper bound $M(x,\bar{x})$. Ignoring the logarithmic terms, we have

$\biggl{|}\frac{\operatorname{\mathbf{Cov}}(T(x),T(\bar{x}))}{\sqrt{\operatorname{\mathbf{Var}}T(x)\cdot\operatorname{\mathbf{Var}}T(\bar{x})}}\biggr{|}\approx\frac{M(x,\bar{x})}{\sqrt{M(x,x)M(\bar{x},\bar{x})}}.$
(55)

Recall that
$M(x,\bar{x})=\operatorname{\mathbf{E}}[\operatorname{\mathbf{E}}(I\mid
X_{1})\operatorname{\mathbf{E}}(\bar{I}\mid X_{1})]$ and decays superlinearly
as $\bar{x}$ moves away from $x$. Using the previous expression (note that
$M(x,x)\approx M(\bar{x},\bar{x})$ due to symmetry between $x$ and $\bar{x}$),
we can bound the correlation from purely geometric considerations. Since the
integrand

$\operatorname{\mathbf{E}}(I\mid X_{1})\operatorname{\mathbf{E}}(\bar{I}\mid
X_{1})$ (56)

decays as $X_{1}$ moves away from $x$ (and $\bar{x}$), we may imagine that its
integral $M(x,x)$ comes from contributions to points $X_{1}$ near $x$, say
those points in a $L_{\infty}$-box of side lengths $d$ with volume $d^{p}$,
i.e., points $\\{y\in[0,1]^{p}:\|x-y\|\leq d/2\\}$. For $M(x,\bar{x})$, the
contribution would come from points that are with $d/2$ of both $x$ _and_
$\bar{x}$, and to a first degree approximation, the volume of these points
$\\{y\in[0,1]^{p}:\|x-y\|_{\infty}\leq d/2,\|\bar{x}-y\|\leq d/2\\}$ is

$(d-z_{1})\dots(d-z_{p})\approx
d^{p}-(z_{1}+\dots+z_{p})d^{p-1},\quad\text{where
$z_{i}=|x_{i}-\bar{x}_{i}|$},$ (57)

where the approximation is accurate if $|z_{i}|\ll 1$. The proportion of the
volume of these points is therefore $1-\frac{1}{d}\|x-\bar{x}\|_{1}$, which
leads to the heuristic

$\biggl{|}\frac{\operatorname{\mathbf{Cov}}(T(x),T(\bar{x}))}{\sqrt{\operatorname{\mathbf{Var}}T(x)\cdot\operatorname{\mathbf{Var}}T(\bar{x})}}\biggr{|}\approx
1-c\|x-\bar{x}\|_{1},\qquad\text{for some constant $c$}.$ (58)

The RHS has the correct scaling when $x=\bar{x}$, i.e., the correlation equals
one when $\|x-\bar{x}\|_{1}=0$. To maintain correct scalinag at the other
extreme with $\|x-\bar{x}\|_{1}=p$, we should take $c=1/p$, so that

$\biggl{|}\frac{\operatorname{\mathbf{Cov}}(T(x),T(\bar{x}))}{\sqrt{\operatorname{\mathbf{Var}}T(x)\cdot\operatorname{\mathbf{Var}}T(\bar{x})}}\biggr{|}\approx
1-\frac{1}{p}\sum_{i=1}^{p}|x_{i}-\bar{x}_{i}|.$ (59)

Of course, this heurstic is “wrong” in that it does not depend on $s$; our
theoretical results show that even for non-diametrically opposed points, the
correlation drops to zero as $s\to\infty$. Therefore, another recommendation
is to use

$\biggl{|}\frac{\operatorname{\mathbf{Cov}}(T(x),T(\bar{x}))}{\sqrt{\operatorname{\mathbf{Var}}T(x)\cdot\operatorname{\mathbf{Var}}T(\bar{x})}}\biggr{|}\approx\min\biggl{(}1-\frac{s^{\epsilon}}{p}\sum_{i=1}^{p}|x_{i}-\bar{x}_{i}|,\;0\biggr{)},$
(60)

for some $\epsilon>0$, where dependence $s^{\epsilon}$ comes from matching the
heuristic $c=1/d$ above with the proof of Proposition 3.

##  5 Simulations

In this section, we conduct numerical experiments of varying sample sizes to
gauge the empirical relevance of our theoretical results and heuristics. In
our experiments, we set $p=2$, so that the covariates $X$ are distributed on
the unit square. Instead of adopting a uniform distribution, for ease of
interpretation, the distribution of $X$ is

$X\sim\begin{cases}\bar{N}(\mu_{1},I_{2})&\text{with probability $1/4$}\\\
\bar{N}(\mu_{2},I_{2})&\text{with probability $1/4$}\\\
\bar{N}(\mu_{3},I_{2})&\text{with probability $1/4$}\\\
\bar{N}(\mu_{4},I_{2})&\text{with probability
$1/4$}\end{cases}\quad\text{where}\quad\begin{cases}\mu_{1}=(0.3,0.3)^{\intercal}\\\
\mu_{2}=(0.3,0.7)^{\intercal}\\\ \mu_{3}=(0.7,0.3)^{\intercal}\\\
\mu_{4}=(0.7,0.7)^{\intercal}\\\ \end{cases}$ (61)

and $\bar{N}$ denotes a truncated multivariate Gaussian distribution on the
unit square.222 That is, $\bar{N}(\mu,\Sigma)$ denotes the conditional
distribution of $x\sim N(\mu,\Sigma)$ on the event $x\in[0,1]^{2}$. Thus, $X$
has a bounded density on the unit square, and has four peaks at
$\mu_{1},\dots,\mu_{4}$. The distribution of $Y$ conditional on
$X=(x_{1},x_{2})$ is

$Y\sim\frac{x_{1}+x_{2}}{2}+\frac{1}{5}N(0,1).$ (62)

The random splitting probability is $\delta=1/2$, and the regularity
parameters are $\alpha=0.01$ and $k=1$, so that the tree is grown to the
fullest extent (i.e., terminal nodes may contain a single observation), with
each terminal node lying on the $101\times 101$ grid of the unit square. For
each sample size $n$, five thousand trees are grown, and the estimates are
aggregated to compute the correlation.

Figure 1 plots the correlation of estimates at $x$ and $\bar{x}$ as a function
of the $L_{1}$ norm $\|x-\bar{x}\|_{1}$. The calculation is performed by first
fixing $x$, then calculating the sample correlation (across five thousand
trees) as $\bar{x}$ ranges over each cell: the correlation is associated with
the $L_{1}$ norm $\|x-\bar{x}\|_{1}$. This process is then repeated by varying
the reference point $x$, and the correlation at $\|x-\bar{x}\|_{1}$ is the
average of the correlations observed.

Figure 1: Correlation as a function of sample size and $L_{1}$ norm.

The figure demonsrates that the linear heuristic (60) given in the previous
section is conservative: it is evident that correlation decreases super-
linearly as $x$ and $\bar{x}$ become separated. Figure 2 plots the correlation
on a logarithmic scale, which shows that that correlation decay is exponential
in a neighborhood of unity. In other words, simulations suggest that the
correct heuristic is of the shape

$\biggl{|}\frac{\operatorname{\mathbf{Cov}}(T(x),T(\bar{x}))}{\sqrt{\operatorname{\mathbf{Var}}T(x)\cdot\operatorname{\mathbf{Var}}T(\bar{x})}}\biggr{|}\approx
e^{-\lambda\|x-\bar{x}\|_{1}}$ (63)

for some suitable $\lambda$. Note, however, that the rate of decay of the
correlation depends heavily on the exact splitting algorithm employed. Thus,
an avenue for future empirical work is to investigate the decay behavior of
popular random forest algorithms used in practice.

Figure 2: Logarithm of correlation as a function of sample size and $L_{1}$
norm.

##  6 Conclusion and Points for Future Work

Regression trees and random forests are popular and effective non-parametric
estimators in practical applications. A recent paper by Athey and Wager shows
that the random forest estimate at any point is asymptotically Gaussian; in
this paper, we extend this result to the multivariate case and show that the
vector of estimates at multiple points is jointly normal. Specifically, the
covariance matrix of the limiting normal distribution is diagonal, so that the
estimates at any two points are independent in sufficiently deep trees.
Moreover, the off-diagonal term is bounded by quantities capturing how likely
two points belong to the same partition of the resulting tree. Certainly
stability properties of the base learner are central to our results, c.f.
Propositions 3 and 4. We test our proposed covariance bound and the associated
coverage rates of confidence intervals in numerical simulations.

This paper provides the a theoretical basis for performing inferences on
functionals of target function (e.g., a heterogeneous treatment effect) when
the functional is based on values of the target function at multiple points in
the feature space. Specifically, we show that the covariance vanishes in the
limit relative to the variance, and provide heuristics on the size of the
correlation in finite samples.

There are a couple interesting avenues for future work: the first is to extend
our framework to cover categorical or discrete-valued features. Here, new
assumptions would be required in order to maintain the guarantee that node
sizes are “not too small.” Another direction is to use our bounds—with
possible improvements—on the covariance matrix of the underlying $U$-statistic
and leverage recent results [6, 7] in order to provite finite sample bounds.
This would provide a more sound theoretical underpinning of our heuristics.

##  Appendix A Proofs

###### Proof of Proposition 1.

For random vectors in $\mathbf{R}^{q}$, define the inner product

$\langle X,Y\rangle\coloneqq\operatorname{\mathbf{E}}(X^{\intercal}MY).$ (64)

For each subset $A\subseteq\\{1,\dots,n\\}$, let $H_{A}$ be the set of square-
integrable random vectors of the form $g(X_{i}:i\in A)$, where $g$ is a
function of $|A|$ arguments. It is easy to see that $H_{A}$ are pairwise
orthogonal, by induction on $r=|A|$, the direct sum $\bigoplus_{B\subseteq
A}H_{B}$ is equal to set of all random variables which are functions of
$\\{X_{i}:i\in A\\}$. The required decomposition follows since $\\{X_{i}\\}$
are IID. ∎

###### Proof of Proposition 2.

We will prove the slightly more general statement that if $A_{n}$ and $B_{n}$
are two sequences of square matrices with bounded entries such that

$B_{ii}>\delta\text{ for some $\delta$ for all $n$}\quad\text{and}\quad
A_{ii}\geq\frac{B_{ii}}{\log n}$ (65)

and $A_{ij}=o(1/\log n)$, then $\operatorname{tr}(A^{-1}B)\to 0$. To prove
this, note that $|\det A|=\prod_{\pi}a_{i\pi_{i}}$ over permutations $\pi$;
since the off-diagonal entries of $A_{ij}$ vanish relative to $A_{ii}$, we
have $|\det A|\sim\prod A_{ii}$, where the notation $a\sim b$ stands for
$c|b|\leq|a|\leq c^{\prime}|b|$ for constants $c$ and $c$’ not depending on
$n$.

Next, recall that Cramer’s Rule yields the formula $(A^{-1})_{ii}=\det
A_{-i}/\det A$, where $A_{-i}$ is the matrix $A$ with its $i$-th row and
$i$-th column removed. By the same token, $|\det A_{-i}|\sim\prod_{j\neq
i}A_{jj}$. Therefore,

$(A^{-1})_{ii}\sim\frac{1}{A_{ii}}.$ (66)

In particular, the $i$-th diagonal entry of the matrix $A^{-1}B$ is given by

$(A^{-1}B)_{ii}=(A^{-1})_{ii}B_{ii}+\sum_{j\neq
i}(A^{-1})_{ij}B_{ji}\sim\frac{B_{ii}}{A_{ii}}\leq\log n.$ (67)

Therefore, whenever ∎

###### Proof of Proposition 3.

Recall that the splitting algorithm has a probability $\delta$ chance of
splitting on the $j$-th axis. Since each terminal node contains a constant
number of points, the number of terminal nodes is equal (up to constant) to
$s$, the number of training samples. Therefore, the number of splits to a
terminal node is equal to $\log_{2}s=\log s$.

Let $c_{j}(x)$ denote the number of splits in coordinate $j$ along the
sequence of splits leading to the terminal node containing $x$. By our
randomization assumption, each split has at least an independent chance
$\delta$ of being chosen, so that

$c_{j}(x)\succeq\frac{1}{p}B\biggl{(}\log\frac{s}{K},\delta\biggr{)}\qquad\text{where
$\succeq$ stands for stochastic dominance}.$ (68)

Recall that each split along the $j$-th axis decreases the length of the
$j$-th dimension of the hyper-rectangle by a factor at least $(1-\alpha)$.
Since $x\neq\bar{x}$, there exists $j$ for which $x_{j}\neq\bar{x}_{j}$ and we
may assume WLOG that $x_{1}\neq\bar{x}_{1}$. For $x$ and $\bar{x}$ to belong
in the same terminal node, it is necessary that the length of the terminal
node along the first axis is at least $|x_{1}-\bar{x}_{1}|$; since we start in
the unit hypercube,

$(1-\alpha)^{c_{1}(x)}\geq|x_{1}-\bar{x}_{1}|\implies
c_{1}(x)\leq\frac{\log|x_{1}-\bar{x}_{1}|}{\log(1-\alpha)}\eqqcolon\rho,$ (69)

where $\rho$ is a constant not depending $s$ nor $\delta$. Therefore,

$\operatorname{\mathbf{P}}(\text{$x$ and $\bar{x}$ lie in the same terminal
node})\leq\operatorname{\mathbf{P}}\biggl{[}B\biggl{(}\log\frac{s}{K},\delta\biggr{)}\leq
p\rho\biggr{]}.$ (70)

Since $p\rho$ is a constant, it is not hard to see that

$\operatorname{\mathbf{P}}\biggl{[}B\biggl{(}\log\frac{s}{K},\delta\biggr{)}\leq
p\rho\biggr{]}\leq(1-\delta+o(1))^{\log
s/K}=\biggl{(}\frac{1}{s}\biggr{)}^{\log\frac{1}{1-\delta+o(1)}}.$ (71)

The base of the logarithm is two since the tree is binary. Therefore, if we
choose $\delta>1/2$, the exponent exeeds $1$ and the proof is complete. ∎

###### Proof of Proposition 4.

The easiest case is the splitting decision in the root node $[0,1]^{p}$, so we
start here. We prove the result by introducing a coupling between the split
decisions with and without conditioning on $X_{1}=x_{1}$

$\begin{split}S&=\operatorname*{arg\,max}_{i}f_{i}\biggl{(}\frac{1}{s}\sum_{i=1}^{s}m_{1}(X_{i}),\dots,\frac{1}{s}\sum_{i=1}^{s}m_{Q}(X_{i})\biggr{)}\eqqcolon
f_{i}\\\
S^{\prime}&=\operatorname*{arg\,max}_{i}f_{i}\biggl{(}\frac{1}{s}\biggl{[}m_{1}(x_{1})+\sum_{i=2}^{s}m_{1}(X_{i})\biggr{]},\dots,\frac{1}{s}\biggl{[}m_{Q}(x_{1})+\sum_{i=2}^{s}m_{Q}(X_{i})\biggr{]}\biggr{)}\eqqcolon
f_{i}^{\prime}.\end{split}$ (72)

Here, $S$ is the split made on the sample $X_{1},\dots,X_{s}$ and $S^{\prime}$
is the split made on the sample conditional on $X_{1}=x_{1}$. Note that we may
WLOG assume that the splits are not randomly chosen, since the splitting
probabilities are trivially equal conditional on that event. Since
$X_{2},\dots,X_{s}$ is independent of $X_{1}$, this is a valid coupling of the
two splitting processes. Clearly, a necessary condition for $S\neq S^{\prime}$
is the existence of a pair $1\leq i\neq j\leq P$ for which

$f_{i}>f_{j}\quad\text{but}\quad f_{j}^{\prime}>f_{i}.$ (73)

Due to our Lipschitz and subExponential assumptions, it is clear that $f_{i}$
and $f_{j}$ their respective limits
$f(\operatorname{\mathbf{E}}m_{1},\dots,\operatorname{\mathbf{E}}m_{Q})$;
hence, whenever
$f_{i}(\operatorname{\mathbf{E}}m_{1},\dots,\operatorname{\mathbf{E}}m_{Q})>f_{j}(\operatorname{\mathbf{E}}m_{1},\dots,\operatorname{\mathbf{E}}m_{Q})$
we will have

$f_{i}-f_{j}>\frac{1}{s}\quad\text{with probability at least $1-O(e^{-cs})$
for some constant $c$.}$ (74)

However, the difference of the arguments of $f_{i}$ in (72) differ by at most
$1/s$ by Lipschitz continuity. Therefore,
$(f_{i}>f_{j},f_{j}^{\prime}>f_{i}^{\prime})$ occurs with probability at most
$1-O(e^{-cs})$, and we finish by noting taking a union over the $\binom{P}{2}$
pairs $(i,j)$.

Note that the above actually proves something stronger, namely that for every
split $\tau$,

$\operatorname{\mathbf{P}}(S\neq s\mid
S^{\prime}=\tau)<e^{-cs}\quad\text{and}\quad\operatorname{\mathbf{P}}(S^{\prime}\neq
s\mid S=\tau)<e^{-cs}.$ (75)

It follows that see that for any $\tau$, the total variation distance between
$(X_{2},\dots,X_{s}\mid S=\tau)$ and $(X_{2},\dots,X_{s}\mid S^{\prime}=\tau)$
is at most $e^{-cs}$. To see this, note that $S$ and $S^{\prime}$ are
functions of $X_{i}$ only, so that the densities of these two distributions
are

$p(x)=\operatorname{\mathbf{1}}(S(x)=\tau)\frac{p(x)}{\operatorname{\mathbf{P}}(S=\tau)}\quad\text{and}\quad
p^{\prime}(x)=\operatorname{\mathbf{1}}(S^{\prime}(x)=\tau)\frac{p(x)}{\operatorname{\mathbf{P}}(S^{\prime}=\tau)}$
(76)

respectively. We may assume WLOG that
$\operatorname{\mathbf{P}}(S^{\prime}=\tau)\geq\operatorname{\mathbf{P}}(S=\tau)$
so that the total variation is

$\begin{split}\int|p(x)-p^{\prime}(x)|&=\int_{S=\tau}p(x)-p^{\prime}(x)+\int_{S\neq\tau,S^{\prime}=\tau}p^{\prime}(x)\\\
&=1-\frac{\operatorname{\mathbf{P}}(S^{\prime}=\tau,S=\tau)}{\operatorname{\mathbf{P}}(S^{\prime}=\tau)}+\frac{\operatorname{\mathbf{P}}(S^{\prime}=\tau,S\neq\tau)}{\operatorname{\mathbf{P}}(S^{\prime}=\tau)}=2\operatorname{\mathbf{P}}(S\neq\tau\mid
S^{\prime}=\tau)<e^{-cs}.\end{split}$ (77)

The upshot is that when considering the splitting probability in the next
node, we can ignore the difference in the distribution of $X_{2},\dots,X_{s}$
when conditioning on $S=\tau$ versus conditioning on $S^{\prime}=\tau$ and pay
a cost $O(e^{-cs})$.

Now consider bounding the difference of the spltting probabilities at the next
split

$\operatorname{\mathbf{P}}(S_{2}=s\mid
S_{1}=\tau)-\operatorname{\mathbf{P}}(S_{2}=s\mid S_{1}=\tau,X_{1}=x_{1}).$
(78)

Again, the strategy is to find a coupling $(S_{2},S_{2}^{\prime})$ such that

$S_{2}\sim(S_{2}\mid S_{1}=\tau)\quad\text{and}\quad
S_{2}^{\prime}\sim(S_{2}\mid S_{1}=\tau,X_{1}=x_{1})$ (79)

with $\operatorname{TV}(S_{2},S_{2}^{\prime})\leq e^{-s|v|}$, $|v|$ being the
hyper-rectangle corresponding to one of the halfspaces produced by $\tau$.
Since the distribution of $X_{1},\dots,X_{s}$ conditional on $S_{1}=\tau$
differs from its unconditional distribution by an amount $e^{-cs}$ in the
total variation distance, we could use the following coupling

$\begin{split}S_{2}&=\operatorname*{arg\,max}f_{i}\biggl{(}\frac{1}{n_{v}}\sum_{X_{i}\in
v}m_{1}(X_{i}),\dots,\frac{1}{s|v|}\sum_{X_{i}\in v}m_{Q}(X_{i})\biggr{)}\\\
S_{2}^{\prime}&=\operatorname*{arg\,max}f_{i}\biggl{(}\frac{1}{n_{v}}\sum_{X_{i}^{\prime}\in
v}m_{1}(X_{i}^{\prime}),\dots,\frac{1}{s|v|}\sum_{X_{i}\in
v}m_{Q}(X_{i}^{\prime})\biggr{)}\end{split}$ (80)

where $X_{i}$ follows the distribution of $(X_{1},\dots,X_{s})$ conditional on
$S_{1}=\tau$ and $X_{i}^{\prime}$ follows the distribution conditional on
$S_{1}=\tau$ and $X_{1}=x_{1}$. By conditioning on $S_{1}=\tau$ instead of
$S_{1}^{\prime}=\tau$, we increase the total variation by an amount $e^{-cs}$
via the triangle inequality.

Now the rest of the proof is the same as in the base case, noting that with
high probability, the number $n_{v}$ of points in $v$ is equal to $s|v|$ up to
an multiplicative constant $(1-\eta)$ with probability $e^{-s\eta^{2}}$. The
previous bounds are applied recursively at each depth $l$ of the DAG. At depth
$l$, we incur an “approximation cost” from the total variation distance
bounded by $e^{-|v|s}$. Since $s|v|\geq\alpha^{l}$, it follows that $l\leq
O(\log s|v|)$, whence the cumulative cost at depth $l$ is $O(\log s|v|\cdot
e^{-c|v|s})$. Putting everything together, we have proven that

$\operatorname{TV}(p,p^{\prime})\leq O(\log(|v|s)\cdot e^{-c|v|s})\leq
o\biggl{(}\frac{1}{s|v|}\biggr{)}^{1+\delta}\quad\text{for some
$\delta>0$}.\qed$ (81)

###### Proof of Proposition 5.

The claim is trivially true (by choosing an appropriate constant) if $v$ is a
terminal node. Fix a node $v$ such that $x_{1}\notin v$ and

$\mathbf{X}=\mathbf{X}_{v}=\\{X_{i}:X_{i}\in v\\}$ (82)

be the set of points landing in $v$, so that
$k\coloneqq|\mathbf{X}|\in\\{1,\dots,n-1\\}$. Clearly, $f$ and $f^{\prime}$
depend on the distribution of $\mathbf{X}$ conditional on $\Pi=v$ and
$\bar{\Pi}=v$, where $\Pi$ and $\Pi^{\prime}$ are the splits with splitting
probabilities $p$ and $p^{\prime}$, respectively. We will show that, for each
$k\in\\{1,\dots,n-1\\}$, the total variation distance of

$(\mathbf{X}\mid|\mathbf{X}|=k,\Pi=v)\quad\text{and}\quad(\mathbf{X}\mid|\mathbf{X}|=k,\Pi^{\prime}=v,X_{1}=x_{1})$
(83)

is bounded by $(\log s)\cdot\Delta(s|v|)$. First, since $x_{1}\notin v$,
$X_{1}$ is not an element, $\mathbf{X}$ and $X_{1}$ are independent
conditional on $\Pi$. Therefore,

$\operatorname{dist}(\mathbf{X}\mid|\mathbf{X}|=k,\Pi^{\prime}=v,X_{1}=x_{1})=\operatorname{dist}(\mathbf{X}\mid|\mathbf{X}|=k,\Pi^{\prime}=v).$
(84)

The depth of $v$ is at most $\log_{2}s$, so by the stability assumption,
$\operatorname{\mathbf{P}}(\Pi=\Pi^{\prime})\geq 1-(\log s)\Delta(s|v|)$.
Since the total variation distance of the distributions in (83) differs by
$\operatorname{\mathbf{P}}(\Pi\neq\Pi^{\prime})$, the result follows. ∎

###### Poor of Proposition 6.

The idea is to recursively expand the formulas $\operatorname{\mathbf{E}}T$
and $\operatorname{\mathbf{E}}(T\mid X_{1})$ in terms of the DAG. We start
with

$\begin{split}|\operatorname{\mathbf{E}}T-\operatorname{\mathbf{E}}(T\mid
X_{1}=x_{1})|&=\biggl{|}\sum_{v}p(e)f(e)-\sum_{v}p^{\prime}(e)f^{\prime}(e)\biggr{|}\\\
&\leq\sum_{v}|p(e)-p^{\prime}(e)|f^{\prime}(e)+\biggl{|}\sum_{v}p^{\prime}(e)(f(e)-f^{\prime}(e))\biggr{|},\end{split}$
(85)

where the sum runs over nodes $v$ after the first split, i.e., $[0,1]^{p}\to
v$. The second summand may be split into two, one over $x_{1}\in v$ and the
other $x_{1}\notin v$. If we assume splitting stability held with function
$\Delta$, then Proposition 5 allows us to bound the second term, so that

$\begin{split}|\operatorname{\mathbf{E}}T-\operatorname{\mathbf{E}}(T\mid
X_{1}=x_{1})|&\leq\sum_{v}|p(e)-p^{\prime}(e)|f^{\prime}(e)+\biggl{|}\sum_{v}p^{\prime}(e)(f(e)-f^{\prime}(e))\biggr{|}\\\
&\leq\Delta(\alpha s)+(\log s)\Delta(\alpha s)+\sum_{x_{1}\in
v}p^{\prime}(e)|f(e)-f^{\prime}(e)|,\end{split}$ (86)

where used the fact that $|v|\geq\alpha$. Now, each of the terms
$|f(e)-f^{\prime}(e)|$ may be bounded by
$\Delta(\alpha^{2}s)+\log(s)\Delta(\alpha^{2}s)+\sum_{x_{1}\in w}(\cdots)$.
Continuing in this way, we have

$|\operatorname{\mathbf{E}}T-\operatorname{\mathbf{E}}(T\mid
X_{1}=x_{1})|\leq\log(s)(\Delta(s)+p_{1}\Delta(\alpha
s)+p_{2}\Delta(\alpha^{2}s)+\dots)=\log
s\sum_{\ell=0}^{\infty}p_{\ell}\Delta(\alpha^{\ell}s),$ (87)

where $p_{\ell}$ is the probability that $x$ and $x_{\ell}$ belong to the same
node after $\ell$ splits. In other words, if we let $L$ be the number of
splits after which $x$ and $x_{1}$ are separated, then

$p_{\ell}=\operatorname{\mathbf{P}}(L\geq\ell).$ (88)

Since $x\neq x_{1}$, we may assume without loss of generality that
$\|x-x_{1}\|>\beta$ for some fixed $\beta$ (c.f. the proof of Proposition 3).
In particular,

$p_{\ell}\leq(1-\delta+o(1))^{\ell}$ (89)

for sufficiently large $\ell$. Moreover,
$\Delta(\alpha^{l}s)=\frac{1}{s^{1+\epsilon}}(\frac{1}{\alpha^{1+\epsilon}})^{\ell}$,
whennce $\delta>1-\alpha^{1+\delta}$ is enough to ensure that infinite series
is less than $\frac{1}{s^{1+\epsilon}}$. In particular, as $s\to 0$, we may
take $\epsilon\to 0$, so that the restriction is satisfied (after suitable
constants) by $\delta>1-\alpha$. This completes the proof. ∎

## References

  * [1] Leo Breiman.  Random forests.  Machine Learning, 2001.
  * [2] Stefan Wager and Susan Athey.  Estimation and inference of heterogeneous treatment effects using random forests.  Journal of the American Statistical Association, 113(523):1228–1242, 2018.
  * [3] S. Wager and Guenther Walther.  Adaptive concentration of regression trees, with application to random forests.  arXiv: Statistics Theory, 2015.
  * [4] Susan Athey, Julie Tibshirani, and Stefan Wager.  Generalized random forests.  Ann. Statist., 47(2):1148–1178, 04 2019.
  * [5] Nino Arsov, Martin Pavlovski, and Ljupco Kocarev.  Stability of decision trees and logistic regression, 2019.
  * [6] Victor Chernozhukov, Denis Chetverikov, and Kengo Kato.  Central limit theorems and bootstrap in high dimensions.  Ann. Probab., 45(4):2309–2352, 07 2017.
  * [7] Xiaohui Chen.  Gaussian and bootstrap approximations for high-dimensional u-statistics and their applications.  Ann. Statist., 46(2):642–678, 04 2018.