# Data-adaptive exposure thresholds for the Horvitz-Thompson estimator of the
Average Treatment Effect in experiments with network interference

Vydhourie Thiyageswaran
Department of Statistics
University of Washington
&Tyler McCormick
Department of Statistics
University of Washington
&Jennifer Brennan
Google Research

###### Abstract

Randomized controlled trials often suffer from interference, a violation of
the Stable Unit Treatment Values Assumption (SUTVA) in which a unit’s
treatment assignment affects the outcomes of its neighbors. This interference
causes bias in naive estimators of the average treatment effect (ATE). A
popular method to achieve unbiasedness is to pair the Horvitz-Thompson
estimator of the ATE with a known exposure mapping: a function that identifies
which units in a given randomization are not subject to interference. For
example, an exposure mapping can specify that any unit with at least
$h$-fraction of its neighbors having the same treatment status does not
experience interference. However, this threshold $h$ is difficult to elicit
from domain experts, and a misspecified threshold can induce bias. In this
work, we propose a data-adaptive method to select the "$h$"-fraction threshold
that minimizes the mean squared error of the Hortvitz-Thompson estimator. Our
method estimates the bias and variance of the Horvitz-Thompson estimator under
different thresholds using a linear dose-response model of the potential
outcomes. We present simulations illustrating that our method improves upon
non-adaptive choices of the threshold. We further illustrate the performance
of our estimator by running experiments on a publicly-available Amazon product
similarity graph. Furthermore, we demonstrate that our method is robust to
deviations from the linear potential outcomes model.

##  1 Introduction

Under network interference, where the Stable Unit Treatment Values Assumption
(SUTVA) is violated, estimating average treatment effects incorporating direct
effects (impact of the treatment on the respondent’s outcome) and indirect
effects (impact of treated peers on the respondent’s outcome) is a central
problem in learning causal effects. We study the setting of randomized
controlled trials. For example, estimating the adoption of a new idea through
random assignment of people to advertisements is complicated by the
dissemination of information through social interactions. In settings where
user behavior informs a prediction algorithm, the behavior of other users
creates spillovers through their impact on the prediction algorithm. As such
discrete interactions can be modeled through network interference models, the
problem of estimating average treatment effects under network interference
becomes a ubiquitous one.

A major challenge arises from the lack of knowledge of the exact interference
structure that contributes to the indirect effect. For some treatments, having
treated peers might have very little impact (e.g. a drug for a noncommunicable
disease that’s only available as part of a clinical trial), whereas for others
the impact could be extremely large (e.g. vaccines). _Exposure maps_ classify
individuals into groups with similar contact patterns with treated peers.

One simple form of exposure map classifies individuals as indirectly exposed
if more than some fixed fraction $h$ of their network is treated, meaning that
the researcher expects contact with a smaller fraction than $h$ to have no
impact on the outcome and a fraction lager than $h$ to have a substantial, but
relatively consistent impact as the fraction increases above the threshold.
Our setting corresponds to what Centola and Macy (2007) refer to as _contested
contagion._ That is, settings where _both_ treated and untreated peers
influence an individual’s outcome, but in different directions, and the size
of the respondents’ network also matters. Centola and Macy (2007) give an
example of cleaning litter from a beach–a person might be willing to pick up
litter alone if 10 people use the beach, but the task would seem futile
without a group of civically-minded individuals if 1,000 people use the beach.
This conception of exposure is particularly salient in studying the adoption
of technology, particularly in settings where the technology becomes more
useful as more members of the network adopt it, known as the “threshold of
adoption" (Acemoglu et al., 2011; Reich, 2016). Adopting a messenger app, for
example, becomes more appealing as a larger fraction of social contacts also
use the app.

Given a known $h$-fractional exposure mapping, the Horvitz-Thompson estimator
under this criterion is commonly used as an unbiased estimator of the average
treatment effect. Without domain knowledge of the interference structure,
however, an estimator with fixed treatment exposure conditions for this
estimand will be high in either bias or variance. One could, for example,
resort to using the Horvitz-Thompson estimator at extreme $h$s for
unbiasedness, such as when the unit in question, and all its neighbors are
treated. However, this incurs the cost of a large variance, as the estimator
scales inversely with the probability of the unit and all its neighbors being
treated. This probability can be very small for high-degree nodes,
particularly, in dense graphs and when treatment fractions are small. Choosing
a small threshold, however, means that many individuals who are in reality not
impacted indirectly by the treatment will be counted amongst the indirectly
exposed, biasing effect estimates. This motivates the study of mean-squared-
error (MSE)-minimizing exposure threshold selection to combine with existing
estimators.

Our contribution to this setting is to propose a simple, data-dependent
approach to estimating the average treatment effect optimally in terms of the
MSE. We compute an approximate rate of change of bias and variance with
respect to the threshold to inform the MSE-optimal threshold for the Horvitz-
Thompson estimator. In particular, we propose using linear regression for an
approximate bias-rate of change measure as a simple and intuitive approach.
Additionally, this study allows us to understand graph structures under which
this approximation is more precise.

###  1.1 Related Works

The exposure mapping framework is a common strategy for quantifying
spillovers. In Aronow and Samii (2017); Ugander et al. (2013); Sussman and
Airoldi (2017); Hardy et al. (2019), and Auerbach and Tabord-Meehan (2021),
the authors expounded upon different exposure mappings, and estimation under
these settings. Eckles et al. (2017); Toulis and Kao (2013) assume a
(weighted) fraction of treated neighbors in characterizing the exposure. We
work under a similar setting. Other common assumptions include exposures
described by the raw count numbers of the treated/control neighbors Ugander et
al. (2013). In Basse and Airoldi (2018); Cai et al. (2015), a (known)
generalized linear form of network effects on neighborhood treatments was
assumed. We do not assume this. Related to our work, in Zhu et al. (2024), the
authors fit functionals on the treatment assignment and exposure.

For off-policy evaluation in reinforcement learning, Su et al. (2020)
investigate adaptive MSE-optimal kernel bandwidth selection (a widely studied
problem, see for example, Fan and Gijbels (1992), Ruppert (1997) Kallus and
Zhou (2018)) using Lepski’s method Goldenshluger and Lepski (2011). To the
best of our knowledge, a reframing of our problem in terms of optimal
bandwidth selection A.3 has not been studied in the average treatment effect
estimation under network interference setting. The authors of Belloni et al.
(2022) have a different, but related, goal of estimating the direct effect
under network interference. Their approach was to first estimate interference
structure by estimating the $m_{i}$-hop influence (beyond the immediate graph
neighborhood), or exposure radius under assumptions of additivity of the main
effects in the potential outcome model.

###  1.2 Notation.

We use $h$ to denote the threshold for “treatment exposure", and $1-h$ for
“control exposure." Therefore, for larger $h$, i.e. $h$ closer to one, we have
a more restrictive setting, where only the subset of treated (resp. control)
units with most of their neighbors treated (resp. control) are counted as
treatment (resp. control) exposed. For smaller $h$, i.e. $h$ close to zero, we
have a less restrictive setting, as a larger subset of treated (resp. control)
units satisfy the threshold condition. Throughout the paper we use $d_{i}$ for
the degree of node $i$, and $d$ when the graph is regular.

##  2 Setup

We consider a fractional exposure mapping. That is, for any treatment
assignment vector $\mathbf{Z}\in\\{0,1\\}^{n}$, the effective influence of the
graph’s treatment assignment on the outcome of node $i$ is equivalent to the
influence of fractional treatment assignments in the neighborhood of node $i$.
This reduces our potential outcome model that so that it takes the following
form:

$\displaystyle
y_{i}(\mathbf{Z})=y_{i}(z_{i},e_{i})=\alpha+g(z_{i})+f(e_{i})+\epsilon_{i}\equiv\alpha+\psi(z_{i},e_{i})+\epsilon_{i},$
(1)

where $z_{i}\in\\{0,1\\}$ is the treatment assignment of unit $i$,
$e_{i}\in[0,1]$ is the fraction of treated neighbors of unit $i$ (not
including $i$ itself), and $\epsilon_{i}$ are independent (and independent of
$(z_{i},e_{i})$) sub-gaussian noise with mean zero, and variance $\sigma^{2}$.
We denote by $Y_{i}$ the observations under this potential outcome model.

We are interested in estimating the average treatment effect (ATE)

$\tau=\frac{1}{n}\sum_{i=1}^{n}Y_{i}(1,1)-\frac{1}{n}\sum_{i=1}^{n}Y_{i}(0,0),.$
(2)

though we could, more generally, take any difference between exposure
categories. We will consider the Horvitz-Thompson estimator for a given
exposure threshold $h$:

$\hat{\tau}_{h}=\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}}{\mathbb{P}\\{z_{i}=1,e_{i}\geq
h\\}}Y_{i}-\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbf{1}\\{z_{i}=0,e_{i}\leq
1-h\\}}{\mathbb{P}\\{z_{i}=0,e_{i}\leq 1-h\\}}Y_{i}$ (3)

Our goal then is to select the threshold $h$ in the Horvitz-Thompson estimator
that minimizes the MSE. More formally, given a candidate set of thresholds
$H$, we select

$h^{*}:=\operatorname*{argmin}_{h\in
H}\text{MSE}(\hat{\tau}_{h})=\operatorname*{argmin}_{h\in
H}\text{Bias}(\hat{\tau}_{h})^{2}+\text{Var}(\hat{\tau}_{h}).$ (4)

The intuition here is that for stronger interference, i.e. when $f(e_{i})$ is
larger for a fixed $e_{i}$, the bias introduced by each edge connecting to a
neighbor with a different treatment status is higher. The MSE-optimal Horvit-
Thompson estimator in that setting is expected to be one that incorporates a
higher exposure threshold. We illustrate this through a display of a bias-
variance tradeoff from a noisy linear model $Y_{i}=\alpha+\beta z_{i}+\gamma
e_{i}+\epsilon_{i}$, that occurs across different thresholds for our adaptive
Horvitz-Thompson estimator 3 for different ratios of $\gamma/\beta$ in the
left panel of Figure 1. In the right panel of Figure 1, we illustrate how our
approach detects the bias via linear regression for the different
$\gamma/\beta$ ratios, and together with the exposure distribution, the MSE-
optimal threshold (matching the left panel) is selected, represented by the
shaded regions.

Figure 1: Left: bias, variance, and RMSE of the thresholding Horvitz-Thompson
estimator across different thresholds for the 1000-node 2nd-power cycle graph
with unit-level randomization, and noisy linear model outcome
$Y_{i}=10+10z_{i}+\gamma e_{i}+\epsilon_{i}$, for
$\epsilon_{i}\sim\mathcal{N}(0,1)$, $\gamma/\beta$ ratios of 0 (red), 0.25
(blue), 1 (orange), and 5 (green). Right: The linear fits (with colors
corresponding to the $\gamma/\beta$ ratios on the left is used to approximate
the induced exposure bias included by the thresholding. The dotted exposure
distribution (yellow) demonstrates the amount of variance reduction by
including more data through the thresholding. The rectangular blocks
illustrate the regions including the data under the MSE-optimal treshold.

We refer to this tradeoff as a general one, as due to the dependence structure
it is not-necessarily the case where the bias and variance are both
monotonically decreasing, and increasing respectively, but the general trend
is decreasing and increasing, respectively.

To illustrate how the bias and the variance of the estimator change across
different thresholds, we draw upon simple examples from Ugander et al. (2013)
in the following subsection. Figure 2 shows the different circulant graphs we
consider under unit-level Bernoulli randomization and cluster randomization.
We extend this to general graphs in the subsequent sections.

###  2.1 Tradeoffs in Circulant graphs: Cycles, and $k$th power Cycle

We first consider unit-level Bernoulli($p$) randomization in the cycle graph
and the $k$th-power cycle graphs. We say a graph is a $k$th power-Cycle graph
if there exists an edge between each node and $2k$ of its nearest neighbors
Ugander et al. (2013). We present results for general graphs subsequently.

######  Proposition 1 (Absolute bias in k-th power cycle graphs under unit-
randomization).

When $p=1/2$ and the potential outcome model is simply linear, i.e.
$Y_{i}=\alpha+\beta z_{i}+\gamma e_{i}$, the absolute bias of the Horvitz-
Thompson estimator for a given threshold $h\equiv l/2k$ for $l=0,2,...,2k$, in
the $k$th-power cycle graph, under unit-level randomization, is equal to
$\gamma\times\left[\frac{\sum_{r=l}^{2k}\left(r/k-1\right)\binom{2k}{r}}{\sum_{r=l}^{2k}\binom{2k}{r}}-1\right].$

######  Proposition 2 (Variance in k-th power cycle graphs under unit-
randomization).

When $p=1/2$ and the potential outcome model is simply linear, i.e.
$Y_{i}=\alpha+\beta z_{i}+\gamma e_{i}$, the variance of the Horvitz-Thompson
estimator for a given threshold $h$, in the $k$-th power cycle graph, under
unit-level randomization is proportional to

$\displaystyle\mathcal{O}\left(\frac{1}{np^{d}}\left[(\alpha+\beta+\gamma
dh)^{2}+(\alpha+\beta+\gamma d(1-h))^{2}-2\gamma h(1-h)d\right]\right).$

Figure 2: Circulant graphs with unit-level randomization and cluster-level
randomization. Blue and red nodes represent treated and control units,
respectively. Left: cycle graph with Ber(0.5) unit randomization. Center: 3rd-
power cycle graph with Benoulli(0.5) unit randomization. Right: 2nd-power
cycle graph with Bernoulli(0.5) cluster randomization, with clusters of size 5
(=2k + 1).

Therefore, the optimal threshold $h^{*}$ depends on $\gamma$, $\beta$, $n$,
$p$, $d$, the graph structure, and the number of other candidate thresholds.
In particular, we note the following general regimes:

  1. 1.

when $\gamma\gg\beta$, then
$\text{Bias}(\hat{\tau}_{h})^{2}\gg\text{Var}(\hat{\tau}_{h})$ for large
enough $n$, i.e. $n\gg p^{d}$. Then, $h^{*}$ would be higher (closer to one).

  2. 2.

when $\gamma\ll\beta$, then
$\text{Bias}(\hat{\tau}_{h})^{2}\ll\text{Var}(\hat{\tau}_{h})$ for not so
large $n$ with fixed $d$, i.e. $n\approx p^{d}$, and
$\text{Bias}(\hat{\tau}_{h})^{2}\gg\text{Var}(\hat{\tau}_{h})$ for larger $n$.
Therefore, we expect that $h^{*}$ decreases as the MSE transitions from the
former to the latter regime.

  3. 3.

when $\gamma\approx\beta$, then for some $n$, we are in the in-between regime,
and $h^{*}$ would be closer to $\frac{1}{2}$. For larger $n$ (while keeping
$d$ fixed), the $\text{Bias}(\hat{\tau}_{h})^{2}\gg\text{Var}(\hat{\tau}_{h})$
, and we expect to optimally compensate for this by increasing $h^{*}$ towards
1.

For $k$th-power cycle graphs, we also consider cluster-randomized design with
cluster size $2k+1$. In Ugander et al. (2013), the authors show that this
clustering size minimizes variance under full neighborhood exposure. We state
the approximate absolute bias and variance for this setting in the following
propositions. Here, define the threshold $h=\frac{l}{d}$ with $l=0,1,...,d$.

######  Proposition 3 (Absolute bias in k-th power cycle graph under cluster-
randomization).

When $p=1/2$, and the potential outcome model is simply linear, i.e.
$Y_{i}=\alpha+\beta z_{i}+\gamma e_{i}$, the bias of the Horvitz-Thompson
estimator for a given threshold $h$ in the k-th power cycle graph under
cluster randomization, with cluster-sizes $2k+1$, is approximately
$2\gamma(h-1)$ for $h\geq 1/2$.

######  Proposition 4 (Variance in the k-th power cycle graph under cluster-
randomization).

When $p=1/2$ and the potential outcome model is simply linear, i.e.
$Y_{i}=\alpha+\beta z_{i}+\gamma e_{i}$, the variance of the Horvitz-Thompson
estimator for a given threshold $h$, in the $k$-th power cycle graph under
cluster-level randomization, with cluster-sizes $2k+1$, is proportional to

$\displaystyle\frac{1}{np^{2}}(3d+1-2dh)[(\beta+\gamma
h)^{2}+(\gamma(1-h))^{2}].$

Under cluster randomization with cluster sizes $2k+1$ for the k-th power cycle
graphs, the variance grows linearly in the degrees of the graph. Therefore,
informally, compared to the unit-randomized design setting, higher node
degrees lead to stronger bias than variance for a fixed exposure function
$f(e_{i}).$

##  3 Estimating the rate of change of bias and variance

In reality, however, one does not have access to the bias and variance. With
the goal of thresholding in mind, we obtain a "bias-signal" computed via a
linear model, to estimate how the bias incurred changes by including more data
under a wider bandwidth corresponding to a smaller threshold $h$. We propose
the following absolute bias-signal estimator:

$\displaystyle\hat{b}(\hat{\tau}_{h})=\frac{1}{n}\sum_{i=1}^{n}\frac{(1-e_{i})\hat{\gamma}\mathbf{1}\\{Z_{i}=1,e_{i}\geq
h\\}}{\mathbf{P}\\{Z_{i}=1,e_{i}\geq
h\\}}+\frac{1}{n}\sum_{i=1}^{n}\frac{e_{i}\hat{\gamma}\mathbf{1}\\{Z_{i}=0,e_{i}\leq
1-h\\}}{\mathbf{P}\\{Z_{i}=0,e_{i}\leq 1-h\\}},$

where $\hat{\gamma}$ is the linear regression coefficient of the outcome on
the exposure variable.

Therefore, if the dependence of the outcome on the treatment exposure is
“strongly positive", $\hat{\gamma}$ would be larger, indicating that including
more data with a smaller threshold $h$ would induce more bias than when the
dependence on the treatment exposure was “weaker".

Next, to estimate the reduction in the variance for a drop in the threshold
$h$, we estimate the variance of the Horvitz-Thompson estimator at a given
threshold $h$ using the variance estimator discussed in Aronow and Samii
(2017). To conserve space, we write this out in Appendix A.1

Since our approach involves "double-dipping" into the data, we need to make
sure that there is no overfitting that occurs. The authors of Chernozhukov et
al. (2018), Belloni et al. (2015), and Kennedy (2022) describe this in more
detail. While sample-splitting would simply take care of this in the case of
independent data, we do not have this under the graph dependence structure we
are concerned with. We propose to use our approach on the whole sample data
and argue this via empirical process theory arguments. In particular, we can
think of our threshold $h$ as a nuisance parameter, and we show that our
estimator for $h$ has a simple enough associated MSE function class. We note
that it is sufficient to show that we are not overfitting too much by showing
that the rate of convergence of the estimated MSE under $\hat{h}$ converges to
the true MSE under the optimal threshold, under the best possible linear fit,
$h_{BL}^{*}$ at least at a $\mathcal{O}(n^{-1/2})$-rate.

######  Proposition 5.

Suppose that Conditions 1, 2, 3, and 4 are satisfied. The corresponding bias
terms in the estimated and true MSE, $\mathbb{M}_{n}^{b}$, and $M_{n}^{b}$,
under “double prediction", satisfy

$\displaystyle\mathbb{E}\left[\sup_{\delta/2<|h_{BL}^{*}-\hat{h}_{n}|<\delta}\sqrt{n}\left(\mathbb{M}_{n}^{b}(\hat{h}_{n})-M_{n}^{b}(\hat{h}_{n})-\mathbb{M}_{n}^{b}(h_{BL}^{*})-M_{n}^{b}(h_{BL}^{*})\right)\right]\to
0,$

as $\delta\to 0$.

Indeed, the terms with differences of the linear regression coefficients from
the best linear coefficient under a random design on $[0,1]$ composed with
uniformly bounded inverse propensity weights have bounded entropy integral.
That is, Donsker conditions are satisfied using Van Der Vaart et al. (1996,
Lemma 3.4.3) under Assumptions 1, 2, and 3, 4, and stochastic equicontinuity
above holds. We also have marginal Central-Limit theorem with respect to the
best linear fit holding for each of the sample-mean-like and $U$-statistics-
like terms in the MSE. Additionally, the latter in the variance parts (see
also Chen et al. (2010)) of uniformly bounded inverse propensity weights
composed with terms of bounded variation paired with sub-exponential noise
(since we have products of sub-gaussians), have bounded entropy integral Van
Der Vaart et al. (1996), which then implies a uniform Central Limit Theorem
for all of the terms. See Appendix A.5 for a more complete proof.

##  4 Theoretical and Empirical Results

We now present theoretical and simulation results in the context of general
graphs. Our exposure maps depend on respondents’ local neighborhoods, meaning
that the only restrictions we need are on the size of those neighborhoods. We
discuss an extension to more general exposure maps in the Discussion. Denote
by $H$ the set of exposure thresholds we consider. Therefore, in the cycle
graphs, $H=\\{\frac{u}{d}:u\in\\{0,1,...,d\\}\\}$ where $d$ is the degree of
the graph.

###  4.1 Theory

######  Condition 1 (Positivity).

For all $i\in[n]$, $r\in\\{0,1\\}$,

$\displaystyle\mathbb{P}(Z_{i}=r)>0$

######  Definition 1 (Bounded statistical leverage, Hsu et al. (2012)).

There exists some finite $b>1$, such that

$\displaystyle\frac{\|\Sigma^{-1/2}x\|}{\sqrt{\mathbb{E}[\|\Sigma^{-1/2}x\|^{2}]}}\equiv\frac{\|\Sigma^{-1/2}x\|}{\sqrt{2}}\leq
b.$

By the positivity assumption 1, bounded statistical leverage holds. More on
statistical leverage can be found in Mahoney et al. (2011); Martinsson and
Tropp (2020); Young (2019); Pilanci and Wainwright (2015). The idea,
informally, is that for homogeneous statistical leverage across the input
data, there are no observations that have too much control over the linear
regression.

In our setting larger spread, as well as large enough data near the boundaries
$\\{0,1\\}$ are important for the robustness to non-linear models. In
particular, having more high-leverage data (at the boundaries) allows for our
linear fit to be more informative in estimating the bias of the average
treatment effect, even in the presence of the non-linearity. That is, our
approach is best when we are able to approximate the slope at the boundaries
well, regardless of the true potential outcome model.

######  Condition 2 (Bounded variation treatment-exposure function).

Let the potential outcome function be
$y_{i}=\alpha+\psi(z_{i},e_{i})+\epsilon_{i}$ as in 1. The function
$\psi(z_{i},e_{i})$ has bounded variation.

######  Definition 2 (Bounded approximation error, Hsu et al. (2012)).

Denote $\left[\beta_{BL},\gamma_{BL}\right]$ by $\theta$, and
$\left[z,e\right]$ by $x$. There exists finite $b_{\text{approx}}\geq 0$ such
that, almost surely,

$\frac{\|\Sigma^{-1/2}x(\psi(x)-\langle\theta,x\rangle)\|}{\sqrt{\mathbb{E}[\|\Sigma^{-1/2}x\|^{2}]}}=\frac{\|\Sigma^{-1/2}x(\psi(x)-\langle\theta,x\rangle)\|}{\sqrt{2}}\leq
b_{approx}.$

That is, together with Condition 1, we assume that there exists finite
$b_{approx}\geq 0$, such that $|(f(x)-\langle\theta,x\rangle)|\leq
b_{approx}/b$. The following then holds:

$\displaystyle\frac{1}{\sqrt{2}}\|\Sigma^{-1/2}x(\psi(x)-\langle\theta,x\rangle)\|$
$\displaystyle\leq\frac{1}{\sqrt{2}}\|\Sigma^{-1/2}x\||(\psi(x)-\langle\theta,x\rangle)|$
$\displaystyle\leq b|(\psi(x)-\langle\theta,x\rangle)|\leq b_{approx}.$

######  Condition 3 (Subgaussian noise).

There exists $0\leq\sigma^{2}<\infty$ such that for all $i\in[n],$

$\displaystyle\mathbb{E}[\exp(\lambda(\epsilon_{i}))]\leq\exp(\lambda^{2}\sigma^{2}/2)$

almost surely for all $\lambda\in\mathbb{R}$.

######  Condition 4 (Bounded first-order interactions).

Let $B_{r}(i)$ denote the ball of radius $r$ around unit $i$. For all
$i\in[n]$,

$|B_{1}(i)|=O(1).$

The condition above is important under unit-randomization design. For a
cluster-randomization design, we would require bounded interactions across
clusters. In particular the graph interference structure should have small
$k_{n}$-conductance (i.e. the $k_{n}$-partition generalization of the Cheeger
constant) for $k_{n}$ growing with $n$. Though restrictive theoretically, it
is consistent with the observation that networks in practice tend to be sparse
and clustered (Chandrasekhar et al., 2020). Violating this condition will
violate exposure positivity (see remark below), and Proposition 5 will not
hold.

######  Remark.

In the weighted graph setting, if practitioners consider equivalence classes
$\tilde{e}_{i}$ of exposures, so that $\tilde{e}_{i}\equiv\tilde{e}_{j}$ if
and only if $|e_{i}-e_{j}|<\varepsilon_{n}$ for $\epsilon_{n}$ specified
below, exposure positivity, i.e. for all $i\in[n]$, $r\in\\{0,1\\}$, and
$s\in[0,1]$, $\mathbf{P}(Z_{i}=r,\tilde{e}_{i}=s)>0$ is satisfied. Then, we
can weaken our bounded degree constraint. Let $B_{r}(i)$ denote the ball of
radius $r$ around unit $i$, such that it decomposes,
$|B_{1}(i)|=|B_{1}^{S}(i)|+|B_{1}^{W}(i)|$ for all $i\in[n]$. $B_{1}^{S}(i)$
is the ball of radius 1 around unit $i$ with strong connections $\eta_{is}$,
and $|B_{1}^{W}(i)|$ be the ball of radius 1 around unit $i$ with weak
connections $\delta_{iw}$ such that
$\sup_{i,s,w}\frac{\eta_{is}}{\delta_{iw}}\to\infty$. For all $i\in[n]$, we
require that $|B_{1}^{S}(i)|=\mathcal{O}(1)$. With an abuse of notation, write
$e_{i}(\mathcal{N})$ to mean the exposure of unit $i$ as a function of the
treated units in the subset $\mathcal{N}\subseteq B_{1}(i)$. Then, we require
also also that $e_{i}(B_{1}^{S}(i))=\Omega(e_{i}(B_{1}(i)))$. Additionally,
defining $v_{i}=z_{i}-\mathbf{E}[z_{i}]$ for random variables $z_{i}$,
$|B_{1}(i)|$ must satisfy the three conditions in Chandrasekhar et al. (2023)
generalizing Aronow and Samii (2017), with an additional condition that the
sum of the affinity set covariances is $\Omega(n)$, for the bias terms as well
as the variance terms in the MSE. Since those conditions nest strong mixing
dependence structure, which nests our current results for the variance terms,
these conditions are sufficient. Then, taking
$\varepsilon_{n}:=\varepsilon(\delta_{n},B_{1}^{W}(i))>0$, gives us uniform
control on the estimated MSE of the estimator with respect to the true MSE of
the estimator, and all the results follow through.

The circulant graphs we consider are the ring graphs that maintain out-of-
synchronization in Wiley et al. (2006). This points to the sparsity of the
circulant graphs. Some examples of graph classes that fall under the
categorization of Condition 4 are expander graphs, and growth-restricted
graphs Alon (1986); Arora et al. (2009); Gkantsidis et al. (2003); Kuhn et al.
(2005); Krauthgamer and Lee (2003); Kowalski (2019). In the causal inference
literature, the effect of the growth-restricted graph setting was studied in
Ugander et al. (2013). For expository purposes, in Section 2.1, we focus our
analysis and simulations on a special case of circulant graphs in growth-
restricted graphs, particularly the cycle graph, and the $k$-th power Cycle
graphs, just as in Ugander et al. (2013).

Define $M_{n}(h),M_{n}(h_{BL}^{*}),M_{n}(h^{*})$ to be the MSE, evaluated at
an arbitrary threshold $h$, at the optimal threshold under the best possible
linear fit, and at the optimal threshold generally, respectively. Then, we
have,

$\displaystyle|M_{n}(h)-M_{n}(h^{*})|$
$\displaystyle\leq|M_{n}(h)-M_{n}(h_{BL}^{*})|+|M_{n}(h_{BL}^{*})-M_{n}(h^{*})|$
(5) $\displaystyle\leq|M_{n}(h)-M_{n}(h_{BL}^{*})|+\text{approx},$ (6)

where

$\displaystyle\text{approx}=|\frac{1}{n}\sum_{i=1}^{n}\sum_{x_{i}\in X_{i}}$
$\displaystyle\left\\{\left(\mathbf{1}\\{x\geq
h_{BL}^{*}\\}\gamma_{BL}(1-x)-\mathbf{1}\\{x\geq
h^{*}\\}(f(1)-f(x))\right)\right.$
$\displaystyle+\left.\left(\mathbf{1}\\{x\leq
1-h_{BL}^{*}\\}\gamma_{BL}x-\mathbf{1}\\{x\leq
1-h^{*}\\}(f(x)-f(0))\right)\right\\}|,$

where $x_{i}$ ranges over the set $X_{i}$ of possible fractions of unit $i$’s
degree.

How well the linear model approximates the outcome-dependence on the exposure
partly describes this quantity. If the potential outcome model does indeed
depend linearly on the exposure, this error is identically zero. Additionally,
the more spread-out the exposure distribution towards the boundaries
$\\{0,1\\}$, the better the approximation across $[0,1-h]$ and $[h,1]$ both of
which characterizes $|h_{BL}^{*}-h^{*}|$, and hence the smaller the
approximation error. Graph structures that are more sparse for unit-level
randomization, or graph structures with inherent strong clustering such that
the edges have high intra-cluster density and much lower inter-cluster
density, with growing number of clusters for cluster-randomization, allow for
more spread-out exposure distributions. That is, in the latter, under cluster-
randomization for graphs with small $k_{n}$-conductance for $k_{n}$ growing
with $n$, we have a more spread-out exposure distribution, as exposure
concentration occurs towards more than one centroid depending on the graph
clusters formed, with some centroids being at the boundaries as a result of
cluster-randomization. See also the note after Condition 4.

In the following theorem, we characterize the probability of choosing the
correct threshold under the best linear fit (i.e. the probability the first
term on the right-hand-side in 5 is zero).

######  Theorem 1.

Suppose that $\mathbb{M}_{n}(h)$ is the estimated MSE, and $M_{n}(h)$ is the
true MSE, associated with a given graph of $n$ nodes and experimental design.
Assume that Conditions 1, 2, 3, and 4 hold. Then, define
$K:=\mathbb{E}[\|\Sigma^{-1/2}x\text{approx}(x)\|^{2}]$,
$\hat{F}:=\left(\frac{1}{n}\sum_{i=1}^{n}\left[\frac{\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}}{\mathbb{P}\\{z_{i}=1,e_{i}\geq
h\\}}e_{i}+\frac{\mathbf{1}\\{z_{i}=0,e_{i}\leq
1-h\\}}{\mathbb{P}\\{z_{i}=0,e_{i}\leq 1-h\\}}(1-e_{i})\right]\right)^{2}$,
$\Delta:=\min_{h,h^{\prime}:h\neq h^{\prime}}|M_{n}(h)-M_{n}(h^{\prime})|$,
$R:=\max_{i}d_{i}$, and $B$ to be the maximum over $i\in[n]$ of the summands
in $\mathbb{M}_{n}(h)$ and $M_{n}(h)$, with $Y_{i}-\epsilon_{i}$ in the
variance terms. $B<\infty$ by Conditions 1,2, and 4. Let $H$ be the set of
candidate thresholds. By Condition 4, $|H|\leq S$ for some constant
$S<\infty$. Let $\gamma_{0}$ be the best possible linear fit. Then,

$\displaystyle\mathbb{P}\\{|M_{n}(h)=M_{n}(h_{BL}^{*})|\\}\geq
1-3SC\exp{\left(-\frac{n}{C_{1}(16K^{2}+2\sigma^{2})}\cdot\left(\frac{\phi^{2}\Delta}{12\gamma_{0}\hat{F}}-C_{2}\right)\right)},$

for some constants $C,C_{1},C_{2}$, and variance $\phi^{2}$ of $e_{i}$.

This gives us the following corollary.

######  Corollary 1.

For the $k$th-power cycle graph of $n$ nodes, under a noisy linear outcome
model, i.e., where $g_{0}(z_{i})=\beta z_{i}$, and
$f_{0}(e_{i})=\gamma_{0}e_{i}$, and $\epsilon_{i}\sim\mathcal{N}(0,1)$ for all
$i\in[n]$,

$\displaystyle\mathbb{P}\\{|M_{n}(h)=M_{n}(h^{*})|\\}\geq
1-3(2k+1)C\exp{\left(-\frac{n}{C_{1}(16K^{2}+2\sigma^{2})}\cdot\left(\frac{\phi^{2}\Delta}{12\gamma_{0}\hat{F}}-C_{2}\right)\right)},$

for some constants $C,C_{1},C_{2}$, and variance $\phi^{2}$ of $e_{i}$.

###  4.2 Simulations

We now compare the performance of our approach to that of the Horvitz-Thompson
estimator with threshold one, and the Horvitz-Thompson estimator with
threshold zero. Additionally, since our objective is essentially an MSE-
optimal bandwidth selection problem for the indicator kernel (see Appendix
A.3, 13), we compare our estimator to the Horvitz-Thompson estimator with the
plugin threshold from Lepski’s method Goldenshluger and Lepski (2011); Su et
al. (2020). We note that Lepski’s method requires monontonicity of the bias
and the variance to work well. In our setting, due to the implicit dependency
structure from the graph in the inverse-propensity weights, the monotonicity
assumption may be violated (see left panel of Figure 1). We write the three
estimators out in Appendix A.2.

In Figure 3, under the linear model with $\alpha=10$, $g(z_{i})=\beta
z_{i}=10z_{i}$, $f(e_{i})=\gamma e_{i}$, and
$\epsilon_{i}\sim\mathcal{N}(0,1),$ we see that our adaptive thresholding
estimator (AdaThresh) consistently performs better than existing estimators,
“interpolating" between the fixed 0/1-threshold Horvitz-Thompson estimators,
and out-performing the Lepski’s method plug-in estimators when the threshold
$\gamma/\beta$ is high.

Figure 3: RMSE (normalized by the ATE) induced by the different Horvitz-
Thompson estimators. Left: Cycle graph under unit-level Ber(0.5)
randomization. Center: 2nd-power cycle graph under unit-level Ber(0.5)
randomization. Right: 2nd-power cycle graph under cluster-level Ber(0.5)
randomization with cluster sizes 5 (=2k + 1). The error bars are two times the
standard deviation.

###  4.3 Real Data

We evaluate the performance of our estimator on the Amazon (DVD) products
similarity network Leskovec et al. (2007). We consider the performance of the
estimator on a subset of (the first) 1000 non-isolated nodes from a total of
19828 nodes, with exposure probabilities computed on the full graph. We
generate simulated data using the linear model with $\alpha=10,g(z_{i})=\beta
z_{i}=10z_{i}$, $f(e_{i})=\gamma e_{i}$, and
$\epsilon_{i}\sim\mathcal{N}(0,1)$ under a unit-level Bernoulli(0.5)
randomization design. Figure 4 shows that we consistently perform better than
existing fixed estimators. See Appendix A.7 for more details.

Figure 4: RMSE (normalized by the ATE) across different thresholds on the
Amazon product network data. The error bars are two times the standard
deviation.

##  5 Discussion

In this paper, we investigated the adaptive Horvitz-Thopson estimator for
symmetric thresholds $h$ and $1-h$ for treatment and control, respectively.
One could also use this approach for a more general estimator with thresholds
$h_{1}$ and $1-h_{0}$, respectively with $h_{1}\neq h_{0}$. Additionally, in
Appendix A.7.4, we demonstrate the performance of this approach using the
Difference-in-Means estimator incorporating exposure thresholds. In Appendix
A.7.3, we illustrate the robustness of our estimator to deviations from
linearity in the potential outcomes model. We note that our framework would
apply in more general settings, including a direct extension to friends-of-
friends (or further connections) or to neighborhoods that are either observed
or learned through a clustering algorithm. An additional interesting future
direction could be to explore an extension of our work to the framework
presented in Chandrasekhar et al. (2023), which allows each individual to have
a unique exposure and then bins exposures that are separated by at least a
threshold.

Furthermore, to improve upon robustness to non-linear settings, we propose an
extension using local regression to estimate the rate of change of bias within
the $[0,1-h],[h,1]$ windows. As long as there is sufficient concentration of
exposures in these windows, robustness is never worse. In Appendix A.7.5, we
display our simulation results in this setting using local linear regression
to estimate the rate of change of bias within the $[0,1-h],[h,1]$ windows. One
could also use other local regression approaches, such as kernel regression,
etc., while maintaining sufficiently small model complexity to avoid needing
to sample-split.

## Acknowledgments and Disclosure of Funding

The authors would like to thank Alex Kokot, and Lars van der Laan for helpful
discussions; and Matthew Eichhorn for suggesting the Amazon product similarity
network dataset.

## References

  * Centola and Macy [2007] Damon Centola and Michael Macy.  Complex Contagions and the Weakenss of Long Ties.  _American Journal of Sociology_ , 113(3):702–734, 2007.
  * Acemoglu et al. [2011] Daron Acemoglu, Asuman Ozdaglar, and Ercan Yildiz.  Diffusion of innovations in social networks.  In _2011 50th IEEE conference on decision and control and European control conference_ , pages 2329–2334. IEEE, 2011.
  * Reich [2016] Bryony Reich.  The diffusion of innovations in social networks.  _Working paper, University College London_ , 2016.
  * Aronow and Samii [2017] Peter M Aronow and Cyrus Samii.  Estimating average causal effects under general interference, with application to a social network experiment.  2017\.
  * Ugander et al. [2013] Johan Ugander, Brian Karrer, Lars Backstrom, and Jon Kleinberg.  Graph cluster randomization: Network exposure to multiple universes.  In _Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 329–337, 2013.
  * Sussman and Airoldi [2017] Daniel L Sussman and Edoardo M Airoldi.  Elements of estimation theory for causal effects in the presence of network interference.  _arXiv preprint arXiv:1702.03578_ , 2017.
  * Hardy et al. [2019] Morgan Hardy, Rachel M Heath, Wesley Lee, and Tyler H McCormick.  Estimating spillovers using imprecisely measured networks.  _arXiv preprint arXiv:1904.00136_ , 2019.
  * Auerbach and Tabord-Meehan [2021] Eric Auerbach and Max Tabord-Meehan.  The local approach to causal inference under network interference.  _arXiv preprint arXiv:2105.03810_ , 2021.
  * Eckles et al. [2017] Dean Eckles, Brian Karrer, and Johan Ugander.  Design and analysis of experiments in networks: Reducing bias from interference.  _Journal of Causal Inference_ , 5(1):20150021, 2017.
  * Toulis and Kao [2013] Panos Toulis and Edward Kao.  Estimation of causal peer influence effects.  In _International conference on machine learning_ , pages 1489–1497. PMLR, 2013.
  * Basse and Airoldi [2018] Guillaume W Basse and Edoardo M Airoldi.  Model-assisted design of experiments in the presence of network-correlated outcomes.  _Biometrika_ , 105(4):849–858, 2018.
  * Cai et al. [2015] Jing Cai, Alain De Janvry, and Elisabeth Sadoulet.  Social networks and the decision to insure.  _American Economic Journal: Applied Economics_ , 7(2):81–108, 2015.
  * Zhu et al. [2024] Hongtao Zhu, Sizhe Zhang, Yang Su, Zhenyu Zhao, and Nan Chen.  Integrating active learning in causal inference with interference: A novel approach in online experiments.  _arXiv preprint arXiv:2402.12710_ , 2024.
  * Su et al. [2020] Yi Su, Pavithra Srinath, and Akshay Krishnamurthy.  Adaptive estimator selection for off-policy evaluation.  In _International Conference on Machine Learning_ , pages 9196–9205. PMLR, 2020.
  * Fan and Gijbels [1992] Jianqing Fan and Irene Gijbels.  Variable bandwidth and local linear regression smoothers.  _The Annals of Statistics_ , pages 2008–2036, 1992.
  * Ruppert [1997] David Ruppert.  Empirical-bias bandwidths for local polynomial nonparametric regression and density estimation.  _Journal of the American Statistical Association_ , 92(439):1049–1062, 1997.
  * Kallus and Zhou [2018] Nathan Kallus and Angela Zhou.  Policy evaluation and optimization with continuous treatments.  In _International conference on artificial intelligence and statistics_ , pages 1243–1251. PMLR, 2018.
  * Goldenshluger and Lepski [2011] Alexander Goldenshluger and Oleg Lepski.  Bandwidth selection in kernel density estimation: oracle inequalities and adaptive minimax optimality.  2011\.
  * Belloni et al. [2022] Alexandre Belloni, Fei Fang, and Alexander Volfovsky.  Neighborhood adaptive estimators for causal inference under network interference.  _arXiv preprint arXiv:2212.03683_ , 2022.
  * Chernozhukov et al. [2018] Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins.  Double/debiased machine learning for treatment and structural parameters.  _The Econometrics Journal_ , 21(1):C1–C68, 01 2018.  ISSN 1368-4221.  doi: 10.1111/ectj.12097.  URL https://doi.org/10.1111/ectj.12097.
  * Belloni et al. [2015] Alexandre Belloni, Victor Chernozhukov, Ivan Fernández-Val, and Christian Hansen.  Program evaluation with high-dimensional data.  Technical report, cemmap working paper, 2015.
  * Kennedy [2022] Edward H Kennedy.  Semiparametric doubly robust targeted double machine learning: a review.  _arXiv preprint arXiv:2203.06469_ , 2022.
  * Van Der Vaart et al. [1996] Aad W Van Der Vaart, Jon A Wellner, Aad W van der Vaart, and Jon A Wellner.  _Weak convergence_.  Springer, 1996.
  * Chen et al. [2010] Louis HY Chen, Larry Goldstein, and Qi-Man Shao.  _Normal approximation by Stein’s method_.  Springer Science & Business Media, 2010.
  * Hsu et al. [2012] Daniel Hsu, Sham M Kakade, and Tong Zhang.  Random design analysis of ridge regression.  In _Conference on learning theory_ , pages 9–1. JMLR Workshop and Conference Proceedings, 2012.
  * Mahoney et al. [2011] Michael W Mahoney et al.  Randomized algorithms for matrices and data.  _Foundations and Trends® in Machine Learning_ , 3(2):123–224, 2011.
  * Martinsson and Tropp [2020] Per-Gunnar Martinsson and Joel A Tropp.  Randomized numerical linear algebra: Foundations and algorithms.  _Acta Numerica_ , 29:403–572, 2020.
  * Young [2019] Alwyn Young.  Channeling fisher: Randomization tests and the statistical insignificance of seemingly significant experimental results.  _The quarterly journal of economics_ , 134(2):557–598, 2019.
  * Pilanci and Wainwright [2015] Mert Pilanci and Martin J Wainwright.  Randomized sketches of convex programs with sharp guarantees.  _IEEE Transactions on Information Theory_ , 61(9):5096–5115, 2015.
  * Chandrasekhar et al. [2020] Arun G Chandrasekhar, Horacio Larreguy, and Juan Pablo Xandri.  Testing models of social learning on networks: Evidence from two experiments.  _Econometrica_ , 88(1):1–32, 2020.
  * Chandrasekhar et al. [2023] Arun G Chandrasekhar, Matthew O Jackson, Tyler H McCormick, and Vydhourie Thiyageswaran.  General covariance-based conditions for central limit theorems with dependent triangular arrays.  _arXiv preprint arXiv:2308.12506_ , 2023.
  * Wiley et al. [2006] Daniel A Wiley, Steven H Strogatz, and Michelle Girvan.  The size of the sync basin.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 16(1), 2006.
  * Alon [1986] Noga Alon.  Eigenvalues, geometric expanders, sorting in rounds, and ramsey theory.  _Combinatorica_ , 6(3):207–219, 1986.
  * Arora et al. [2009] Sanjeev Arora, Satish Rao, and Umesh Vazirani.  Expander flows, geometric embeddings and graph partitioning.  _Journal of the ACM (JACM)_ , 56(2):1–37, 2009.
  * Gkantsidis et al. [2003] Christos Gkantsidis, Milena Mihail, and Amin Saberi.  Conductance and congestion in power law graphs.  In _Proceedings of the 2003 ACM SIGMETRICS International Conference on Measurement and modeling of computer systems_ , pages 148–159, 2003.
  * Kuhn et al. [2005] Fabian Kuhn, Thomas Moscibroda, and Rogert Wattenhofer.  On the locality of bounded growth.  In _Proceedings of the twenty-fourth annual ACM symposium on Principles of distributed computing_ , pages 60–68, 2005.
  * Krauthgamer and Lee [2003] Robert Krauthgamer and James R Lee.  The intrinsic dimensionality of graphs.  In _Proceedings of the thirty-fifth annual ACM symposium on Theory of computing_ , pages 438–447, 2003.
  * Kowalski [2019] Emmanuel Kowalski.  _An introduction to expander graphs_.  Société mathématique de France Paris, 2019.
  * Leskovec et al. [2007] Jure Leskovec, Lada A Adamic, and Bernardo A Huberman.  The dynamics of viral marketing.  _ACM Transactions on the Web (TWEB)_ , 1(1):5–es, 2007.
  * Van der Vaart [2000] Aad W Van der Vaart.  _Asymptotic statistics_ , volume 3.  Cambridge university press, 2000.
  * Sancetta [2009] Alessio Sancetta.  Maximal inequalities for u-processes of strongly mixing random variables.  _Probability and Mathematical Statistics_ , 29, 2009.
  * Beutner and Zähle [2012] Eric Beutner and Henryk Zähle.  Deriving the asymptotic distribution of u-and v-statistics of dependent data using weighted empirical processes.  2012\.
  * Yoshihara [1976] Ken-ichi Yoshihara.  Limiting behavior of u-statistics for stationary, absolutely regular processes.  _Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete_ , 35(3):237–252, 1976.
  * Arcones and Yu [1994] Miguel Angel Arcones and Bin Yu.  Central limit theorems for empirical and u-processes of stationary mixing sequences.  _Journal of Theoretical Probability_ , 7(1):47–71, 1994.
  * Denker and Keller [1986] Manfred Denker and Gerhard Keller.  Rigorous statistical procedures for data from dynamical systems.  _Journal of Statistical Physics_ , 44:67–93, 1986.
  * Chatterjee [2005] Sourav Chatterjee.  _Concentration inequalities with exchangeable pairs_.  Stanford University, 2005.

##  Appendix A Appendix / supplemental material

###  A.1 Variance estimator

The variance estimator, for the Horvitz-Thompson estimator, we use is the one
proposed by Aronow and Samii [2017, Eqn. 7,9].

$\begin{split}\hat{Var}(\hat{\tau}_{h})&=\sum_{i=1}^{n}\frac{\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}Y_{i}^{2}}{n^{2}\pi_{i}^{1}}(\frac{1}{\pi_{i}^{1}}-1)\\\
&+\sum_{i=1}^{n}\sum_{\mathclap{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}}^{n}\frac{\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}\mathbf{1}\\{z_{j}=1,e_{j}\geq
h\\}Y_{i}Y_{j}}{n^{2}\pi_{ij}^{11}}(\frac{\pi_{ij}^{11}}{\pi_{i}^{1}\pi_{j}^{1}}-1)\\\
&+\sum_{i=1}^{n}\frac{\mathbf{1}\\{z_{i}=0,e_{i}\leq
1-h\\}Y_{i}^{2}}{n^{2}\pi_{i}^{0}}(\frac{1}{\pi_{i}^{0}}-1)\\\
&+\sum_{i=1}^{n}\sum_{\mathclap{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}}^{n}\frac{\mathbf{1}\\{z_{i}=0,e_{i}\leq
1-h\\}\mathbf{1}\\{z_{j}=0,e_{j}\leq
1-h\\}Y_{i}Y_{j}}{n^{2}\pi_{ij}^{00}}(\frac{\pi_{ij}^{00}}{\pi_{i}^{0}\pi_{j}^{0}}-1)\\\
&-\frac{2}{n^{2}}\sum_{i=1}^{n}\sum_{j\in[n];\pi_{ij}^{10}>0}(\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}\mathbf{1}\\{z_{j}=0,e_{j}\leq
1-h\\}Y_{i}Y_{j})(\frac{1}{\pi_{i}^{1}\pi_{j}^{0}}-\frac{1}{\pi_{ij}^{10}})\\\
&+\frac{2}{n^{2}}\sum_{i=1}^{n}\sum_{j\in[n];\pi_{ij}^{10}=0}(\frac{\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}Y_{i}^{2}}{2\pi_{i}^{1}}+\frac{\mathbf{1}\\{z_{j}=0,e_{j}\leq
1-h\\}Y_{j}^{2}}{2\pi_{j}^{0}})\end{split}$ (7)

Here, $\pi_{i}^{x},\pi_{ij}^{xy}$ are defined with respect to the threshold
$h$. That is, $\pi_{i}^{1}=\mathbb{P}\\{z_{i}=1,e_{i}\geq h\\}$,
$\pi_{i}^{0}=\mathbb{P}\\{z_{i}=0,e_{i}\leq 1-h\\}$,
$\pi_{ij}^{10}=\mathbb{P}\\{z_{i}=1,e_{i}\geq h,z_{j}=0,e_{i}\leq 1-h\\}$,
$\pi_{ij}^{11}=\mathbb{P}\\{z_{i}=1,e_{i}\geq h,z_{j}=1,e_{i}\geq h\\}$,
$\pi_{ij}^{01}=\mathbb{P}\\{z_{i}=0,e_{i}\leq 1-h,z_{j}=1,e_{i}\geq h\\}$,
$\pi_{ij}^{00}=\mathbb{P}\\{z_{i}=0,e_{i}\leq 1-h,z_{j}=0,e_{i}\leq 1-h\\}$.

###  A.2 Horvitz-Thompson estimators with existing approaches

We write the vanilla Horvitz-Thompson estimator (at threshold one) 8, the
Horvitz-Thompson estimator at threshold zero 9, and the Horvitz-Thompson
estimator with the threshold selected via Lepski’s method 12 below.

####  A.2.1 Vanilla Horvitz-Thompson estimator (at threshold one)

$\displaystyle\hat{\tau}_{\text{HT}_{1}}=\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbf{P}\\{Z_{i}=1,e_{i}=1\\}}{\mathbf{1}\\{Z_{i}=1,e_{i}=1\\}}Y_{i}-\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbf{1}\\{Z_{i}=0,e_{i}=0\\}}{\mathbf{P}\\{Z_{i}=0,e_{i}=0\\}}Y_{i}$
(8)

####  A.2.2 Horvitz-Thompson estimator at threshold zero

$\displaystyle\hat{\tau}_{\text{HT}_{0}}=\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbf{P}\\{Z_{i}=1\\}}{\mathbf{1}\\{Z_{i}=1\\}}Y_{i}-\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbf{1}\\{Z_{i}=0\\}}{\mathbf{P}\\{Z_{i}=0\\}}Y_{i}$
(9)

####  A.2.3 Lepski’s method

As described in Su et al. [2020], we first take

$\displaystyle
I(h):=[\hat{\tau}_{\text{HT}_{h}}-2\hat{\text{SDEV}}(\hat{\tau}_{\text{HT}_{h}}),\hat{\tau}_{\text{HT}_{h}}+2\hat{\text{SDEV}}(\hat{\tau}_{\text{HT}_{h}})].$
(10)

Then, take

$\displaystyle\hat{h}_{\text{Lepski}}:=\min\\{h\in
H:\cap_{h=1}^{k}I(h)\neq\emptyset\\},$ (11)

and

$\displaystyle\hat{\tau}_{\text{LepskiHT}}=\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbf{1}\\{Z_{i}=1,e_{i}\geq\hat{h}_{\text{Lepski}}\\}}{\mathbf{P}\\{Z_{i}=1,e_{i}\geq\hat{h}_{\text{Lepski}}\\}}Y_{i}-\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbf{1}\\{Z_{i}=0,e_{i}\leq
1-\hat{h}_{\text{Lepski}}\\}}{\mathbf{P}\\{Z_{i}=0,e_{i}\leq
1-\hat{h}_{\text{Lepski}}\\}}Y_{i}$ (12)

###  A.3 An equivalent formulation of the Horvitz-Thompson estimator for
exposure threshold $h$

We can rewrite the display in 3 in the following form:

$\hat{\tau}_{h}=\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbf{1}\\{|z_{i}-e_{i}|\leq
h\\}}{\mathbb{P}\\{|z_{i}-e_{i}|\leq h\\}}\tilde{z}_{i}Y_{i}$ (13)

where $\tilde{z}_{i}=2z_{i}-1$, so that $\tilde{z}_{i}\in\\{-1,1\\}.$ It is
then clear that the threshold $h$ controls how much dissimilarity between the
treatment status of units and their neighbors, we allow in our estimation.
This allows us to reframe the problem as an optimal bandwidth selection one.

###  A.4 Proofs to Propositions 1, 2, 3, 4

###### Proof to Proposition 1.

When $p=1/2$ and the potential outcome model is simply linear, i.e.
$Y_{i}=\alpha+\beta z_{i}+\gamma e_{i}$, the absolute bias of the Horvitz-
Thompson estimator for a given threshold $h\equiv l/2k$ for $l=0,2,...,2k$, in
the $k$th-power cycle graph under unit-randomization is equal to is equal to

$\displaystyle\frac{1}{\mathbb{P}\\{z_{i}=1,e_{i}\geq h\\}}\sum_{x_{i}\in
X}\mathbf{1}\\{x_{i}\geq
h\\}\mathbb{E}[\mathbf{1}\\{z_{i}=1,e_{i}=x_{i}\\}]y_{i}(x_{i})$
$\displaystyle-\frac{1}{\mathbb{P}\\{z_{i}=0,e_{i}\leq 1-h\\}}\sum_{x_{i}\in
X}\mathbf{1}\\{x_{i}\leq
1-h\\}\mathbb{E}[\mathbf{1}\\{z_{i}=0,e_{i}=x_{i}\\}]y_{i}(x_{i})$
$\displaystyle-\left(\beta+\gamma\right)$
$\displaystyle=\frac{1}{\sum_{r=l}^{2k}\binom{2k}{r}p^{2k}}\sum_{x_{i}\in
X}\mathbf{1}\\{x_{i}\geq
h\\}\mathbb{E}[\mathbf{1}\\{z_{i}=1,e_{i}=x_{i}\\}]y_{i}(x_{i})$
$\displaystyle-\frac{1}{\sum_{r=l}^{2k}\binom{2k}{r}p^{2k}}\sum_{x_{i}\in
X}\mathbf{1}\\{x_{i}\leq
1-h\\}\mathbb{E}[\mathbf{1}\\{z_{i}=0,e_{i}=x_{i}\\}]y_{i}(x_{i})$
$\displaystyle-(\beta+\gamma)$
$\displaystyle=\frac{p^{2k}}{\sum_{r=l}^{2k}\binom{2k}{r}p^{2k}}\left(\sum_{r=l}^{2k}\binom{2k}{r}\beta+\binom{2k}{r}(r/k-1)\gamma\right)-(\beta+\gamma)$
$\displaystyle=\frac{1}{\sum_{r=l}^{2k}\binom{2k}{r}p^{2k}}\left(\sum_{r=l}^{2k}\binom{2k}{r}(r/k-1)\gamma\right)-\gamma$

where $X=\\{0,1/2k,2/2k,...,1\\}.$ ∎

###### Proof to Proposition 2.

When $p=1/2$ and the potential outcome model is simply linear, i.e.
$Y_{i}=\alpha+\beta z_{i}+\gamma e_{i}$, it is not difficult to see that the
variance of the Horvitz-Thompson estimator for a given threshold $h\equiv
l/2k$, in the $k$-th power cycle graph under unit-level randomization is
proportional to

$\displaystyle\frac{1}{n}\sum_{l=0}^{2k}\binom{2k}{l}\mathbf{1}\\{l/2k\geq
h\\}\left[(\alpha+\beta+\gamma
l/2k)^{2}(\frac{1}{\sum_{l=0}^{2k}\binom{2k}{l}\mathbf{1}\\{l/2k\geq
h\\}p^{2k}}-1)\right.$ $\displaystyle+\left.2(\beta+\gamma
l/(2k))(\gamma(1-l/2k))(2k+1-\frac{2k}{\sum_{l=0}^{2k}\binom{2k}{l}\mathbf{1}\\{l/2k\geq
h\\}p^{2k}})\right.$
$\displaystyle+\left.(\alpha+\gamma(1-l/2k))^{2}(\frac{1}{\sum_{l=0}^{2k}\binom{2k}{l}\mathbf{1}\\{l/2k\geq
h\\}p^{2k}}-1)\right].$

Considering the dominating terms gives us the result. ∎

###### Proof to Proposition 3.

When $p=1/2$ and the potential outcome model is simply linear, i.e.
$Y_{i}=\alpha+\beta z_{i}+\gamma e_{i}$, the absolute bias of the Horvitz-
Thompson estimator for a given threshold $h\equiv l/2k$ for $l=k,k+1,...,2k$,
in the $k$th-power cycle graph under cluster-randomization is equal to is
equal to

$\displaystyle\frac{1}{\mathbb{P}\\{z_{i}=1,e_{i}\geq h\\}}\sum_{x_{i}\in
X}\mathbf{1}\\{x_{i}\geq
h\\}\mathbb{E}[\mathbf{1}\\{z_{i}=1,e_{i}=x_{i}\\}]y_{i}(x_{i})$
$\displaystyle-\frac{1}{\mathbb{P}\\{z_{i}=0,e_{i}\leq 1-h\\}}\sum_{x_{i}\in
X}\mathbf{1}\\{x_{i}\leq
1-h\\}\mathbb{E}[\mathbf{1}\\{z_{i}=0,e_{i}=x_{i}\\}]y_{i}(x_{i})$
$\displaystyle-\left(\beta+\gamma\right)$
$\displaystyle=\frac{(p/(2k+1)+2k/(2k+1)\cdot
p^{2})(\beta+\gamma)+2p/(2k+1)\sum_{r=1}^{k}\mathbf{1}\\{r\geq
d-l\\}(\beta+\gamma\cdot(2l-d)/d)}{p/(2k+1)+p^{2}\cdot
2k/(2k+1)+2p/(2k+1)\sum_{r=1}^{k}\mathbf{1}\\{r\geq d-l\\}}$
$\displaystyle-(\beta+\gamma)$ $\displaystyle=\mathcal{O}(2\gamma(h-2))$

where $X=\\{0,1/2k,2/2k,...,1\\}.$ When $l<k$, the bias scales is the same as
when $l=k.$ ∎

###### Proof to Proposition 4.

When $p=1/2$ and the potential outcome model is simply linear, i.e.
$Y_{i}=\alpha+\beta z_{i}+\gamma e_{i}$, it is not difficult to see that the
variance of the Horvitz-Thompson estimator for a given threshold $h\equiv
l/2k$, in the $k$-th power cycle graph under cluster-randomization, with
cluster-size $2k+1$, is proportional to

$\displaystyle\sum_{u=0}^{d}\mathbf{1}\\{u/d\geq
h\\}\binom{d}{u}\sum_{i=1}^{n}(\frac{(\beta+\gamma
u/d)^{2}}{n^{2}})[(3d+1-2l)(\frac{1}{p^{2}}-1)+(2l+1)(\frac{1}{p}-1)]$
$\displaystyle+\sum_{u=0}^{d}\mathbf{1}\\{u/d\geq
h\\}\binom{d}{u}\sum_{i=1}^{n}(\frac{(\gamma(\frac{d-u}{d}))^{2}}{n^{2}})[(3d+1-2l)(\frac{1}{p^{2}}-1)+(2l+1)(\frac{1}{p}-1)]$
$\displaystyle-\sum_{u=0}^{d}\mathbf{1}\\{u/d\geq
h\\}\binom{d}{u}\frac{2}{n}(\beta+\gamma\cdot u/d)(\gamma\cdot\frac{d-u}{d}).$

Considering the dominating terms gives us the result. ∎

###  A.5 Proof to Proposition 5

Suppose that Conditions 1, 2, 3, and 4 are satisfied. We aim to show that the
corresponding biases in the MSE terms under “ double prediction" satisfy

$\displaystyle\mathbb{E}\left[\sup_{\delta/2<|h_{BL}^{*}-h_{n}|<\delta}\sqrt{n}\left(\mathbb{M}_{n}^{b}(h_{n})-M_{n}^{b}(h_{n})-\mathbb{M}_{n}^{b}(h_{BL}^{*})-M_{n}^{b}(h_{BL}^{*})\right)\right]\to
0,$

as $\delta\to 0$, where $x_{i}$ ranges over the set $X_{i}$ of possible
fractions of degree $i$.

###### Proof.

We begin by considering the empirical process in question. We write out

$\displaystyle\mathbb{M}_{n}(h)$
$\displaystyle=\left(\frac{1}{n}\sum_{i=1}^{n}\frac{\hat{\gamma}_{h}(1-e_{i})\mathbf{1}\\{Z_{i}=1,e_{i}\geq
h\\}}{\mathbf{P}\\{Z_{i}=1,e_{i}\geq
h\\}}+\frac{1}{n}\sum_{i=1}^{n}\frac{\hat{\gamma}_{h}(1-e_{i})\mathbf{1}\\{Z_{i}=1,e_{i}\geq
h\\}}{\mathbf{P}\\{Z_{i}=1,e_{i}\geq h\\}}\right)^{2}$
$\displaystyle+\sum_{i=1}^{n}\frac{\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}Y_{i}^{2}}{n^{2}\pi_{i}^{1}}(\frac{1}{\pi_{i}^{1}}-1)$
$\displaystyle+\left.\sum_{i=1}^{n}\sum_{\mathclap{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}}^{n}\frac{\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}\mathbf{1}\\{z_{j}=1,e_{j}\geq
h\\}Y_{i}Y_{j}}{n^{2}\pi_{ij}^{11}}(\frac{\pi_{ij}^{11}}{\pi_{i}^{1}\pi_{j}^{1}}-1)\right.$
$\displaystyle+\left.\sum_{i=1}^{n}\frac{\mathbf{1}\\{z_{i}=0,e_{i}\leq
1-h\\}Y_{i}^{2}}{n^{2}\pi_{i}^{0}}(\frac{1}{\pi_{i}^{0}}-1)\right.$
$\displaystyle+\left.\sum_{i=1}^{n}\sum_{\mathclap{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}}^{n}\frac{\mathbf{1}\\{z_{i}=0,e_{i}\leq
1-h\\}\mathbf{1}\\{z_{j}=0,e_{j}\leq
1-h\\}Y_{i}Y_{j}}{n^{2}\pi_{ij}^{00}}(\frac{\pi_{ij}^{00}}{\pi_{i}^{0}\pi_{j}^{0}}-1)\right.$
$\displaystyle-\frac{2}{n^{2}}\sum_{i=1}^{n}\sum_{j\in[n];\pi_{ij}^{10}>0}(\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}\mathbf{1}\\{z_{j}=0,e_{j}\leq
1-h\\}Y_{i}Y_{j})(\frac{1}{\pi_{i}^{1}\pi_{j}^{0}}-\frac{1}{\pi_{ij}^{10}})$
$\displaystyle+\frac{2}{n^{2}}\sum_{i=1}^{n}\sum_{j\in[n];\pi_{ij}^{10}=0}(\frac{\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}Y_{i}^{2}}{2\pi_{i}^{1}}+\frac{\mathbf{1}\\{z_{j}=0,e_{j}\leq
1-h\\}Y_{j}^{2}}{2\pi_{j}^{0}}).$

The absolute Horvitz-Thompson bias estimate at threshold $h$ (and $1-h$) is

$\displaystyle\hat{b}(h)=\frac{1}{n}\sum_{i=1}^{n}\frac{\hat{\gamma}(1-e_{i})\mathbf{1}\\{Z_{i}=,e_{i}\geq
h\\}}{\mathbf{P}(Z_{i}=1,e_{i}\geq
h)}+\frac{1}{n}\sum_{i=1}^{n}\frac{\hat{\gamma}(e_{i})\mathbf{1}\\{Z_{i}=0,e_{i}\leq
1-h\\}}{\mathbf{P}(Z_{i}=0,e_{i}\leq 1-h)}$ (14)

The true-bias is

$\displaystyle b^{*}(h)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\sum_{x_{i}\in
X_{i}}\mathbf{1}\\{x_{i}\geq
h\\}y_{i}(1,x_{i})-\frac{1}{n}\sum_{i=1}^{n}\sum_{x_{i}\in
X_{i}}\mathbf{1}\\{Z_{i}=0,e_{i}\leq 1-h\\}y_{i}(0,x_{i})$
$\displaystyle-\frac{1}{n}\sum_{i=1}^{n}y_{i}(1,1)+\frac{1}{n}\sum_{i=1}^{n}y_{i}(0,0),$

where $x_{i}$ ranges over the set $X_{i}$ of possible fractions of degree $i$,

Then, for every $h$, using the identity $x^{2}-y^{2}=(x+y)(x-y)$, we have

$\displaystyle\mathbb{M}_{n}^{b}(h)-M_{n}^{b}(h)\leq
3b^{*}(h)\left(\hat{b}(h)-b^{*}(h)\right)$

We can decompose the second term on the right-hand-side above into the true
bias induced by the best possible linear fit, and the residual true bias:

$\displaystyle\hat{b}(h)-b^{*}(h)$
$\displaystyle=(\hat{b}(h)-b_{BL}^{*}(h))+(b_{BL}^{*}(h)-b^{*}(h))$

where

$\displaystyle\hat{b}(h)-b_{BL}^{*}(h)$
$\displaystyle=\left(\frac{1}{n}\sum_{i=1}^{n}\frac{\hat{\gamma}\mathbf{1}\\{Z_{i}=1,e_{i}\geq
h\\}}{\mathbf{P}(Z_{i}=1,e_{i}\geq
h)}-\left(\frac{1}{n}\sum_{i=1}^{n}y_{i}^{BL}(1,1)-\frac{1}{n}\sum_{i=1}^{n}y_{i}^{BL}(0,0)\right)\right)$
$\displaystyle-\left(\frac{1}{n}\sum_{i=1}^{n}\frac{\hat{\gamma}e_{i}\mathbf{1}\\{Z_{i}=1,e_{i}\geq
h\\}}{\mathbf{P}(Z_{i}=1,e_{i}\geq h)}-\frac{1}{n}\sum_{i=1}^{n}\sum_{x_{i}\in
X_{i}}\mathbf{1}\\{x_{i}\geq h\\}y_{i}^{BL}(x_{i})\right)$
$\displaystyle+\left(\frac{1}{n}\sum_{i=1}^{n}\frac{\hat{\gamma}e_{i}\mathbf{1}\\{Z_{i}=0,e_{i}\leq
1-h\\}}{\mathbf{P}(Z_{i}=0,e_{i}\leq
1-h)}-\frac{1}{n}\sum_{i=1}^{n}\sum_{x_{i}\in X_{i}}\mathbf{1}\\{x_{i}\leq
1-h\\}y_{i}^{BL}(x_{i})\right),$

where $y_{i}^{BL}$ is the projection of $y_{i}$ onto the best possible linear
fit. Next,

$\displaystyle b_{BL}^{*}(h)-b^{*}(h)$
$\displaystyle=\left(\frac{1}{n}\sum_{i=1}^{n}y_{i}^{BL}(1,1)-\frac{1}{n}\sum_{i=1}^{n}y_{i}(1,1)\right)$
$\displaystyle-\left(\frac{1}{n}\sum_{i=1}^{n}\sum_{x_{i}\in
X_{i}}\mathbf{1}\\{x_{i}\geq
h\\}y_{i}^{BL}(x_{i})-\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbf{1}\\{Z_{i}=1,e_{i}\geq
h\\}}{\mathbf{P}\\{Z_{i}=1,e_{i}\geq h\\}}Y_{i}\right)$
$\displaystyle+\left(\frac{1}{n}\sum_{i=1}^{n}\sum_{x_{i}\in
X_{i}}\mathbf{1}\\{x_{i}\leq
1-h\\}y_{i}^{BL}(x_{i})-\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbf{1}\\{Z_{i}=0,e_{i}\leq
1-h\\}}{\mathbf{P}\\{Z_{i}=0,e_{i}\leq 1-h\\}}Y_{i}\right)$
$\displaystyle-\left(\frac{1}{n}\sum_{i=1}^{n}y_{i}^{BL}(0,0)-\frac{1}{n}\sum_{i=1}^{n}y_{i}(0,0)\right).$

We first consider the terms in $b_{BL}^{*}(h)-b^{*}(h)$. By Conditions 4, 2,
3, each of these terms are sample-means of sub-gaussians corresponding to the
approximation error induced by a linear fit, with bounded dependency.

We now proceed to consider each of the terms in $\hat{b}(h)-b_{BL}^{*}(h)$
above. Each of the three parentheses satisfy Donsker conditions. Indeed, under
Conditions 1, 2, 3, and 4, they are sample-mean terms with uniformly bounded
coefficients, and $\hat{\gamma}$ converges at a rate of
$\mathcal{O}(n^{-1/2})$, which we compose with $3b^{*}(h)$ (see Section
3.4.3.2. in Van Der Vaart et al. [1996]). Thus, the terms in
$\mathbb{M}_{n}^{b}(h)-M_{n}^{b}(h)$ have bounded entropy integral. Under our
bounded-variation potential-outcome model (Condition 2), we have that $\log
N_{[]}(\epsilon,\Theta_{n}^{b},L_{2}(\mathbf{P})\leq
K\left(\frac{1}{\epsilon}\right)$, which gives us
$\tilde{J}(\delta,\Theta_{n}^{b},L_{2}(\mathbf{P})\leq\delta^{1/2}$ Van Der
Vaart et al. [1996] (see also Example 19.11 in Van der Vaart [2000]).

Therefore, putting these together, we have the following maximal inequality
for the bias terms

$\displaystyle\mathbb{E}\left[\sup_{\delta/2<|h_{BL}^{*}-\hat{h}_{n}|<\delta}\sqrt{n}\left(\mathbb{M}_{n}^{b}(h_{n})-M_{n}^{b}(h_{n})-(\mathbb{M}_{n}^{b}(h_{BL}^{*})-M_{n}^{b}(h_{BL}^{*}))\right)\right]$
(15)
$\displaystyle\leq\tilde{J}(\delta,\Theta_{n}^{b},L_{2}(\mathbf{P})\left(1+\frac{\tilde{J}(\delta,\Theta_{n}^{b},L_{2}(\mathbf{P})}{\delta\sqrt{n}}\right),$
(16)

with $\tilde{J}(\delta,\Theta_{n}^{b},L_{2}(\mathbf{P})\leq\delta^{1/2}$.
Additionally, in the following we consider the variance parts of the MSE
terms. The variance estimator we use following Aronow and Samii [2017] is

$\displaystyle\hat{v}(h)$
$\displaystyle=\sum_{i=1}^{n}\frac{\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}Y_{i}^{2}}{n^{2}\pi_{i}^{1}}(\frac{1}{\pi_{i}^{1}}-1)$
$\displaystyle+\sum_{i=1}^{n}\sum_{\mathclap{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}}^{n}\frac{\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}\mathbf{1}\\{z_{j}=1,e_{j}\geq
h\\}Y_{i}Y_{j}}{n^{2}\pi_{ij}^{11}}(\frac{\pi_{ij}^{11}}{\pi_{i}^{1}\pi_{j}^{1}}-1)$
$\displaystyle+\sum_{i=1}^{n}\frac{\mathbf{1}\\{z_{i}=0,e_{i}\leq
1-h\\}Y_{i}^{2}}{n^{2}\pi_{i}^{0}}(\frac{1}{\pi_{i}^{0}}-1)$
$\displaystyle+\sum_{i=1}^{n}\sum_{\mathclap{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}}^{n}\frac{\mathbf{1}\\{z_{i}=0,e_{i}\leq
1-h\\}\mathbf{1}\\{z_{j}=0,e_{j}\leq
1-h\\}Y_{i}Y_{j}}{n^{2}\pi_{ij}^{00}}(\frac{\pi_{ij}^{00}}{\pi_{i}^{0}\pi_{j}^{0}}-1)$
$\displaystyle-\frac{2}{n^{2}}\sum_{i=1}^{n}\sum_{j\in[n];\pi_{ij}^{10}>0}(\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}\mathbf{1}\\{z_{j}=0,e_{j}\leq
1-h\\}Y_{i}Y_{j})(\frac{1}{\pi_{i}^{1}\pi_{j}^{0}}-\frac{1}{\pi_{ij}^{10}})$
$\displaystyle+\frac{2}{n^{2}}\sum_{i=1}^{n}\sum_{j\in[n];\pi_{ij}^{10}=0}(\frac{\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}Y_{i}^{2}}{2\pi_{i}^{1}}+\frac{\mathbf{1}\\{z_{j}=0,e_{j}\leq
1-h\\}Y_{j}^{2}}{2\pi_{j}^{0}}).$

The true variance is

$\displaystyle v^{*}(h)$
$\displaystyle=\sum_{i=1}^{n}\frac{y_{i}(h^{+})^{2}}{n^{2}}(\frac{1}{\pi_{i}^{1}}-1)$
$\displaystyle+\sum_{i=1}^{n}\sum_{\mathclap{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}}^{n}\frac{y_{i}(h^{+})y_{j}(h^{+})}{n^{2}}(\frac{\pi_{ij}^{11}}{\pi_{i}^{1}\pi_{j}^{1}}-1)$
$\displaystyle+\sum_{i=1}^{n}\frac{y_{i}(h^{-})^{2}}{n^{2}\pi_{i}^{0}}(\frac{1}{\pi_{i}^{0}}-1)$
$\displaystyle+\sum_{i=1}^{n}\sum_{\mathclap{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}}^{n}\frac{y_{i}(h^{-})y_{j}(h^{-})}{n^{2}}(\frac{\pi_{ij}^{00}}{\pi_{i}^{0}\pi_{j}^{0}}-1)$
$\displaystyle-\frac{2}{n^{2}}\sum_{i=1}^{n}\sum_{j\in[n];\pi_{ij}^{10}>0}y_{i}(h^{+})y_{j}(h^{-})(\frac{\pi_{ij}^{10}}{\pi_{i}^{1}\pi_{j}^{0}}-1)$
$\displaystyle+\frac{2}{n^{2}}\sum_{i=1}^{n}\sum_{j\in[n];\pi_{ij}^{10}=0}y_{i}(h^{+})^{2}y_{j}(h^{-})^{2}.$

Then, we have that

$\displaystyle\sup_{h}\hat{v}(h)-v^{*}(h)$
$\displaystyle=\sup_{h}\left[\left(\sum_{i=1}^{n}\frac{\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}Y_{i}^{2}}{n^{2}\pi_{i}^{1}}(\frac{1}{\pi_{i}^{1}}-1)-\sum_{i=1}^{n}\frac{y_{i}(h^{+})^{2}}{n^{2}}(\frac{1}{\pi_{i}^{1}}-1)\right)\right.$
$\displaystyle+\left.\left(\sum_{i=1}^{n}\sum_{\mathclap{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}}^{n}\frac{\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}\mathbf{1}\\{z_{j}=1,e_{j}\geq
h\\}Y_{i}Y_{j}}{n^{2}\pi_{ij}^{11}}(\frac{\pi_{ij}^{11}}{\pi_{i}^{1}\pi_{j}^{1}}-1)\right.\right.$
$\displaystyle-\left.\left.\sum_{i=1}^{n}\sum_{\mathclap{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}}^{n}\frac{y_{i}(h^{+})y_{j}(h^{+})}{n^{2}}(\frac{\pi_{ij}^{11}}{\pi_{i}^{1}\pi_{j}^{1}}-1)\right)\right.$
$\displaystyle+\left.\left(\sum_{i=1}^{n}\frac{\mathbf{1}\\{z_{i}=0,e_{i}\leq
1-h\\}Y_{i}^{2}}{n^{2}\pi_{i}^{0}}(\frac{1}{\pi_{i}^{0}}-1)-\sum_{i=1}^{n}\frac{y_{i}(h^{-})^{2}}{n^{2}\pi_{i}^{0}}(\frac{1}{\pi_{i}^{0}}-1)\right)\right.$
$\displaystyle+\left.\left(\sum_{i=1}^{n}\sum_{\mathclap{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}}^{n}\frac{\mathbf{1}\\{z_{i}=0,e_{i}\leq
1-h\\}\mathbf{1}\\{z_{j}=0,e_{j}\leq
1-h\\}Y_{i}Y_{j}}{n^{2}\pi_{ij}^{00}}(\frac{\pi_{ij}^{00}}{\pi_{i}^{0}\pi_{j}^{0}}-1)\right.\right.$
$\displaystyle\left.\left.-\sum_{i=1}^{n}\sum_{\mathclap{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}}^{n}\frac{y_{i}(h^{-})y_{j}(h^{-})}{n^{2}}(\frac{\pi_{ij}^{00}}{\pi_{i}^{0}\pi_{j}^{0}}-1)\right)\right.$
$\displaystyle-\left.\left(\frac{2}{n^{2}}\sum_{i=1}^{n}\sum_{j\in[n];\pi_{ij}^{10}>0}(\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}\mathbf{1}\\{z_{j}=0,e_{j}\leq
1-h\\}Y_{i}Y_{j})(\frac{1}{\pi_{i}^{1}\pi_{j}^{0}}-\frac{1}{\pi_{ij}^{10}})\right.\right.$
$\displaystyle+\left.\left.\frac{2}{n^{2}}\sum_{i=1}^{n}\sum_{j\in[n];\pi_{ij}^{10}>0}y_{i}(h^{+})y_{j}(h^{-})(\frac{\pi_{ij}^{10}}{\pi_{i}^{1}\pi_{j}^{0}}-1)\right)\right.$
$\displaystyle+\left.\left(\frac{2}{n^{2}}\sum_{i=1}^{n}\sum_{j\in[n];\pi_{ij}^{10}=0}(\frac{\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}Y_{i}^{2}}{2\pi_{i}^{1}}+\frac{\mathbf{1}\\{z_{j}=0,e_{j}\leq
1-h\\}Y_{j}^{2}}{2\pi_{j}^{0}})\right.\right.$
$\displaystyle\left.\left.-\frac{2}{n^{2}}\sum_{i=1}^{n}\sum_{j\in[n];\pi_{ij}^{10}=0}y_{i}(h^{+})^{2}y_{j}(h^{-})^{2}\right)\right].$

We consider each of the parentheses above. The first and third terms are both
$o(n^{-1/2})$ since the $\frac{\mathbf{1}\\{\cdot\\}}{\mathbf{P}\\{\cdot\\}}$
terms in the summands are uniformly bounded. The second, fourth, and fifth
terms are V-statistics terms. Together with the centering terms from
$M_{n}(h)$, from Sancetta [2009], Beutner and Zähle [2012], Yoshihara [1976],
Arcones and Yu [1994] and Denker and Keller [1986], we see that under bounded
(up to fourth) moments, they have uniform $\mathcal{O}(n^{-1/2})$ convergence
rates (see also Chen et al. [2010, Theorems 9.1, 9.2]). Finally, we consider
the sixth term, and notice that this is $o(n^{-1})$, as there are a bounded
number of terms $j$, under Condition 4, for which $\pi_{ij}=0$ for each $i$.
Additionally, the smaller the difference between a fully treated outcome and a
fully control outcome, i.e. the smaller the range of the outcome, the smaller
this term is.

Similar to the bias calculations, under Condition 4, the conditions of
Sancetta [2009, Theorem 2.1], have been met, and hence, we have the following
maximal inequality from the variance terms:

$\displaystyle\mathbb{E}\left[\sup_{\delta/2<|h_{BL}-\hat{h}_{n}|<\delta}\sqrt{n}\left(\mathbb{M}_{n}^{v}(h_{n})-M_{n}^{v}(h_{n})-(\mathbb{M}_{n}^{v}(h_{BL}^{*})-M_{n}^{v}(h_{BL}^{*}))\right)\right]$
(17)
$\displaystyle\leq\tilde{J}(\delta,\Theta_{n}^{v},L_{2}(\mathbf{P})\left(1+\frac{\tilde{J}(\delta,\Theta_{n}^{v},L_{2}(\mathbf{P})}{\delta\sqrt{n}}\right),$
(18)

with $\tilde{J}(\delta,\Theta_{n}^{v},L_{2}(\mathbf{P})\leq\delta^{1/2}.$

Therefore, by a triangle inequality, the following stochastic equicontinuity
result holds:

$\displaystyle\mathbb{E}\left[\sup_{\delta/2<|h_{BL}-\hat{h}_{n}|<\delta}\sqrt{n}\left(\mathbb{M}_{n}(h_{n})-M_{n}(h_{n})-(\mathbb{M}_{n}(h_{BL}^{*})-M_{n}(h_{BL}^{*}))\right)\right]\to
0,$ (19)

as $\delta\to 0.$ ∎

###  A.6 Proof of Theorem 1

###### Proof.

Defining $K:=\mathbb{E}[\|\Sigma^{-1/2}x\text{approx}(x)\|^{2}]$,
$\hat{F}:=\left(\frac{1}{n}\sum_{i=1}^{n}\left[\frac{\mathbf{1}\\{z_{i}=1,e_{i}\geq
h\\}}{\mathbf{P}\\{z_{i}=1,e_{i}\geq
h\\}}e_{i}+\frac{\mathbf{1}\\{z_{i}=0,e_{i}\leq
1-h\\}}{\mathbf{P}\\{z_{i}=0,e_{i}\leq 1-h\\}}(1-e_{i})\right]\right)^{2}$,
$\Delta:=\min_{h,h^{\prime}:h\neq h^{\prime}}|M_{n}(h)-M_{n}(h^{\prime})|$,
$R=\max_{i}\text{deg}_{i}$, and $B$ the maximum over $i\in[n]$ of the summands
in $\mathbb{M}_{n}(h)$ and $M_{n}(h)$ with $Y_{i}-\epsilon_{i}$ in the
variance terms. Note that $B<\infty$ by Condition 4, and 2. For any
$\delta>0$:

$\displaystyle\mathbb{P}\\{|M_{n}(h)-M_{n}(h^{*})|\geq\delta\\}$
$\displaystyle\leq\mathbb{P}\\{\max_{h}|\mathbb{M}_{n}(h,\hat{\gamma})-M_{n}(h,\gamma_{0})|\geq\frac{\Delta}{2}\\}$
$\displaystyle\leq\sum_{h\in
H}\mathbb{P}\\{|\mathbb{M}_{n}(h,\hat{\gamma})-\mathbb{M}_{n}(h,\gamma_{0})+\mathbb{M}_{n}(h,\gamma_{0})-M_{n}(h,\gamma_{0})|\geq\frac{\Delta}{2}\\}$
$\displaystyle\leq
S\left(\mathbb{P}\\{|\mathbb{M}_{n}(h,\hat{\gamma})-\mathbb{M}_{n}(h,\gamma_{0})|\geq\frac{\Delta}{4}\\}+\mathbb{P}\\{|\mathbb{M}_{n}(h,\gamma_{0})-M_{n}(h,\gamma_{0})|\geq\frac{\Delta}{4}\\}\right)$
$\displaystyle=S\left(\mathbb{P}\\{|\hat{b}^{2}(h,\hat{\gamma})-\hat{b}^{2}(h,\gamma_{0})|\geq\frac{\Delta}{4}\\}+\mathbb{P}\\{|\mathbb{M}_{n}(h,\gamma_{0})-M_{n}(h,\gamma_{0})|\geq\frac{\Delta}{4}\\}\right)$
$\displaystyle\leq
S\left(\mathbb{P}\\{|\hat{\gamma}-\gamma_{0}|\geq\frac{\Delta}{12\gamma_{0}\hat{F}}\\}+\mathbb{P}\\{|\mathbb{M}_{n}(h,\gamma_{0})-M_{n}(h,\gamma_{0})|\geq\frac{\Delta}{4}\\}\right)$
$\displaystyle\leq
S\left(3\exp{\left(-\frac{n}{C_{1}(16K^{2}+2\sigma^{2})}\cdot\left(\frac{\phi^{2}\Delta}{12\gamma_{0}\hat{F}}-C_{2}\right)\right)}+\mathbb{P}\\{|\mathbb{M}_{n}(h,\gamma_{0})-M_{n}(h,\gamma_{0})|\geq\frac{\Delta}{4}\\}\right),$

where $\phi^{2}$ is the variance of $e_{i}.$ The last inequality is obtained
from Hsu et al. [2012, Theorem 11], via a simple change of variable argument.

Let $M_{n}^{f}$ and $\mathbb{M}_{n}^{f}$ denote the true and estimated mean-
squared errors with variances involving the de-noised version of $Y_{i}$, and
$M_{n}^{\epsilon}$ and $\mathbb{M}_{n}^{\epsilon}$ denote the remaining terms
in the true and estimated mean-squared errors with variances involving only
the $\epsilon_{i}$ parts. By a combinatorial Bernstein’s inequality Chatterjee
[2005], since we have bounded degree (dependence) by $R$,

$\displaystyle\mathbb{P}\\{|\mathbb{M}_{n}(h,\gamma_{0})-M_{n}(h,\gamma_{0})|\geq\frac{\Delta}{4}\\}$
$\displaystyle\leq\mathbb{P}\\{|\mathbb{M}_{n}^{f}(h,\gamma_{0})-M_{n}^{f}(h,\gamma_{0})|\geq\frac{\Delta}{16}\\}+\mathbb{P}\\{|\mathbb{M}_{n}^{\epsilon}(h,\gamma_{0})-M_{n}^{\epsilon}(h,\gamma_{0})|\geq\frac{\Delta}{16}\\}$
$\displaystyle\leq\exp{\left(-\frac{\Delta^{2}}{16^{2}\left(12RB^{2}+4\sqrt{2}B\Delta/16\right)}\right)}+\mathbb{P}\\{|\mathbb{M}_{n}^{\epsilon}(h,\gamma_{0})-M_{n}^{\epsilon}(h,\gamma_{0})|\geq\frac{\Delta}{16}\\}$

where the first inequality is obtained via the crude upper-bound
$(a+b)^{2}\leq 2a^{2}+2b^{2}$, followed by a union-bound. Finally, we upper-
bound the second term on the RHS above using the usual Bernstein inequality.
When $\Delta/16\in[0,\sqrt{2}\sigma^{2}]$, we have that

$\displaystyle\mathbb{P}\\{|\mathbb{M}_{n}^{\epsilon}(h,\gamma_{0})-M_{n}^{\epsilon}(h,\gamma_{0})|\geq\frac{\Delta}{16}\\}\leq\exp{\left(-\frac{n^{2}\Delta^{2}}{c16^{2}\cdot
2\sigma^{4}}\right)},$

and when $\Delta/16>\sqrt{2}\sigma^{2}$, we have that

$\displaystyle\mathbb{P}\\{|\mathbb{M}_{n}^{\epsilon}(h,\gamma_{0})-M_{n}^{\epsilon}(h,\gamma_{0})|\geq\frac{\Delta}{16}\\}\leq\exp{\left(-\frac{n\Delta}{c32\cdot\sigma^{2}}\right)},$

Therefore, the dominant terms are the sub-exponential decay terms, giving us
that for any $\delta>0$,

$\displaystyle\mathbb{P}\\{|M_{n}(h)-M_{n}(h^{*})|\geq\delta\\}\leq
3SC\exp{\left(-\frac{n}{C_{1}(16K^{2}+2\sigma^{2})}\cdot\left(\frac{\phi^{2}\Delta}{12\gamma_{0}\hat{F}}-C_{2}\right)\right)},$

for some constant $C$. ∎

To prove Corollary 1, we note that for $k$th power cycle graphs, $|H|=2k+1$.
Additionally, under the linear model, $h^{*}=h_{BL}^{*}.$

###  A.7 Simulations

####  A.7.1 Experimental details on the Amazon product similarity graph data

The simulations were run on a CPU. Our experiments focus on a subset of the
19828-node DVD graph. In particular, we considered the 17924-node subgraph by
removing all isolated nodes. To compute the exposure probabilities, we used
$10^{6}$ simulations. We then selected the first 1000 nodes of the 19828 to
analyze. 200 replicates were run, generating random treatment assignments and
the corresponding outcomes. For each of the replicates, 1000 separate
simulation runs were generated to compute the oracle MSEs. In total, this took
approximately 2 hours. The graph data is available at:
https://snap.stanford.edu/data/amazon-meta.html

####  A.7.2 Experimental details on simulated graphs

All synthetic graph simulations were run on a machine of Intel® Xeon®
processors with 48 CPU cores, and 50GB of RAM. We simulated 1000 replicates,
generating random treatment assignments and the corresponding outcomes, with
each oracle MSE computed using 1000 separate simulation runs. The exposure
probabilities under each threshold were computed as proposed in Aronow and
Samii [2017] using 1000 simulation iterations. In total, this took
approximately 30 minutes on average (across the different potential outcomes,
and graph settings). The code is available at:
https://anonymous.4open.science/r/AdaThresh-DDEE/

####  A.7.3 Simulations for non-linear potential outcomes models

We investigate the robustness of our estimator to potential outcomes models
that are non-linear in the exposure. In particular, we consider simulations
from the sigmoid, and sine exposure functions. In Figure 5, we display the
performance of our estimator under a sigmoid (left) and sine (right)
interference function, respectively. Our adaptive threshold Horvitz-Thompson
estimator (AdaThresh) improves upon other existing Horvitz-Thompson
estimators.

Figure 5: RMSE (normalized by the ATE) induced by the different Horvitz-
Thompson estimators. Left: 2nd-power cycle graph under unit-level Ber(0.5)
randomization under sigmoid $f(e_{i})=\gamma/(1+\exp{(-e_{i})})$. Right: 2nd-
power cycle graph under unit-level Ber(0.5) randomization with
$f(e_{i})=\gamma(1-\sin(\pi\cdot e_{i})))$ being a sine function. The error
bars are two times the standard deviation.

####  A.7.4 Simulations for Difference-in-Means estimator with adaptive
exposure thresholds

We consider our approach using the Difference-in-Means estimator incorporating
exposure thresholds:

$\displaystyle\hat{\tau}_{\text{DiM}_{h}}=\frac{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=1,e_{i}\geq
h\\}Y_{i}}{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=1,e_{i}\geq
h\\}}-\frac{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=0,e_{i}\leq
1-h\\}Y_{i}}{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=0,e_{i}\leq 1-h\\}}.$ (20)

We use the following bias estimator:

$\displaystyle\hat{b}(\hat{\tau}_{h})=\sum_{i=1}^{n}\frac{(1-e_{i})\hat{\gamma}\mathbf{1}\\{Z_{i}=1,e_{i}\geq
h\\}}{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=1,e_{i}\geq
h\\}}+\sum_{i=1}^{n}\frac{e_{i}\hat{\gamma}\mathbf{1}\\{Z_{i}=0,e_{i}\leq
1-h\\}}{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=0,e_{i}\leq 1-h\\}},$

where $\hat{\gamma}$ is the linear regression coefficient for the exposure
variable.

We use the following variance estimator:

$\displaystyle\hat{Var}(\hat{\tau})=\frac{2}{n-1}\left(\hat{Var}_{T}+\hat{Var}_{C}\right)$

where

$\displaystyle\hat{Var}_{T}=\frac{1}{n_{1}}\sum_{i=1}^{n}\left(Y_{i}\mathbf{1}\\{Z_{i}=1,e_{i}\geq
h\\}-\frac{1}{n_{1}}\sum_{i=1}^{n}Y_{i}\mathbf{1}\\{Z_{i}=1,e_{i}\geq
h\\}\right)^{2}$

where $n_{1}:=\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=1,e_{i}\geq h\\}$, and

$\displaystyle\hat{Var}_{C}=\frac{1}{n_{0}}\sum_{i=1}^{n}\left(Y_{i}\mathbf{1}\\{Z_{i}=0,e_{i}\leq
1-h\\}-\frac{1}{n_{0}}\sum_{i=1}^{n}Y_{i}\mathbf{1}\\{Z_{i}=0,e_{i}\leq
1-h\\}\right)^{2}$

where $n_{0}:=\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=0,e_{i}\leq 1-h\\}$.

We compare the performance of our adaptive estimator to the vanilla
difference-in-means estimator, the difference-in-means analogue of the vanilla
Horvitz-Thompson estimator, and the difference-in-means estimator with a
threshold plugin via Lepski’s method. We write these out below:

  * •

vanilla difference-in-means estimator

$\displaystyle\hat{\tau}_{\text{DiM}_{0}}=\frac{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=1\\}Y_{i}}{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=1\\}}-\frac{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=0\\}Y_{i}}{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=0\\}}$
(21)

  * •

difference-in-means estimator at threshold $h=1$

$\displaystyle\hat{\tau}_{\text{DiM}_{1}}=\frac{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=1,e_{i}=1\\}Y_{i}}{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=1,e_{i}=1\\}}-\frac{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=0,e_{i}=0\\}Y_{i}}{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=0,e_{i}=0\\}}$
(22)

  * •

difference-in-means estimator at threshold $\hat{h}_{\text{Lepski}}$ where

$\displaystyle\hat{h}_{\text{Lepski}}:=\min\\{h\in
H:\cap_{h=1}^{k}I(h)\neq\emptyset\\},$ (23)

with

$\displaystyle
I(h):=[\hat{\tau}_{\text{DIM}_{h}}-2\hat{\text{SDEV}}(\hat{\tau}_{\text{DIM}_{h}}),\hat{\tau}_{\text{DIM}_{h}}+2\hat{\text{SDEV}}(\hat{\tau}_{\text{DIM}_{h}})],$
(24)

and

$\displaystyle\hat{\tau}_{\text{LepskiDiM}}=\sum_{i=1}^{n}\frac{\mathbf{1}\\{Z_{i}=1,e_{i}\geq\hat{h}_{\text{Lepski}}\\}}{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=1,e_{i}\geq\hat{h}_{\text{Lepski}}\\}}Y_{i}-\sum_{i=1}^{n}\frac{\mathbf{1}\\{Z_{i}=0,e_{i}\leq
1-\hat{h}_{\text{Lepski}}\\}}{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=0,e_{i}\leq
1-\hat{h}_{\text{Lepski}}\\}}Y_{i}$ (25)

In Figures 6 and 7, we display the performance of our adaptive threshold
(AdaThresh) Difference-in-Means estimators in comparison to other estimators.
AdaThresh improves upon fixed threshold Difference-in-Means estimators.

Figure 6: RMSE (normalized by the ATE) under the linear model with
$\alpha=10$, $g(z_{i})=\beta z_{i}=10z_{i}$, $f_{0}(e_{i})=\gamma e_{i}$, and
$\epsilon_{i}\sim\mathcal{N}(0,1),$ induced by the different Difference-in-
Means estimators. Left: Cycle graph under unit-level Ber(0.5) randomization.
Center: 2nd-power cycle graph under unit-level Ber(0.5) randomization. Right:
2nd-power cycle graph under cluster-level Ber(0.5) randomization with cluster
sizes 5 (=2k + 1). The error bars are two times the standard deviation. Figure
7: RMSE (normalized by the ATE) induced by the different Difference-in-Means
estimators. Left: 2nd-power cycle graph under unit-level Ber(0.5)
randomization under sigmoid $f(e_{i})=\gamma/(1+\exp{(-e_{i})})$. Right: 2nd-
power cycle graph under unit-level Ber(0.5) randomization with
$f(e_{i})=\gamma(1-\sin(\pi\cdot e_{i})))$ being a sine function. The error
bars are two times the standard deviation.

####  A.7.5 Simulations using local linear regression

We extend our global linear regression approach to a local linear regression
one to estimate the rate of change of the bias. We write the new bias
estimators for the Horvitz-Thompson and Difference-in-Means estimators. We
illustrate the performance of the local linear regression in the settings
above, as well as settings that significantly deviates from linearity.

We write the bias estimator for the Horvitz-Thompson estimator as:

$\displaystyle\hat{b}(\hat{\tau}_{h})=\frac{1}{n}\sum_{i=1}^{n}\frac{(1-e_{i})\hat{\gamma}_{h}\mathbf{1}\\{Z_{i}=1,e_{i}\geq
h\\}}{\mathbf{P}\\{Z_{i}=1,e_{i}\geq
h\\}}+\frac{1}{n}\sum_{i=1}^{n}\frac{e_{i}\hat{\gamma}_{1-h}\mathbf{1}\\{Z_{i}=0,e_{i}\leq
1-h\\}}{\mathbf{P}\\{Z_{i}=0,e_{i}\leq 1-h\\}},$

where $\hat{\gamma}_{h}$ and $\hat{\gamma}_{1-h}$ are the local linear
regression coefficients for the exposure variable in the intervals $[h,1]$,
and $[0,1-h]$, respectively.

Similarly, write the bias estimator for the Difference-in-Means estimator as:

$\displaystyle\hat{b}(\hat{\tau}_{h})=\sum_{i=1}^{n}\frac{(1-e_{i})\hat{\gamma}_{h}\mathbf{1}\\{Z_{i}=1,e_{i}\geq
h\\}}{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=1,e_{i}\geq
h\\}}+\sum_{i=1}^{n}\frac{e_{i}\hat{\gamma}_{1-h}\mathbf{1}\\{Z_{i}=0,e_{i}\leq
1-h\\}}{\sum_{i=1}^{n}\mathbf{1}\\{Z_{i}=0,e_{i}\leq 1-h\\}},$

where $\hat{\gamma}_{h}$ and $\hat{\gamma}_{1-h}$ are the local linear
regression coefficients for the exposure variable in the intervals $[h,1]$,
and $[0,1-h]$, respectively.

Figure 8: RMSE (normalized by the ATE) under the linear model with
$\alpha=10$, $g(z_{i})=\beta z_{i}=10z_{i}$, $f(e_{i})=\gamma e_{i}$, and
$\epsilon_{i}\sim\mathcal{N}(0,1),$ induced by the different (local) Horvitz-
Thompson estimators. Left: Cycle graph under unit-level Ber(0.5)
randomization. Center: 2nd-power cycle graph under unit-level Ber(0.5)
randomization. Right: 2nd-power cycle graph under cluster-level Ber(0.5)
randomization with cluster sizes 5 (=2k + 1). The error bars are two times the
standard deviation. Figure 9: RMSE (normalized by the ATE) under the linear
model with $\alpha=10$, $g(z_{i})=\beta z_{i}=10z_{i}$, $f_{(}e_{i})=\gamma
e_{i}$, and $\epsilon_{i}\sim\mathcal{N}(0,1),$ induced by the different
(local) Difference-in-Means estimators. Left: Cycle graph under unit-level
Ber(0.5) randomization. Center: 2nd-power cycle graph under unit-level
Ber(0.5) randomization. Right: 2nd-power cycle graph under cluster-level
Ber(0.5) randomization with cluster sizes 5 (=2k + 1). The error bars are two
times the standard deviation.

Figures 8, 10, 9, and 11, display how local linear regression bias estimates
perform.

Figure 10: RMSE (normalized by the ATE) under the linear model with
$\alpha=10$, $g(z_{i})=\beta z_{i}=10z_{i}$,
$\epsilon_{i}\sim\mathcal{N}(0,1),$ induced by the different (local) Horvitz-
Thompson estimators. Left: 2nd-power cycle graph under unit-level Ber(0.5)
randomization under sigmoid $f(e_{i})=\gamma/(1+\exp{(-e_{i})})$. Right: 2nd-
power cycle graph under unit-level Ber(0.5) randomization with
$f(e_{i})=\gamma(1-\sin(\pi\cdot e_{i})))$ being a sine function. The error
bars are two times the standard deviation. Figure 11: RMSE (normalized by the
ATE) induced by the different (local) Difference-in-Means estimators. Left:
2nd-power cycle graph under unit-level Ber(0.5) randomization under sigmoid
$f(e_{i})=\gamma/(1+\exp{(-e_{i})})$. Right: 2nd-power cycle graph under unit-
level Ber(0.5) randomization with $f(e_{i})=\gamma(1-\sin(\pi\cdot e_{i})))$
being a sine function. The error bars are two times the standard deviation.


Gravity and Cosmology in Kaniadakis Statistics
N. Sadeghnezhad111Email<EMAIL_ADDRESS>
Research Institute for Astronomy and Astrophysics of Maragha (RIAAM)
University of Maragheh, P.O. Box 55136-553, Maragheh, Iran

###### Abstract

By using the Kaniadakis statistics, we discuss the modifications of Newtonian
gravity and radial velocity profile in the light of Verlinde’s formalism for
gravitational entropy. After considering the implications of
$\kappa$-statistics on the gravitational potential, it is shown that an
accelerated universe may be obtained by considering the Friedmann first
equation in this non- extensive statistics.

##  1 Introduction

The very deep and significant connection between gravitation, thermodynamics
and quantum theory, have emerged since the notable discovery of black hole
thermodynamics [1-4]. Jacobson [5] derived the Einstein field equations
starting from the Clausius relation at the horizon and using local Rindler
observers.
In order to agree a black hole with the laws of thermodynamics, it is
necessary to assign it an amount of entropy known as black hole entropy (or
Bekenstein- Hawking entropy) [2,6-7]. In this way, the four laws of black hole
mechanics [1-3,8] similar to those of thermodynamics, describe physical
properties of black holes.
These laws and subsequent efforts showed that there is a deep connection
between Einstein’s gravitational field equations and first law of
thermodynamics [5, 9]. Description of horizon area as a form of entropy and
also surface gravity as a temperature, are evidence of this connection.
Accelerated expansion of the late universe as well as primary inflationary era
can be described by the so- called fundamental relationship between the two
seemingly distinct theories.
Black hole entropy, according to the ”area theorem” [10] and Bekenstein’s
conjecture [2,6,9], should be a monotonic function of area

$S_{BH}=\frac{k_{B}A}{4l_{p}^{2}}=k_{B}\frac{Ac^{3}}{4G\hbar}.$ (1.1)

Here $A$ is the surface area of a black hole (area of the event horizon),
$k_{B}$ is the Boltzmann constant and $l_{p}$ stands for the Planck length
$\sqrt{\frac{G\hbar}{c^{3}}}$ while $G$, $\hbar$ and $c$ denote Newton’s
gravity constant, the Planck- Dirac constant and the speed of light,
respectively.
For a stationary and spherically symmetric black hole, we can write

$A=4\pi r_{h}^{2}=16\pi(\frac{GM}{c^{2}})^{2},$ (1.2)

where $r_{h}=\frac{2GM}{c^{2}}$ denotes the horizon’s radius and $M$ is the
black hole’s mass.
Bekenstein- Hawking (or Black- Hole) entropy in (1.1) is also the maximal
entropy that can be obtained by the Bekenstein bound [11]. Black hole physics
and Bekenstein bound are strong supporting evidence for the holographic
principle which states that microscopic degrees of freedom can be represented
holographically either on the boundary of space- time or on horizons. Also,
holographic bound limits how much entropy can be contained in matter and
energy occupying a specified volume of space.
Considering a microscopic structure for gravity and assuming it as an emergent
phenomenon, are the results that derived from the appeared close relationship
between gravitation and thermodynamics. In this way, entropic gravity (or
emergent gravity) was founded as a controversial theory which describes
gravity as an entropic force that is not a fundamental interaction. In support
of this approach, Verlinde regarded the gravity as a manifestation of entropy
in the universe and derived Newton’s gravitation law by this entropic force
[12]. Staring from first principles, namely equipartition law of energy and
using only space independent concepts, he showed that Newton’s law of
gravitation appear naturally and also, a relativistic generalization of the
arguments directly leads to the Einstein field equations of general
relativity.
Soon after Padmanabhan introduced the statistical origin of gravity and
connection between field equations of gravity and equipartion of energy [13],
Verlinde addressed, in a holographic scenario, the non- fundamental beginning
of gravity and suggested that gravity and space- time can indeed be explained
as emergent phenomena [12]. In other words, the tendency of systems to raise
their entropy is the origin of this emergence.
In this framework, gravity is a macroscopic effect that arises from the
changes of information regarding the positions of material bodies in the
microscopic level, which is due to the tendency of systems in increasing their
entropy [12]. In this regard, the confirmed relevance of general relativity
and thermodynamics in the black hole thermodynamics can be considered as an
example asserting this idea. Verlinde’s hypothesis has been so far a
motivation for some works on cosmology, dark energy, cosmological inflation
and cosmological acceleration.
According to Verlinde, the magnitude of the entropic force on a test particle
which moves apart from a holographic screen is given by

$F\Delta x=T\Delta S,$ (1.3)

where $\Delta x=\frac{\hbar}{mc}$, $T$ and $\Delta S=2\pi k_{B}$ are,
respectively, displacement of test particle from the screen, temperature and
entropy change on the surface. Thus, a non-vanishing temperature is
correlative with a non-zero force which, in turn, leads to a non-zero
acceleration. According to Unruh [14], a temperature equal to

$T=\frac{1}{2\pi k_{B}}\frac{\hbar a}{c},$ (1.4)

will be experienced by an observer in a frame with acceleration $a$. $T$ in
(1.4) which is the necessary temperature to cause an acceleration $a$, can be
considered as the temperature associated with the bits on the holographic
screen.
Supposing a sphere with radius $R$ as a storage device for information, the
number of used bits follows from

$N=\frac{Ac^{3}}{G\hbar},$ (1.5)

where $G$ and $A$ are Newton’s constant and area of the holographic surface,
respectively.
Assuming that the total energy of the system, $E$, is divided over the bits,
temperature of the system can be determined by the equipartition rule

$E=\frac{1}{2}Nk_{B}T.$ (1.6)

The distributed energy on the screen $E$ is a reason for the presence of some
mass $M=\frac{E}{c^{2}}$ surrounded by the screen. Thus, Verlinde recovered
Newton’s law of gravitation, practically from first principles [12]

$F=G\frac{mM}{R^{2}}.$ (1.7)

Boltzmann- Gibbs (BG) statistical mechanics and its related thermodynamics
have certainly a fundamental position in the theoretical physics. Although
statistical mechanics, as a formalism founded on the Boltzmann- Gibbs- Shannon
entropy, can successfully describe the simple and short- range equilibrium
systems, its validity and capability for complex physical, natural and
artificial systems such as non-equilibrium systems and gravitational systems
(which are long- ranged) is doubtful. So, it is useful to generalize the
statistical mechanics to the non-extensive statistical mechanics which has the
ability to overcome the limitations of BG entropy and describe some of the
complex phenomena. Study of self- gravitating systems [15], solar neutrino
deficit [16-19], peculiar velocities of galaxy clusters [20], flux of cosmic
rays [21] and some cosmological effects [22] are samples for application of
non- extensive statistical mechanics in astronomy and astrophysics.
The well- known Kaniadakis statistics [23-25], also known as
$\kappa$-statistics among the other models, generalizes the standard BG
statistics and can describe a large class of experimentally observed
phenomena. This framework, which emerges from the relativistic formulation of
statistical mechanics, has succeeded when applied in various cases for
example, analyzing the spectrum of the cosmic rays [26-27] and cosmic effects
[28].
In this paper, we apply the $\kappa$\- statistics in Verlinde’s framework of
gravitational entropy. Thereinafter, we consider modifications to the
Newtonian gravity and implications of $\kappa$-statistics on the radial
velocity profile and gravitational potential. We will see that it is possible
to obtain an accelerated universe by considering the Friedmann first equation
in the $\kappa$\- framework.
This paper is outlined as follows. In section 2, we will provide a brief
review of Kaniadakis formulation. In section 3, we employ the $\kappa$\-
statistics to derive the modified Newtonian second law of motion. In section
4, we will investigate some cosmological aspects of $\kappa$\- statistics. The
conclusions will be depicted in the last section.

##  2 Kaniadakis Statistics

Formally, the $\kappa$\- framework is defined by one- parameter deformation of
the exponential and logarithm functions as

$\exp_{\kappa}(f)=(\sqrt{1+\kappa^{2}f^{2}}+\kappa f)^{\frac{1}{\kappa}},$
(2.1) $\ln_{\kappa}(f)=\frac{f^{\kappa}-f^{-\kappa}}{2\kappa},$ (2.2)

where deformation parameter $\kappa$ is an unknown parameter [23-25]. The
constructed statistical mechanics recovers the ordinary Boltzmann- Gibbs
statistical mechanics at the $\kappa\rightarrow 0$ limit.
The statistical distribution and entropy associated to this $\kappa$\-
framework can be written as

$f_{\kappa}=\exp_{\kappa}(-\beta[E+\mu]),$ (2.3)
$S_{\kappa}=-k_{B}\sum_{i}^{W}\frac{p_{i}^{1+\kappa}-p_{i}^{1-\kappa}}{2\kappa},$
(2.4)

which reduce to the standard results in the limit $\kappa=0$. Specially,
$\kappa$\- entropy in (2.4) leads to Gibbs entropy

$S=-k_{B}\sum_{i=1}^{W}p_{i}\ln(p_{i}),$ (2.5)

at the $\kappa\rightarrow 0$ limit while $p_{i}$ is the probability for
occupying the $i$th state and $W$ is the number of states in distribution. Non
specified parameter $\beta$ in (2.3) contains all the information about the
temperature of the system.
Using the kinetic foundation of $\kappa$\- statistics [23-24, 29-30], the
$\kappa$\- equipartition theorem can be calculated through the relation

$E_{\kappa}=\frac{1}{2}Nf(\kappa)k_{B}T,$ (2.6)

where $N$ is the total number of the degrees of freedom for the system and

$f(\kappa)=\frac{(1+\frac{\kappa}{2})}{2\kappa(1+\frac{3\kappa}{2})}\frac{\Gamma(\frac{1}{2\kappa}-\frac{3}{4})\Gamma(\frac{1}{2\kappa}+\frac{1}{4})}{\Gamma(\frac{1}{2\kappa}+\frac{3}{4})\Gamma(\frac{1}{2\kappa}-\frac{1}{4})}\quad
0\leq\kappa<\frac{2}{3}.$ (2.7)

It is easy to see that the classical equipartition theorem for each
microscopic degrees of freedom is recovered in the $\kappa=0$ limit
($\lim_{\kappa\rightarrow
0}S_{\kappa}=S_{BH},E_{\kappa=0}=\frac{1}{2}Nk_{B}T$), while the expression
(2.7) diverges for $\kappa=\frac{2}{3}$.
In a microcanonical ensemble, which all the states have the same probability
$p_{i}=\frac{1}{W}$ , Gibbs entropy in (2.5) leads to Boltzmann entropy

$S=k_{B}\ln W,$ (2.8)

and using Bekenstein entropy (1.1) we have

$S_{BH}=k_{B}\frac{Ac^{3}}{4G\hbar}=k_{B}\ln W.$ (2.9)

In the so called ensemble, on the other hand, $\kappa$\- entropy reduces to

$S_{\kappa}=k_{B}\frac{W^{\kappa}-W^{-\kappa}}{2\kappa},$ (2.10)

where in combining with (2.9) and assuming $k_{B}=1$ , we have

$S_{\kappa}=\frac{1}{\kappa}\sinh(\kappa
S_{BH})=\frac{1}{\kappa}\sinh(\frac{A\kappa c^{3}}{4G\hbar}),$ (2.11)

found out previously in the [31].

##  3 Dynamics in Kaniadakis Gravity

Now, in the $\kappa$\- statistics scenarios, the net applied force to particle
$m$ by source $M$

$F^{\kappa}=T\frac{dS_{\kappa}}{dA}\frac{\Delta A}{\Delta x},$ (3.1)

can be written as the following form by using the Eq.s (1.4), (1.5) and (2.6)

$F^{\kappa}=ma_{\kappa}\cosh(\frac{\kappa}{f(\kappa)}\frac{\pi
Mc^{3}}{\hbar}\frac{1}{a_{\kappa}}).$ (3.2)

Considering that the total energy on the surface in (2.6) is in fact the
relativistic rest mass of the source mass $E=Mc^{2}$ and using Eq.s (1.4) and
(1.5) we easily obtain

$a^{\kappa}=G_{\kappa}\frac{M}{R^{2}}:G_{\kappa}=\frac{G}{f(\kappa)},$ (3.3)

where $a^{\kappa}$ and $G_{\kappa}$ are the modified gravitational
acceleration of test particle and effective gravitational constant in
$\kappa$\- framework, respectively.
Comparing $a^{\kappa}$ with Newtonian acceleration $a^{N}=\frac{GM}{R^{2}}$ ,
we see that they have the same sign because of the positivity of $f(\kappa)$
in the range of $\kappa$ . So, these two different frameworks expect a similar
behavior for the acceleration.
Now, we can compare the modified gravitational force in $\kappa$\- framework
(3.2), with Newtonian force and obtain

$\frac{F^{\kappa}}{F^{N}}=\cosh(GM\frac{\kappa}{R^{2}}).$ (3.4)

This expression tends to unity at large distances, without any divergence, in
the allowed range of $\kappa$. In addition, $F^{\kappa}$ also does not vanish
in this range.
In a circular motion, and using the acceleration $a^{\kappa}$ in (3.3), we
easily reach at

$v=\sqrt{\frac{G_{\kappa}M}{r}},$ (3.5)

where $r$ and $v$ are respectively the radius and velocity of motion.
$G_{\kappa}$ in (3.5) must be positive, to have a real value for velocity,
which that is the case for positive $f(\kappa)$ in the range
$0<\kappa<\frac{2}{3}$.
The source mass $M$, also, has a gravitational field which applies the force
$\frac{GM}{r^{2}}$ on the test particle $m$ . Equating (3.2) with this force
yields

$\frac{GM}{r}=v^{2}\cosh(d\frac{r}{v^{2}}):d=\frac{\kappa}{f(\kappa)}\frac{\pi
Mc^{3}}{\hbar},$ (3.6)

which by expanding the $\cosh(d\frac{r}{v^{2}})$ and approximating it up to
first order, we obtain the velocity profile as

$v^{2}\cong\frac{GM}{r}+\frac{d^{2}}{2GM}r^{3}.$ (3.7)

The second term in rhs of (3.7) causes an increase for the velocity of test
particle in comparison to the Newtonian one.

##  4  Cosmological Aspects in Kaniadakis Framework

Now, we want to calculate the Friedmann first equation in the acquired
$\kappa$\- gravity, by following the used recipe in [31]. Using the Eq.s (3.2)
and (3.3), one finds for the gravitational potential

$\phi(r)=\int\frac{F^{\kappa}}{m}dr=MG_{\kappa}[-\frac{1}{r}+\frac{1}{6}(\frac{\kappa\pi
c^{3}}{\hbar G})^{2}r^{3}+O(r^{6})],$ (4.1)

which can be approximated up to $r^{3}$ as

$\phi(r)\cong-\frac{MG_{\kappa}}{r}(1-\frac{\alpha^{2}}{6}\frac{1}{r^{4}}),$
(4.2)

where $\alpha=\frac{\kappa\pi c^{3}}{\hbar G}$ and $r=\frac{1}{H}$ is the
location of apparent horizon. Following the approach [32] one can find the
modified Friedemann equation. Since we intend to regard the test particle on
the edge of a flat FRW universe, we use the (1.6) in [32] to find

$H^{2}=\frac{8\pi G_{\kappa}}{3}\rho(1-\frac{\alpha^{2}}{6H^{4}}),$ (4.3)

where $\rho$ is the energy density of cosmic fluid. This modified Friedmann
equation converts to standard one $H^{2}=\frac{8\pi G}{3}\rho$ in the limit
$\kappa\rightarrow 0(G_{\kappa}\rightarrow G)$.

##  5 Summary

In this work, by considering the Kaniadakis statistics ($\kappa$-statistics)
in the framework of Verlinde’s hypothesis, we studied some modifications of
Newtonian gravity. Applying of this non- extensive statistics to radial
velocity profile shows an increase for the velocity of test particle in
comparison with Newtonian one. Next, we explored the modification of
gravitational potential in $\kappa$-statistics. Finally, we obtained the
Friedmann first equation in a classical approach which results an accelerated
universe and converts to standard one when the parameter of generalized
entropy tends to zero.

## References

  * [1] J. M. Bardeen, B. Carter, and S. W. Hawking, Commun. Math. Phys. 31, 161 (1973).
  * [2] J. D. Bekenstein, Phys. Rev. D 7, 2333 (1973).
  * [3] S. W. Hawking, Commun. Math. Phys. 43, 199 (1975).
  * [4] S.W. Hawking, Commun. Math. Phys. 46, 206 (1976).
  * [5] T. Jacobson, Phys. Rev. Lett. 75, 1260 (1995).
  * [6] Bekenstein J.D., Lett. Nuovo Cimento 4, 737 (1972)).
  * [7] S. W. Hawking, Nature 248, 30 (1974).
  * [8] J. D. Bekenstein, Phys. Rev. D 9, 3292 (1974).
  * [9] S. A. Hayward, Class. Quantum Grav. 15, 3147 (1998).
  * [10] S. W. Hawking, Phys. Rev. Lett. 26, 1344 (1971).
  * [11] J. D. Bekenstein, Phys. Rev. D 23, 287 (1981).
  * [12] E. Verlinde, JHEP 1104, 029 (2011).
  * [13] T. Padmanabhan, Mod. Phys. Lett. A 25, 1129 (2010).
  * [14] W.G. Unruh, Phys. Rev. D 14, 870 (1976).
  * [15] A. R. Plastino and A. Plastino, Phys. Lett. A 174, 384 (1993).
  * [16] G. Kaniadakis, A. Lavagno, M. Lissia and P. Quarati, Proc. 5th International Workshop on Relativistic Aspects of Nuclear Physics, Rio de Janeiro-Brazil, 1997,World Scientific, Singapore, p. 193.
  * [17] G. Kaniadakis, A. Lavagno and P. Quarati, Phys. Lett. B 369, 308 (1996).
  * [18] G. Kaniadakis, A. Lavagno and P. Quarati, Astrophysics and Space Science 258, 145 (1998).
  * [19] A. Lavagno and P. Quarati, Nucl. Phys. B, Proc. Suppl. 87, 209 (2000).
  * [20] Lavagno, G. Kaniadakis, M. Rego-Monteiro, P. Quarati and C. Tsallis, Astrophysical Letters and Communications 35, 449 (1998).
  * [21] C. Tsallis, J.C. Anjos and E.P. Borges, Phys. Lett. A 310, 372 (2003).
  * [22] V.H. Hamity and D.E. Barraco, Phys. Rev. Lett. 76, 4664 (1996).
  * [23] G. Kaniadakis, Physica A 296, 405 (2001).
  * [24] G. Kaniadakis, Phys. Rev. E 66, 056125 (2002).
  * [25] G. Kaniadakis, Entropy 15, 3983 (2013).
  * [26] G. Kaniadakis, P. Quarati and A. Scarfone, Physica A 305, 76 (2002).
  * [27] G. Kaniadakis and A. M. Scarfone, Physica A 305, 69 (2002).
  * [28] E. M. C. Abreu, J. A. Neto, E. M. Barboza and R. C. Nunes, Physica A 441, 141 (2016).
  * [29] E.P. Bento , J.R.P. Silva , R. Silva, Physica A 392, 666 (2013).
  * [30] G. Kaniadakis, Phys. Lett. A 288, 283 (2001).
  * [31] H. Moradpour, A. H. Ziaie, and M. Kord Zangeneh, Eur. Phys. J. C 80: 732 (2020).
  * [32] M. Kord Zangeneh, H. Moradpour and N. Sadeghnezhad, Mod. Phys. Lett. A 34, 1950168 (2019).