# Sensitivity of Bayesian Casual Forests to Modeling Choices:
A Re-analysis of the 2022 ACIC Data Challenge

Ajinkya H. Kokandakar, Hyunseung Kang11footnotemark: 1, and Sameer K.
Deshpande11footnotemark: 1 Dept. of Statistics, University of
Wisconsin–Madison. Correspondence to<EMAIL_ADDRESS>
###### Abstract

We demonstrate how Hahn et al.’s Bayesian Causal Forests model (BCF) can be
used to estimate conditional average treatment effects for the longitudinal
dataset in the 2022 American Causal Inference Conference Data Challenge.
Unfortunately, existing implementations of BCF do not scale to the size of the
challenge data. Therefore, we developed flexBCF — a more scalable and flexible
implementation of BCF — and used it in our challenge submission. We
investigate the sensitivity of our results to two ad hoc modeling choices we
made during our initial submission: (i) the choice of propensity score
estimation method and (ii) the use of sparsity-inducing regression tree
priors. While we found that our overall point predictions were not especially
sensitive to these modeling choices, we did observe that running BCF with
flexibly estimated propensity scores often yielded better-calibrated
uncertainty intervals.

##  1 Introduction

Given pairs of observed predictor vectors $\bm{x}\in\mathbb{R}^{p}$ and scalar
outputs $y\in\mathbb{R},$ Chipman et al., (2010)’s Bayesian additive
regression trees (BART) model approximates the regression function
$f(\bm{x}):=\mathbb{E}[y|\bm{x}]$ with a sum of binary regression trees. BART
often delivers accurate point estimates and well-calibrated uncertainty
intervals around evaluations $f(\bm{x})$ without requiring users to (a) pre-
specify the function al form of $f$ or (b) tune any hyperparameters. BART’s
ease-of-use and generally excellent tuning-free estimation, prediction, and
uncertainty quantification have made BART a popular “off-the-shelf” modeling
tool.

As Hill, (2011) noted, these advantages make BART an attractive choice for
estimating heterogeneous treatment effects. To wit, conditional average
treatment effects may be non-linear and may depend on complicated interactions
between confounders, covariates, and treatment. Correctly specifying such non-
linearity and interactions is difficult — if not impossible — in most
practical settings, making BART particularly compelling. Given observed
triplets $(\bm{x},y,z)$ of covariates $\bm{x}$, outcome $y$, and treatment
indicator $z,$ Hill, (2011) used BART to obtain an estimate
$\hat{f}(\bm{x},z)$ of the response surface $\mathbb{E}[y|\bm{x},z]$ before
estimating treatment effects as $\hat{f}(\bm{x},1)-\hat{f}(\bm{x},0).$ In past
iterations of the ACIC Data Challenge, BART-based models have consistently
ranked among the top-performing methods. One method in particular — Hahn et
al., (2020)’s Bayesian causal forests (BCF) model — stands out as particularly
adept at estimating heterogeneous treatment effects. In this paper, we
describe our experience deploying BCF to the 2022 ACIC Data Challenge.

At a high-level, BCF models the expected outcome as

$\mathbb{E}[Y|\bm{x},z]=\mu(\bm{x})+z\times\tau(\bm{x})$ (1)

and further approximates the prognostic function $\mu$ and the treatment
effect function $\tau$ using separate regression tree ensembles. When fitting
the tree ensembles, BCF includes an estimate of the propensity score as an
additional covariate. Although these ensembles are a priori independent, they
are dependent a posteriori, as BCF leverages all observations (both treated
and control) to learn $\mu.$

At first glance, BCF appears ill-suited to the 2022 ACIC Data Challenge: BCF
was initially proposed for cross-sectional data and the Data Challenge
features longitudinal data with time-varying covariates and outcomes. However,
as we show in Section 2.2, under a mild assumption about the data generating
process, the Data Challenge admits an additive decomposition similar to the
one in Equation (1). Also, while our analysis plan enables one to use BCF to
estimate the treatment effect function — and subsequently the desired subgroup
average treatment effects — we found that the implementation available in the
R package bcf did not scale to the challenge data. We therefore re-implemented
a version of BCF that can scale to the Data Challenge.

When writing our implementation and subsequently deploying it during the
competition, we made a number of essentially ad hoc decisions. First, we
included an estimate of Imai and Ratkovic, (2014)’s covariate balancing
propensity score (CBPS) as an additional covariate in the expected outcome
model. Second, in our regression tree prior, we selected splitting variables
using Linero, (2018)’s sparsity-inducing Dirichlet prior instead of the
uniform prior that is used in the bcf package. Briefly, Linero, (2018)’s prior
encourages regression trees to split on a small number of covariates, thereby
facilitating a type of automatic variable selection. In the context of
treatment effect estimation, such a sparsity-inducing prior can potentially
identify the main drivers of effect heterogeneity. We assess the sensitivity
of our contest submissions to these modeling choices. Specifically, we study
the extent to which our results change when we use different estimates of the
propensity score and when we use Chipman et al., (2010)’s uniform prior on the
splitting variables in the trees.

Here is an outline for the rest of the paper. We briefly introduce the
competition dataset in Section 2 before presenting our identification analysis
and discussing our new implementation of BCF. Then, in Section 3, we compare
the estimates for several configurations of the BCF model. Finally, we discuss
the implications of our findings in Section 4.

##  2 Identification and estimation

In this section, we review the 2022 ACIC Data Challenge’s study design and
identification strategy. We then outline our estimation strategy.

###  2.1 Review: Study design, notation, and definitions

The 2022 ACIC Data Challenge involved analyzing a total of 3,400 synthetic
datasets consisting of 200 independent realizations of 17 different data
generating processes. Each synthetic dataset simulated observations of
beneficiary-level Medicare expenditure data over a period of four years. Each
beneficiary received care at one of about 500 different medical practices. In
addition to simulated time-invariant beneficiary- and practice-level
covariates, each dataset also included average measurements of patient
covariates within each practice. In each data generating process, treatment
was assigned to some practices after two years. The main task of the Challenge
was to assess how the treatment deployed at the practice level affected future
medical expenditures of beneficiaries in the treated practices.

For each beneficiary $i$ in practice $j$, let $Y_{ijt}$ denote their outcome
(Medicare expenditure) in year $t$. Also, for each practice $j$, let
$Z_{j}\in\\{0,1\\}$ denote its treatment status where $Z_{j}=1$ indicates that
practice $j$ was assigned treatment between year 2 and year 3 and let
$Z_{j}=0$ indicates that practice $j$ was not treated during that period. As
the notation suggests, the treatment is assigned only once, in between year 2
and year 3, and before assigning treatment, all beneficiaries are untreated.
Also, treatment is assigned at the practice level and all beneficiaries in a
practice received the same treatment condition. Hence, we omit subscripts $i$
and $t$ in $Z_{j}$ for notational simplicity111We can also define treatment
status to be $Z_{jt}\in\\{0,1\\}$ where $Z_{jt}=1$ if practice $j$ was treated
in year $t$ and $Z_{jt}=0$ if practice $j$ was not treated in year $t$. We can
then encode the aforementioned constraints on the treatment status by
assuming: (i) for every practice $j$ and year $t=1,2$,
$\mathbb{P}(Z_{jt}=0)=1$ and (ii) for every practice $j$,
$\mathbb{P}(Z_{j3}=Z_{j4})=1$.. Finally, for each beneficiary $i$ in practice
$j$, let $\bm{X}_{ij}\in\mathbb{R}^{p}$ be all observed, “pre-treatment
covariates,” i.e. variables measured at years $1$ and $2$ before treatment is
assigned. These covariates can include time-invariant or time-varying
variables in year 1 and 2, either at the practice level or at the beneficiary
level. Again, for notational simplicity, we omit the subscript $t$ in
$\bm{X}_{ij}$222Again, we can also define pre-treatment covariates by
partitioning them into time-varying variables and time-invariant variables,
similar to the notation used to define causal effects in longitudinal settings
(Hernan and Robins,, 2020, Chapter 19). Specifically, for each beneficiary $i$
in practice $j$, let $\bm{X}_{ijt}$ denote their vector of time-varying
variables in year $t=1,2$. We also let $\bm{X}_{ij0}$ denote their vector of
time-invariant variables. Then, we can let
$\overline{\bm{X}}_{ij2}=(\bm{X}_{ij0},\bm{X}_{ij1},\bm{X}_{ij2})$ be the
collection of all pre-treatment covariates up to (and including) year 2..

We use potential outcomes to define causal effects (Splawa-Neyman et al.,,
1990; Rubin,, 1974). For each beneficiary $i$ in practice $j$, let
$Y_{ijt}(1)$ denote their potential outcome at time $t$ when their practice
was treated between years 2 and 3. Similarly, let $Y_{ijt}(0)$ denote their
potential outcome at time $t$ when their practice was untreated between years
2 and 3333The potential outcomes notation makes some implicit assumptions,
notably (i) treatment status of practice $j$ cannot affect the potential
outcome of beneficiaries in practice $j^{\prime}\neq j$, (ii) there are no
different versions of treatment (or non-treatment), (iii) treatment status
does not vary amongst beneficiaries in the same practice, and (iv) treatment
status does not change once assigned in between years 2 and 3. Conditions (i)
and (ii) are often known as the Stable Unit Treatment Value Assumption
(Rubin,, 1980, 1986).

We focus on the conditional average treatment effect on the treated (CATT),
defined as

$\textrm{CATT}(\bm{x},t)=\mathbb{E}\left[Y_{ijt}(1)-Y_{ijt}(0)|Z_{j}=1,\bm{X}_{ij}=\bm{x}\right],\quad{}t\in\\{3,4\\}.$
(2)

Note that taking the average with respect to the distribution of the
covariates $\bm{X}_{ij}$ gives us the overall average treatment effect on the
treated (ATT), i.e. ${\rm
ATT}(t)=\mathbb{E}\left[\textrm{CATT}(\bm{X}_{ij},t)\right]$. For both
estimands (and contrary to the Data Challenge), the expectation is taken over
a hypothetical super-population of practices and beneficiaries at each time
$t$ in order to better align with the Bayesian paradigm we pursue in the
paper444The super-population we have in mind is based on “model-based”
sampling (Särndal,, 1978) where for each time $t$, the entire collection of
random variables
$(Y_{11t}(1),Y_{11t}(0),Z_{1},\bm{X}_{11},\ldots,Y_{n_{jt}Jt}(1),Y_{n_{jt}Jt}(0),Z_{J},\bm{X}_{n_{jt}J})$
follows some distribution $F$ and the expectation is over the distribution
$F$. For example, suppose at each time $t$, we have
$(Y_{ijt}(1),Y_{ijt}(0),Z_{j},\bm{X}_{ij})\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}G$
where $G$ can be a parametric distribution or a mixture of parametric
distributions, say a mixture of $J$ parametric densities where each density
represents one of the $J$ practices. Then $F=\prod_{ij}G$ and the expectation
$\mathbb{E}\left[Y_{ijt}(1)-Y_{ijt}(0)|Z_{j}=1,\bm{X}_{ij}=\bm{x}\right]$,
which is over $F$, is equivalent to the expectation over $G$. Historically,
model-based sampling has been associated with Bayesian paradigms as it allows
investigators to encode subject-matter expertise by imposing constraints on
$F$; see, e.g., Ericson, (1969), Särndal, (1978, §5), Gelman and Hill, (2006,
Chapter 21.2), and and Schochet, (2013) for further discussions. Under
appropriate conditions, an estimator (or a test) of the super-population
estimand is also reasonable for estimating the finite-sample analog of the
estimand; see Lehmann and Romano, (2005, pg. 188) for one example and Ding et
al., (2017)..

###  2.2 Review: Causal Identification

To identify the CATT from the observed data, the Data Challenge made the
following three assumptions for all beneficiaries and practices:

  * (A1)

For all $t$, $Y_{ijt}=Y_{ijt}(0)\times\mathbbm{1}\left(t\leq
2\right)+Y_{ijt}(Z_{j})\times\mathbbm{1}\left(t>2\right)$.

  * (A2)

For each $s\in\\{1,2\\}$ and $t\in\\{3,4\\},$
$(Y_{ijt}(0)-Y_{ijs}(0))\perp\kern-6.00006pt\perp Z_{j}\mid\bm{X}_{ij}$.

  * (A3)

There is a $\delta>0$ such that $\mathbb{P}(Z_{j}=1)>\delta$ and
$\mathbb{P}(Z_{j}=1|\bm{X}_{ij}=\bm{x})<1-\delta$ for all $\bm{x}$.

Under assumptions (A1)–(A3), for all $s=\\{1,2\\}$ and $t\in\\{3,4\\},$ we can
write $\textrm{CATT}(\bm{x},t)$ as a function of the observed data using
standard arguments from difference-in-differences (e.g. Heckman et al.,, 1997;
Abadie,, 2005):

$\displaystyle\begin{split}\textrm{CATT}(\bm{x},t)\ =\
&\mathbb{E}\left[Y_{ijt}\mid
Z_{j}=1,\bm{X}_{ij}=\bm{x}\right]-\mathbb{E}\left[Y_{ijt}\mid
Z_{j}=0,\bm{X}_{ij}=\bm{x}\right]\\\ &-\Big{\\{}\mathbb{E}\left[Y_{ijs}\mid
Z_{j}=1,\bm{X}_{ij}=\bm{x}\right]-\mathbb{E}\left[Y_{ijs}\mid
Z_{j}=0,\bm{X}_{ij}=\bm{x}\right]\Big{\\}}.\end{split}$ (3)

We provide a full proof in Appendix A. Additionally, using Equation (3)
“twice” under different values of $s=1,2$, we observe

$\displaystyle\mathbb{E}\left[Y_{ij2}-Y_{ij1}\mid
Z_{j}=1,\bm{X}_{ij}=\bm{x}\right]$
$\displaystyle=\mathbb{E}\left[Y_{ij2}-Y_{ij1}\mid
Z_{j}=0,\bm{X}_{ij}=\bm{x}\right]$ (4)

In words, Equation (4) roughly states that conditional on covariates, the
outcome trend before year 3 is identical between the treated and the untreated
groups. This additional structure implied by the assumptions from the Data
Challenge is the basis for our modeling assumption in Section 2.3 that
eventually allows us to utilize a faster version of BCF.

###  2.3 Estimation & Implementation

We will additionally assume that for every $(t,\bm{x},z)$, we can write

$\mathbb{E}\left[Y_{ijt}|\bm{X}_{ij}=\bm{x},Z_{j}=z\right]=\mu(\bm{x},t)+z\times\mathbbm{1}\left(t>2\right)\times\tau(\bm{x},t).$
(5)

Under Equation (5) and the assumptions introduced in Section 2.2, we have
$\textrm{CATT}(\bm{x},t)=\tau(\bm{x},t).$ Also under Equation (5), we have
$\mathbb{E}\left[Y_{ij2}-Y_{ij1}|Z_{j}=z,\bm{X}_{ij}=\bm{x}\right]=\mu(\bm{x},2)-\mu(\bm{x},1)$
for $z\in\\{0,1\\}$, satisfying Equation (4). In contrast to Assumption (A2),
which was about trends in outcomes, Equation (5) makes an assumption directly
about the outcome at each time $t.$ A benefit of our additional assumption is
that the functional forms of Equations (1) and (5) are now similar, making it
possible to use BCF to estimate the $\textrm{CATT}.$ We note that BCF allows
$\mu$ and $\tau$ to depend on different covariates. We included the average of
the patient-level covariates at both pre-treatment time points as covariates
for $\mu$ and similar averages at all four time points as covariates for
$\tau.$ Although patient-level covariates are time-invariant, practice
composition generally changed over time as patients attrited from the study.
Consequently, our conditioning on these post-treatment covariates can lead to
biased estimation of the CATT unless patient attrition is not impacted by
treatment (Rosenbaum,, 1984). For our initial submission and in our subsequent
re-analysis, we assumed that patient attrition — and therefore practice
composition — was not impacted by treatment.

Unfortunately, we encountered several problems when using the bcf package to
fit the model in Equation (5) to the beneficiary-level data. First, the
implementation in that package begins by internally copying all of the
covariate data into a new matrix and allocating a matrix to hold posterior
samples of evaluations of $\tau$ for all observations. When modeling the
beneficiary-level data, we included the practice label as an additional
categorical covariate. bcf converted the practice label into several hundred
binary indicators (one per practice). Consequently, bcf allocated around 5GB
of memory to store a copy of each beneficiary-level dataset. Further, bcf
attempted to allocate around 20GB of memory for a matrix storing 2000
posterior samples of the effect of treatment on each beneficiary. These memory
requirements were well beyond the scope of our personal machines. When we
attempted to run BCF on the full data on a high-memory computer, we found that
the sampler was extremely slow. On further inspection, we discovered that the
posterior updates of each regression tree involved several redundant
computations. Basically, while updating a single tree, bcf loops over the
entire dataset several times to identify the leaf to which every observation
is assigned. Because the underlying Gibbs sampler only changes at most two
leaves at a time, the map from observations to tree leaves does not change
much iteration to iteration, rendering many of these loops redundant.

We developed a more efficient implementation of BCF that (i) avoids making
unnecessary copies of the data; (ii) does not make redundant calculations when
updating the trees; (iii) does not instantiate a matrix to hold all posterior
draws of $\tau(\bm{x},t)$; and (iv) provides users the option to use Linero,
(2018)’s sparsity-inducing prior on the splitting variables used by each tree.
Our sampler returns a list of character vectors whose elements contain string
representations of each regression tree in the ensemble used to approximate
$\tau.$ We wrote a new predict routine that takes these character vectors to
output posterior samples of evaluations of $\tau.$ Our implementation is
available in an R package that we call flexBCF. Compared to bcf, flexBCF is
much faster and much less memory-intensive: our analyses ran in a few hours on
shared computers in a high-throughput computing cluster and we needed only
about 4GB of memory to analyze each beneficiary-level dataset with flexBCF.

##  3 Re-analysis of Challenge Data

In our submission to the ACIC Data Challenge, we ran flexBCF using a sparsity-
inducing prior on splitting variables and with an estimated covariate
balancing propensity score. These choices were, admittedly, rather ad hoc. To
see whether our submitted results change substantially under different
choices, we re-analyzed the competition data using different propensity score
estimation methods and using (or not) the sparsity-inducing prior on splitting
variables. In addition to the covariate-balancing propensity score estimate
(hereafter CBPS), we considered propensity scores estimated with BART
(hereafter BART), $L_{1}$-regularized logistic regression (hereafter LASSO),
and gradient boosting (hereafter GBM). We computed the four propensity score
estimates using the R packages CBPS (Fong et al.,, 2022), BART (Sparapani et
al.,, 2021), glmnet (Simon et al.,, 2011), and GBM (Greenwell et al.,, 2022).
Note that CBPS and LASSO assume that the log-odds of treatment are linear in
the covariates $\bm{x}$ while BART and GBM do not make any a priori
assumptions about the functional form of the log-odds of treatment.
Accordingly, we refer to CBPS and LASSO as parametric and BART and GBM as non-
parametric propensity score estimators. We will refer to the model used in our
initial Challenge submission as CBPS(S) where the suffix (S) refers to the
sparsity-inducing prior. We will use the suffix (D) to refer to the uniform
prior on splitting variables.

We re-analyzed the practice-level data using all eight configurations of
propensity score and sparsity-inducing prior and we re-analyzed the
beneficiary-level data using only CBPS(S), CPBS(D), GBM(S), and GBM(D). We
compared the root mean square error (RMSE) and two uncertainty interval
diagnostics (coverage and relative interval length to CBPS(S)) for (i) the
overall ATT (i.e. the estimated CATT’s averaged over all beneficiaries and
practices) and (ii) twelve subgroups effects defined as the estimated CATT
averaged over different pre-specified subgroups defined by pre-treatment
practice covariates. Throughout, we report these performance measures within
several categories of DGPs defined by the amount of confounding and impact
heterogeneity. Specifically, we consider five categories of DGPs: (i) no
confounding but large impact heterogeneity; (ii) weak confounding with small
impact heterogeneity; (iii) weak confounding with large impact heterogeneity;
(iv) strong confounding with small impact heterogeneity; and (v) strong
confounding with large impact heterogeneity.

###  3.1 Practice level results

Overall ATT. Figure 1 compares the RMSEs of all methods for the overall ATT
for each DGP broken down by category. In general, as the amount of confounding
and treatment effect heterogeneity increased, the RMSE of all methods
increased. Further, we did not observe substantial differences between the
methods within each category of DGP.

Figure 1: Distribution of RMSEs for the overall ATT for each DGP category
using practice-level data

Interestingly, in the absence of confounding, using a nonparametric propensity
score estimate (i.e. BART or GBM) yielded slightly larger RMSE than using
parametric estimates (CBPS or L1-regularized logistic regression). However, in
the presence of even weak confounding, the nonparametric propensity score
estimates fared slightly better. That said, the difference in performance is
quite small and is almost certainly within the range of Monte Carlo
variability. Additionally, using the sparsity-inducing prior did not seem to
improve the overall RMSE.

For each of the 17 DGPs, we computed the coverage of the 90% posterior
credible interval of the ATT. Figure 2 reveals that all methods tended to have
lower than nominal coverage, with coverage decreasing as the amount of
confounding and effect heterogeneity increased. In the presence of
confounding, the use of parametric propensity score estimates appeared to
yield diminished coverage.

Figure 2: Coverage of the 90% posterior credible interval for the overall ATT
for each DGP category using practice-level data. Note the unconfounded
category contains only one DGP while the other four categories contain 4 DGPs.

Figure 3 compares the lengths of the 90% posterior credible intervals produced
by each method relative to CBPS(S). Values larger than one indicate that the
method yielded longer uncertainty intervals.

Figure 3: Lengths of 90% posterior credible interval for the overall ATT
relative to CBPS(S) for each DGP category using practice-level data

Figures 1–3 suggest that in the presence of confounding, fitting BCF with
flexible propensity score estimate (i.e. BART or GBM) might not improve
estimation but could yield wider uncertainty intervals with much better
frequentist coverage than parametric propensity score estimates (i.e. CBPS or
LASSO).

Subgroup effects. We similarly compared the estimation performance and
uncertainty quantification on the eight models for each subgroup effect. For
brevity, we present performance results averaged over all twelve subgroups in
the main text and defer subgroup-by-subgroup comparisons to the Appendix B.
Figure 4 compares the RMSE averaged across all twelve subgroup effects.
Similar to the overall ATT, we did not observe substantial differences in the
quality of estimated subgroup effects across the different methods.

Further, we found that nonparametric propensity score estimators yielded
slightly longer posterior credible intervals that displayed higher frequentist
coverage than parametric propensity score estimators. We observed similar
results for each subgroup; see Figures 9–14 in Appendix B.

Figure 4: RMSE averaged across the twelve subgroup effects for each DGP
category using practice-level data

To summarize, the practice level analyses do not indicate a clear preference
for any of the eight models we compared in terms of point estimation. However,
we did observe that more flexible propensity score estimates (i.e. BART and
GBM) had better-calibrated uncertainty intervals than parametric propensity
score estimates (i.e. CBPS and LASSO). Finally, the use of sparsity-inducing
priors did not seem to affect performance much at all.

###  3.2 Beneficiary-level results

In this section we compare the performance of four of methods fit to the full
beneficiary-level data.

Overall ATT. Figure 5 shows the RMSE for the overall ATT across all
beneficiaries in each category of DGP for the dense and sparse models based on
a CBPS- and GBM-based propensity score estimate.

Figure 5: RMSE for overall ATT for each DGP category using beneficiary-level
data

In the absence of confounding, there is very little discernible difference
between the RMSE for the four models. Similar to findings in the previous
subsection, in the presence of confounding, the GBM-based models have somewhat
lower RMSE as compared to the CBPS based models. The RMSE for the beneficiary
level models are slightly higher when compared to the practice level
counterparts in the presence of confounding, but are essentially the same when
there is no confounding.

Figures 6 and 7 compare the coverage and relative lengths of the 90% posterior
credible intervals for the ATT. As was the case in the practice level models,
using a CBPS-based propensity score estimate yielded shorter intervals with
worse coverage than using a GBM-based propensity score estimate.

Figure 6: Coverage for overall ATT for each DGP category using beneficiary-
level data. Figure 7: Relative interval lengths for overall ATT for each DGP
category using beneficiary-level data.

Interestingly, the uncertainty intervals in the beneficiary level analysis
have much lesser coverage as compared to practice level versions of the same
models. The difference in coverage is quite large in the presence of
confounding.

Subgroup effects. Figure 8 shows boxplots of the RMSE averaged across the 12
subgroup effects for each category of DGP. We found that the RMSE for the
subgroup effects are relatively insensitive to the choice of propensity score
estimate; see also Figures 15 and 16 in Appendix B.2 for comparisons of
performance for each subgroup.

Figure 8: RMSE averaged across all subgroup effects for each DGP category
using beneficiary-level data

Interestingly, we found that when we averaged the RMSEs over the 12 different
subgroup, using the sparsity-inducing prior yielded slightly smaller average
RMSE. Although using the sparsity-inducing prior also tended to produce
somewhat shorter uncertainty intervals, the difference in relative interval
lengths was quite small; see Figures 17–20 in Appendix B.

##  4 Discussion

In this work, we showed that this year’s ACIC Data Challenge admitted a
decomposition similar to the one used by Hahn et al., (2020)’s Bayesian casual
forest model. Because the bcf package could not scale to the size of the data
in the competition, we implemented a faster and less resource-intensive
version, which is available in the flexBCF package. We investigated the
sensitivity of our competition submission to two modeling choices - (i) the
choice of propensity score method, and (ii) enforcing (or not) sparsity in the
BART priors, for both practice-level and beneficiary-level modeling. We found
that the beneficiary-level modeling did not give any more accurate results
than the practice level modeling. Although we did not observe much sensitivity
to our modeling choices on the Challenge datasets, our re-analysis does
suggest that fitting BCF-like models with flexible propensity score estimates
may yield slightly longer uncertainty intervals with better frequentist
coverage.

## Acknowledgements

We are grateful to Mariel Finucane, Jennifer Starling, and Dan Thal from
Mathematica Policy Research for several helpful discussions following the Data
Challenge.

This research was performed using the computing resources and assistance of
the UW–Madison Center For High Throughput Computing (CHTC) in the Department
of Computer Sciences. The CHTC is supported by UW–Madison, the Advanced
Computing Initiative, the Wisconsin Alumni Research Foundation, the Wisconsin
Institutes for Discovery, and the National Science Foundation, and is an
active member of the OSG Consortium, which is supported by the National
Science Foundation and the U.S. Department of Energy’s Office of Science.

Support for A.H.K. and S.K.D. was provided by the University of
Wisconsin–Madison, Office of the Vice Chancellor for Research and Graduate
Education with funding from the Wisconsin Alumni Research Foundation.

## References

  * Abadie, (2005) Abadie, A. (2005).  Semiparametric difference-in-differences estimators.  The Review of Economic Studies, 72(1):1–19.
  * Chipman et al., (2010) Chipman, H. A., George, E. I., and McCulloch, R. E. (2010).  BART: Bayesian additive regression trees.  The Annals of Applied Statistics, 4(1):266–298.
  * Ding et al., (2017) Ding, P., Li, X., and Miratrix, L. W. (2017).  Bridging finite and super population causal inference.  Journal of Causal Inference, 5(2).
  * Ericson, (1969) Ericson, W. A. (1969).  Subjective Bayesian models in sampling finite populations.  Journal of the Royal Statistical Society: Series B, 31(2):195–233.
  * Fong et al., (2022) Fong, C., Ratkovic, M., and Imai, K. (2022).  CBPS: Covariate Balancing Propensity Score.  R package version 0.23.
  * Gelman and Hill, (2006) Gelman, A. and Hill, J. (2006).  Data analysis using regression and multilevel/hierarchical models.  Cambridge University Press.
  * Greenwell et al., (2022) Greenwell, B., Boehmke, B., Cunningham, J., and Developers, G. (2022).  gbm: Generalized Boosted Regression Models.  R package version 2.1.8.1.
  * Hahn et al., (2020) Hahn, P. R., Murray, J. S., and Carvalho, C. M. (2020).  Bayesian regression tree models for causal inference: Regularization, confounding, and heterogeneous effects.  Bayesian Analysis, 15(3):965–1056.
  * Heckman et al., (1997) Heckman, J. J., Ichimura, H., and Todd, P. E. (1997).  Matching as an econometric evaluation estimator: Evidence from evaulating a job training programme.  The Review of Economic Studies, 64(4):605–654.
  * Hernan and Robins, (2020) Hernan, M. A. and Robins, J. M. (2020).  Causal Inference: What If? CRC Press.
  * Hill, (2011) Hill, J. L. (2011).  Bayesian nonparametric modeling for causal inference.  Journal of Computational and Graphical Statistics, 20(1):217–240.
  * Imai and Ratkovic, (2014) Imai, K. and Ratkovic, M. (2014).  Covariate balancing propensity score.  Journal of the Royal Statistical Society: Series B, 76(1):243–263.
  * Lehmann and Romano, (2005) Lehmann, E. L. and Romano, J. P. (2005).  Testing Statistical Hypotheses.  Springer.
  * Linero, (2018) Linero, A. R. (2018).  Bayesian regression trees for high-dimensional prediction and variable selection.  Journal of the American Statistical Association, 113(522):626–636.
  * Rosenbaum, (1984) Rosenbaum, P. R. (1984).  The consequences of adjustment for a concomitant variable that has been affected by the treatment.  Journal of the Royal Statistical Society: Series A, 147(5):656–666.
  * Rubin, (1974) Rubin, D. B. (1974).  Estimating causal effects of treatments in randomized and nonrandomized studies.  Journal of Educational Psychology, 66(5):688–701.
  * Rubin, (1980) Rubin, D. B. (1980).  Randomization analysis of experimental data: The Fisher randomization test.  Journal of the American Statistical Association, 75(371):575–582.
  * Rubin, (1986) Rubin, D. B. (1986).  Comment.  Journal of the American Statistical Association, 81(396):961–962.
  * Särndal, (1978) Särndal, C.-E. (1978).  Design-based and model-based inference in survey sampling.  Scandinavian Journal of Statistics, 5(1):27–52.
  * Schochet, (2013) Schochet, P. Z. (2013).  Estimators for clustered education rcts using the Neyman model for causal inference.  Journal of Educational and Behavioral Statistics, 38(3):219–238.
  * Simon et al., (2011) Simon, N., Friedman, J., Hastie, T., and Tibshirani, R. (2011).  Regularization paths for cox’s proportional hazards model via coordinate descent.  Journal of Statistical Software, 39(5):1–13.
  * Sparapani et al., (2021) Sparapani, R., Spanbauer, C., and McCulloch, R. (2021).  Nonparametric machine learning and efficient computation with Bayesian additive regression trees: The BART R package.  Journal of Statistical Software, 97(1):1–66.
  * Splawa-Neyman et al., (1990) Splawa-Neyman, J., Dabrowska, D. M., and Speed, T. P. (1990).  On the application of probability theory to agricultural experiments. essay on principles. section 9.  Statistical Science, 5(4):465–472.

Appendix

##  Appendix A Identification

By assumption (A3), we have $\delta<\mathbb{P}(Z_{j}=1)<1-\delta$ so that the
population consists of a non-trivial mixture of treated and untreated
practices. Also, there are some values of $\bm{x}$ where
$\delta<\mathbb{P}(Z_{j}=1\mid\bm{X}_{ij}=\bm{x})<1-\delta$.

Let $t$ be a post-intervention time-point (i.e. $t\in\\{3,4\\}$) and let $s$
be a pre-intervention time-point (i.e. $s\in\\{1,2\\}$). Then, for any $s$ and
$t$ where $\delta<\mathbb{P}(Z_{j}=1\mid\bm{X}_{ij}=\bm{x})<1-\delta$, we have
the following:

$\displaystyle\mathbb{E}\left[Y_{ijt}\mid
Z_{j}=1,\bm{X}_{ij}=\bm{x}\right]-\mathbb{E}\left[Y_{ijt}\mid
Z_{j}=0,\bm{X}_{ij}=\bm{x}\right]$
$\displaystyle\,=\mathbb{E}\left[Y_{ijt}(1)\mid
Z_{j}=1,\bm{X}_{ij}=\bm{x}\right]-\mathbb{E}\left[Y_{ijt}(0)\mid
Z_{j}=0,\bm{X}_{ij}=\bm{x}\right]$ $\displaystyle(\text{by A1})$
$\displaystyle\,=\underbrace{\mathbb{E}\left[Y_{ijt}(1)-Y_{ijt}(0)\mid
Z_{j}=1,\bm{X}_{ij}=\bm{x}\right]}_{\operatorname{CATT}\left(\bm{x},t\right)}$
$\displaystyle\,\quad+\mathbb{E}\left[Y_{ijt}(0)\mid
Z_{j}=1,\bm{X}_{ij}=\bm{x}\right]-\mathbb{E}\left[Y_{ijt}(0)\mid
Z_{j}=0,\bm{X}_{ij}=\bm{x}\right]$
$\displaystyle\,=\operatorname{CATT}(\bm{x},t)$
$\displaystyle\,\quad+\mathbb{E}\left[Y_{ijt}(0)-Y_{ijs}(0)\mid
Z_{j}=1,\bm{X}_{ij}=\bm{x}\right]-\mathbb{E}\left[Y_{ijt}(0)-Y_{ijs}(0)\mid
Z_{j}=0,\bm{X}_{ij}=\bm{x}\right]$
$\displaystyle\,\quad+\mathbb{E}\left[Y_{ijs}(0)\mid
Z_{j}=1,\bm{X}_{ij}=\bm{x}\right]-\mathbb{E}\left[Y_{ijs}(0)\mid
Z_{j}=0,\bm{X}_{ij}=\bm{x}\right]$
$\displaystyle\,=\operatorname{CATT}(\bm{x},t)$
$\displaystyle\,\quad+\underbrace{\mathbb{E}\left[Y_{ijt}(0)-Y_{ijs}(0)\mid\bm{X}_{ij}=\bm{x}\right]-\mathbb{E}\left[Y_{ijt}(0)-Y_{ijs}(0)\mid\bm{X}_{ij}=\bm{x}\right]}_{=0}$
$\displaystyle\,\quad+\mathbb{E}\left[Y_{ijs}(0)\mid
Z_{j}=1,\bm{X}_{ij}=\bm{x}\right]-\mathbb{E}[Y_{ijs}(0)\mid
Z_{j}=0,\bm{X}_{ij}=\bm{x}]$ $\displaystyle(\text{by A2})$
$\displaystyle\,=\operatorname{CATT}(\bm{x},t)+\mathbb{E}\left[Y_{ijs}(0)\mid
Z_{j}=1,\bm{X}_{ij}=\bm{x}\right]-\mathbb{E}\left[Y_{ijs}(0)\mid
Z_{j}=0,\bm{X}_{ij}=\bm{x}\right]$
$\displaystyle\,=\operatorname{CATT}(\bm{x},t)+\mathbb{E}\left[Y_{ijs}|Z_{j}=1,\bm{X}_{ij}=\bm{x}\right]-\mathbb{E}\left[Y_{ijs}\mid
Z_{j}=0,\bm{X}_{ij}=\bm{x}\right]$ $\displaystyle(\text{by A1})$

Rearranging we can express the estimand of interest in terms of expectations
of observables as follows:

$\displaystyle\begin{split}\textrm{CATT}(\bm{x},t)\ =\
&\mathbb{E}\left[Y_{ijt}\mid
Z_{j}=1,\bm{X}_{ij}=\bm{x}\right]-\mathbb{E}\left[Y_{ijt}\mid
Z_{j}=0,\bm{X}_{ij}=\bm{x}\right]\\\ &-\Big{\\{}\mathbb{E}\left[Y_{ijs}\mid
Z_{j}=1,\bm{X}_{ij}=\bm{x}\right]-\mathbb{E}\left[Y_{ijs}\mid
Z_{j}=0,\bm{X}_{ij}=\bm{x}\right]\Big{\\}},\end{split}$

which is the same as Equation (3).

##  Appendix B Additional figures

In Section B.1, we present the subgroup effects i.e. CATT’s for both the
practice-level and beneficiary-level models. Since covariates $X_{1},X_{3}$,
and $X_{5}$ are binary, we present their results together, and do the same for
the categorical covariates $X_{2}$ and $X_{4}$. Figures 9 and 10 present the
RMSE for the CATTs defined by binary and categorical covariates respectively;
Figures 11 and 12 show the frequentist coverage of the uncertainty intervals
for the CATTs defined by binary and categorical covariates respectively, while
Figures 13 and 14 show interval lengths relative to the CBPS(S) model.

We follow the same template for the beneficiary-level modeling results in
Section B.2: Figures 15 and 16 show the RMSE for the beneficiary-level model
CATT estimates, Figures 17 and 18 show the frequentist coverage and Figures 19
and 20 present the relative interval lengths of the uncertainty intervals (on
a logarithmic scale) obtained from the beneficiary-level model.

###  B.1 Practice level modeling

Figure 9: RMSE for binary covariates’ CATT for each DGP category using
practice-level data Figure 10: RMSE for categorical covariates’ CATT for each
DGP category using practice-level data Figure 11: Interval coverage for binary
covariates’ CATT for each DGP category using practice-level data Figure 12:
Interval coverage for categorical covariates’ CATT for each DGP category using
practice-level data Figure 13: Relative interval length for binary covariates’
CATT for each DGP category using practice-level data Figure 14: Relative
interval length categorical covariates’ CATT for each DGP category using
practice-level data

###  B.2 Beneficiary-level modeling

Figure 15: RMSE for binary covariates’ CATT for each DGP category using
beneficiary-level data Figure 16: RMSE for categorical covariates’ CATT for
each DGP category using beneficiary-level data Figure 17: Interval coverage
for binary covariates’ CATT for each DGP category using beneficiary-level data
Figure 18: Interval coverage for categorical covariates’ CATT for each DGP
category using beneficiary-level data Figure 19: Relative interval length for
binary covariates’ CATT for each DGP category using beneficiary-level data.
Note: The lengths are displayed on a logarithmic scale Figure 20: Relative
interval length categorical covariates’ CATT for each DGP category using
beneficiary-level data. Note: The lengths are displayed on a logarithmic scale


# EMPIRICAL UPSCALING OF POINT-SCALE SOIL MOISTURE MEASUREMENTS FOR SPATIAL
EVALUATION OF MODEL SIMULATIONS AND SATELLITE RETRIEVALS

###### Abstract

The evaluation of modelled or satellite-derived soil moisture (SM) estimates
is usually dependent on comparisons against in-situ SM measurements. However,
the inherent mismatch in spatial support (i.e., scale) necessitates a cautious
interpretation of point-to-pixel comparisons. The upscaling of the in-situ
measurements to a commensurate resolution to that of the modelled or retrieved
SM will lead to a fairer comparison and statistically more defensible
evaluation. In this study, we presented an upscaling approach that combines
spatiotemporal fusion with machine learning to extrapolate point-scale SM
measurements from 28 in-situ sites to a 100 m resolution for an agricultural
area of 100 km × 100 km. We conducted a four-fold cross-validation, which
consistently demonstrated comparable correlation performance across folds,
ranging from 0.6 to 0.9. The proposed approach was further validated based on
a cross-cluster strategy by using two spatial subsets within the study area,
denoted as cluster A and B, each of which equally comprised of 12 in-situ
sites. The cross-cluster validation underscored the capability of the
upscaling approach to map the spatial variability of SM within areas that were
not covered by in-situ sites, with correlation performance ranging between 0.6
and 0.8. In general, our proposed upscaling approach offers an avenue to
extrapolate point measurements of SM to a spatial scale more akin to climatic
model grids or remotely sensed observations. Future investigations should
delve into a further evaluation of the upscaling approach using independent
data, such as model simulations, satellite retrievals or field campaign data.

Index Terms—  Soil moisture; Upscaling; Machine learning; Extreme Gradient
Boosting; Spatiotemporal fusion; Cross-validation

##  1 Introduction

Soil moisture (SM) is vital in global water, energy cycle and biogeochemical
flux system [1]. Spatiotemporally continuous SM datasets have been extensively
developed and explored throughout remote sensing techniques and model
simulations [2, 3, 4]. However, SM estimates commonly display discrepancies
among different satellite platforms and models [5, 6], even when subjected to
identical meteorological forcings [7, 8]. A thorough evaluation of these
datasets usually depends on comparisons against in-situ SM measurements, while
the inherent mismatch in spatial support (i.e., scale) necessitates a cautious
interpretation of point-to-pixel comparisons.

At the agricultural scale, there is a necessity to enhance our comprehension
and application of in-situ SM measurements. This knowledge is often
accumulated and refined through the deployment of dense SM sensors. Existing
in-situ SM monitoring networks and databases [9, 10, 11] always play a pivotal
role in evaluating SM derived from model simulations or satellite retrievals.
However, the representativeness of point measurements is constrained to a
limited area and their distribution is often uneven [12, 13], leading to gaps
in coverage and potential biases in the validation process. Researchers have
persistently raised questions about the effective alignment of SM footprints
between in-situ data and grids from climatic models or satellite platforms,
with upscaling techniques emerging as an appealing solution [13]. In this
study, we presented an upscaling approach that integrates spatiotemporal
fusion with machine learning to extrapolate point-scale SM measurements from
28 in-situ sites to a 100 m resolution. Subsequently, multiple iterations of
validation were conducted to thoroughly assess the performance of the proposed
approach.

##  2 Study area and data

###  2.1 Study area and in-situ measurements

We selected the Yanco agricultural region as the study area (Fig. 1). Yanco is
situated within the Murrumbidgee Catchment, characterised by a semi-arid
climate that is representative for most agricultural regions in southeast
Australia. The annual total precipitation in this area averages approximately
400 mm [14]. We collected the in-situ SM data from the OzNet Hydrological
Monitoring Network [10] for the specific purposes of upscaling and cross-
validation. The locations of OzNet sites are shown as black squares in Fig. 1.
We categorised the OzNet sites into four groups to facilitate a four-fold
cross-validation. Moreover, we chose two subsets where OzNet sites have a
relatively dense distribution, referred to as the cluster A and cluster B
(Fig. 1). We also employed cluster A and B reciprocally for a cross-cluster
validation. Specifically, sites within cluster A were utilised for training
and predicting within cluster B, and vice versa.

Fig. 1: The location (highlighted in red rectangle within the NSW State) and
the land cover information in 2020 of Yanco agricultural region. The purple
circle, red triangle and black dots represent the locations of in-situ sties
from CosmOz, OzFlux and OzNet, respectively. The cluster A spans the region of
146.06-146.16 °E and 34.62-34.77 °S (10 × 15 km2) and the cluster B spans the
region of 146.25-146.35 °E and 34.92-35.02 °S (10 × 10 km2).

###  2.2 Geospatial predictors for upscaling

We collected a diverse set of geospatial variables to serve as predictors in
the upscaling process (Fig. 2). These predictors encompass indices and surface
temperature data from the MODerate Resolution Imaging Spectroradiometer
(MODIS) [15, 16] and Landsat-resolution products [17, 18], climatic variables
from the ANUClimate 2.0 [19], Smoothed Digital Elevation Model (DEM-S) [20],
and soil data from the Soil and Landscape Grid Australia [21, 22].

##  3 Methodology

Fig. 2 presents the experimental design employed herein, consisting of two key
steps. Firstly, we utilised a nearest neighbour method to resample all
predictors to 0.01° (about 100 m) resolution grids under the World Geodetic
System 1984 (WGS84) datum. Additionally, we performed spatiotemporal fusion
between MODIS and Landsat data to downscale MODIS albedo, NDVI, and LST to a
100 m resolution, which is more comparable to the agricultural foorprints.
Secondly, we conducted machine learning (ML)-based upscaling using in-situ SM
between 2016-2019 as the response, with the collected geospatial data serving
as predictors. To thoroughly assess the performance of upscaling, a four-fold
cross-validation and a cross-cluster validation were carried out, evaluating
the transferability of the upscaling approach while accounting for regional
nuances.

Fig. 2: The experimental design herein. Fig. 3: Density scatterplots of (a)
albedo, (b) NDVI and (c) LST between 01/Jan/2016 and 31/Dec/2019 (excluding
training dates of spatiotemporal fusion). The density of scatters was
calculated using a binning approach.

###  3.1 Spatiotemporal fusion

The Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model (ESTARFM)
[23] has been extensively applied to generate fine-resolution surface
reflectance [23, 24]. It has an exceptional performance in fusing surface
reflectance, while an unbiased variant (ubESTARFM) [25] has demonstrated more
effectiveness in fusing land surface temperature (LST) data. Hence, we
implemented ESTARFM to fuse MODIS and Landsat surface reflectance data (i.e.,
albedo and NDVI), and ubESTARFM to fuse MODIS and Landsat LST data, to
generate daily 100 m estimates of these predictors.

###  3.2 ML-based upscaling

We employed the eXtreme Gradient Boosting (XGBoost) [26] model to perform the
upscaling of in-situ SM. ML models have showcased capability in establishing
robust regression relationships between SM estimates and geospatial predictors
across different scales. However, the majority of ML-based studies focused on
downscaling coarse-resolution SM grids to finer resolutions [27, 28, 29],
while fewer attempts have been made to explore the upscaling space [30, 31].

##  4 Results and discussion

###  4.1 Evaluations of downscaled predictors

Fig. 3 depicts scatterplots illustrating the downscaled performance of (a)
albedo, (b) NDVI, and (c) LST in comparison to MODIS data during 01/Jan/2016 -
31/Dec/2019, encompassing 33,361 samples and using a normalised binning
density. The downscaled albedo exhibited a bias of 0, ubRMSE of 0.03, and R of
0.77, with the majority of samples clustering around the value of 0.2. For the
downscaled NDVI, a slight negative bias of -0.02 was observed, with an ubRMSE
of 0.08 and an R of 0.88. Most NDVI values were concentrated between 0.2 and
0.3. The downscaled LST had a bias of 0.04 K, ubRMSE of 1.25 K and R of 0.99.
The predominant LST concentration occurred at approximately 285 K, with a
secondary concentration between 305 and 315 K.

Fig. 4 provides an illustrative example on 02/Apr/2017 showcasing the spatial
comparison between MODIS data and downscaled predictors. Fig. 4 (g-i) present
downscaled patterns of predictors across the entire study area, demonstrating
consistency with MODIS data but exhibiting sharpened features in comparison
(Fig. 4 a-c). Upon closer examination of a zoomed area (Fig. 4 j-l), the
downscaled predictors revealed enhanced details that are more visually
discernible compared to the MODIS data (Fig. 4 d-f). In general, the
spatiotemporal fusion enabled a better capture of spatial details of
predictors at the field scale, offering improved representation of the
agricultural landscapes.

Fig. 4: Spatial comparison between MODIS and downscaled predictors on
02/Apr/2017. (a-c) are MODIS albedo, NDVI and LST for the study area,
respectively; (d-f) are zoomed windows of MODIS predictors spanning 146.25 to
146.35 °E and -34.95 to -35.05 °N; (g-i) are the downscaled albedo, NDVI and
LST for the study area, respectively; and (j-l) are zoomed windows of
downscaled predictors covering the same area with (d-f).

###  4.2 Evaluations of upscaled in-situ SM

Fig. 5 presents two normalised Taylor diagrams illustrating the performance of
the XGBoost model in upscaling SM using (a) four-fold cross-validation and (b)
cross-cluster validation. The distribution of markers from each fold suggested
variability in performance metrics across the four folds using the XGBoost
model, which consistently demonstrated a correlation performance falling
within the range of 0.6 to 0.9 (Fig. 5 a). Similarly, the distribution of
markers from each cluster also indicated variability in performance metrics,
with the majority of markers demonstrating a correlation within the range of
0.6 to 0.8 (Fig. 5 b). In summary, these diagrams underscore the potential
transferability of the upscaling approach, showcasing its ability to maintain
consistent performance across different folds and spatial clusters.

Fig. 5: Normalised Taylor diagrams of upscaled SM using XGBoost during
2016-2019 for (a) four-fold cross-validation and (b) cross-cluster validation.

Fig. 6 depicts the SHapley Additive exPlanations (SHAP) values of the top six
predictors of the XGBoost model using violin plots from 2016 to 2019 for (a)
all sites; (b) cluster A; and (c) cluster B. For all sites, the order of
predictors was vapour pressure deficit (VPD), albedo, NDVI, DEM, LST, and
Evapotranspiration (ET). In cluster A, the order slightly changed, with NDVI
at the top followed by DEM; while in cluster B, a new predictor, the incoming
solar radiation (Srad), was introduced and was ranked as the 6th. The SHAP
values in all three segments ranged from -0.1 to 0.2, with the majority
concentrated between -0.05 and 0.05. In general, while the top predictors
within the two clusters exhibited similarity, their respective rankings
differed. This discrepancy suggests the model’s reliance on regional land
cover and landscapes in shaping the importance and impact of specific
predictors.

Fig. 6: SHAP values of top 6 predictors of the XGBoost model using violin
plots during 2016-2019 for (a) all sites; (b) cluster A; and (c) cluster B.

Fig. 7 illustrates the spatial distribution of upscaled SM utilising XGBoost
in the study area, with zoomed areas focusing on two clusters at two distinct
dates: 01/Feb/2016 (austral summer) and 30/Jul/2016 (austral winter). Each
cluster was further compared using two different training strategies,
including the global training and cross-cluster training. Notably, the cross-
cluster training strategy revealed more pronounced variations in SM compared
to the global training strategy. This visualisation provides insights into the
effectiveness of different training strategies in upscaling SM to a 100 m
resolution. The observed variations in SM across the landscape underscore the
importance of considerations of regional nuances in training strategies for
accurate and context-specific predictions.

Fig. 7: Spatial examples of upscaled SM using XGBoost on 01/Feb/2016 and
30/Jul/2016 for the Yanco agricultural region, zoomed areas based on global
training, and zoomed areas based on cross-cluster training, respectively.

##  5 Conclusion

In this study, we presented an upscaling approach that integrates
spatiotemporal fusion with machine learning to extrapolate point-scale SM
measurements from 28 OzNet in-situ sites to a 100 m resolution. The
effectiveness of the proposed approach was assessed through multiple
iterations of validation, revealing its capability to map the spatial
variability of SM using sparsely distributed in-situ data, as demonstrated by
the correlation performance ranging between 0.6 and 0.9 and visualised
comparisons. This approach provides a pathway to extrapolate point
measurements of SM to a spatial scale more comparable to climatic model grids
and satellite retrievals. Future investigations should consider incorporating
independent data (e.g., field campaign data) for further assessment.

## References

  * [1] Sonia I. Seneviratne, Thierry Corti, Edouard L. Davin, Martin Hirschi, Eric B. Jaeger, Irene Lehner, Boris Orlowsky, and Adriaan J. Teuling,  “Investigating soil moisture–climate interactions in a changing climate: A review,”  Earth-Science Reviews, vol. 99, no. 3, pp. 125–161, 2010.
  * [2] Gianpaolo Balsamo, Anna Agusti-Panareda, Clement Albergel, Gabriele Arduini, Anton Beljaars, Jean Bidlot, Eleanor Blyth, Nicolas Bousserez, Souhail Boussetta, Andy Brown, Roberto Buizza, Carlo Buontempo, Frédéric Chevallier, Margarita Choulga, Hannah Cloke, Meghan F. Cronin, Mohamed Dahoui, Patricia De Rosnay, Paul A. Dirmeyer, Matthias Drusch, Emanuel Dutra, Michael B. Ek, Pierre Gentine, Helene Hewitt, Sarah P.E. Keeley, Yann Kerr, Sujay Kumar, Cristina Lupu, Jean-François Mahfouf, Joe McNorton, Susanne Mecklenburg, Kristian Mogensen, Joaquín Muñoz-Sabater, Rene Orth, Florence Rabier, Rolf Reichle, Ben Ruston, Florian Pappenberger, Irina Sandu, Sonia I. Seneviratne, Steffen Tietsche, Isabel F. Trigo, Remko Uijlenhoet, Nils Wedi, R. Iestyn Woolway, and Xubin Zeng,  “Satellite and in situ observations for advancing global earth surface modelling: A review,”  Remote Sensing, vol. 10, no. 12, 2018.
  * [3] Jianzhi Dong, Wade T. Crow, Kenneth J. Tobin, Michael H. Cosh, David D. Bosch, Patrick J. Starks, Mark Seyfried, and Chandra Holifield Collins,  “Comparison of microwave remote sensing and land surface modeling for surface soil moisture climatology estimation,”  Remote Sensing of Environment, vol. 242, pp. 111756, 2020.
  * [4] Jian Peng, Clement Albergel, Anna Balenzano, Luca Brocca, Oliver Cartus, Michael H. Cosh, Wade T. Crow, Katarzyna Dabrowska-Zielinska, Simon Dadson, Malcolm W.J. Davidson, Patricia de Rosnay, Wouter Dorigo, Alexander Gruber, Stefan Hagemann, Martin Hirschi, Yann H. Kerr, Francesco Lovergine, Miguel D. Mahecha, Philip Marzahn, Francesco Mattia, Jan Pawel Musial, Swantje Preuschmann, Rolf H. Reichle, Giuseppe Satalino, Martyn Silgram, Peter M. van Bodegom, Niko E.C. Verhoest, Wolfgang Wagner, Jeffrey P. Walker, Urs Wegmüller, and Alexander Loew,  “A roadmap for high-resolution satellite soil moisture applications – confronting product characteristics with user requirements,”  Remote Sensing of Environment, vol. 252, pp. 112162, 2021.
  * [5] Rolf H Reichle, Randal D Koster, Jiarui Dong, and Aaron A Berg,  “Global soil moisture from satellite observations, land surface models, and ground data: Implications for data assimilation,”  Journal of Hydrometeorology, vol. 5, no. 3, pp. 430–442, 2004.
  * [6] Randal D Koster, Zhichang Guo, Rongqian Yang, Paul A Dirmeyer, Kenneth Mitchell, and Michael J Puma,  “On the nature of soil moisture in land surface models,”  Journal of Climate, vol. 22, no. 16, pp. 4322–4335, 2009.
  * [7] Paul A Dirmeyer, Zhichang Guo, and Xiang Gao,  “Comparison, validation, and transferability of eight multiyear global soil wetness products,”  Journal of Hydrometeorology, vol. 5, no. 6, pp. 1011–1033, 2004.
  * [8] Yong Chen and Huiling Yuan,  “Evaluation of nine sub-daily soil moisture model products over china using high-resolution in situ observations,”  Journal of Hydrology, vol. 588, pp. 125054, 2020.
  * [9] Steffen Zacharias, Heye Bogena, Luis Samaniego, Matthias Mauder, Roland Fuß, Thomas Pütz, Mark Frenzel, Mike Schwank, Cornelia Baessler, Klaus Butterbach-Bahl, Oliver Bens, Erik Borg, Achim Brauer, Peter Dietrich, Irena Hajnsek, Gerhard Helle, Ralf Kiese, Harald Kunstmann, Stefan Klotz, Jean Charles Munch, Hans Papen, Eckart Priesack, Hans Peter Schmid, Rainer Steinbrecher, Ulrike Rosenbaum, Georg Teutsch, and Harry Vereecken,  “A network of terrestrial environmental observatories in germany,”  Vadose Zone Journal, vol. 10, no. 3, pp. 955–973, 2011.
  * [10] A. B. Smith, J. P. Walker, A. W. Western, R. I. Young, K. M. Ellett, R. C. Pipunic, R. B. Grayson, L. Siriwardena, F. H. S. Chiew, and H. Richter,  “The murrumbidgee soil moisture monitoring network data set,”  Water Resources Research, vol. 48, no. 7, 2012.
  * [11] Wouter Dorigo, Peter van Oevelen, Wolfgang Wagner, Matthias Drusch, Susanne Mecklenburg, Alan Robock, and Thomas Jackson,  “A new international network for in situ soil moisture data,”  Eos, Transactions American Geophysical Union, vol. 92, no. 17, pp. 141–142, 2011.
  * [12] L. Brocca, R. Morbidelli, F. Melone, and T. Moramarco,  “Soil moisture spatial variability in experimental areas of central italy,”  Journal of Hydrology, vol. 333, no. 2, pp. 356–373, 2007.
  * [13] Wade T. Crow, Aaron A. Berg, Michael H. Cosh, Alexander Loew, Binayak P. Mohanty, Rocco Panciera, Patricia de Rosnay, Dongryeol Ryu, and Jeffrey P. Walker,  “Upscaling sparse ground-based soil moisture observations for the validation of coarse-resolution satellite soil moisture products,”  Reviews of Geophysics, vol. 50, no. 2, 2012.
  * [14] Sabah Sabaghy, Jeffrey P. Walker, Luigi J. Renzullo, Ruzbeh Akbar, Steven Chan, Julian Chaubell, Narendra Das, R. Scott Dunbar, Dara Entekhabi, Anouk Gevaert, Thomas J. Jackson, Alexander Loew, Olivier Merlin, Mahta Moghaddam, Jian Peng, Jinzheng Peng, Jeffrey Piepmeier, Christoph Rüdiger, Vivien Stefan, Xiaoling Wu, Nan Ye, and Simon Yueh,  “Comprehensive analysis of alternative downscaled soil moisture products,”  Remote Sensing of Environment, vol. 239, pp. 111586, 2020.
  * [15] Crystal Schaaf and Zhuosen Wang,  “Mcd43a4 modis/terra+aqua nadir brdf-adjusted reflectance daily l3 global - 500m [dataset],”  NASA LP DAAC.
  * [16] Zhengming Wan,  “New refinements and validation of the collection-6 modis land-surface temperature/emissivity product,”  Remote Sensing of Environment, vol. 140, pp. 36–45, 2014.
  * [17] Fuqin Li, David L. B. Jupp, Shanti Reddy, Leo Lymburner, Norman Mueller, Peter Tan, and Anisul Islam,  “An evaluation of the use of atmospheric and brdf correction to standardize landsat data,”  IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 3, no. 3, pp. 257–270, 2010.
  * [18] Juan P. Guerschman, Tim R. McVicar, Jamie Vleeshower, Thomas G. Van Niel, Jorge L. Peña-Arancibia, and Yun Chen,  “Estimating actual evapotranspiration at field-to-continent scales by calibrating the cmrset algorithm with modis, viirs, landsat and sentinel-2 data,”  Journal of Hydrology, vol. 605, pp. 127318, 2022.
  * [19] M. F. Hutchinson, T. Xu, J. L. Kesteven, I. J. Marang, and B. J. Evans,  “Anuclimate v2.0 [dataset],”  NCI Australia, Canberra, 2021.
  * [20] GA,  “Digital elevation model (dem) of australia derived from lidar 5 metre grid,”  Geoscience Australia, Canberra, 2015.
  * [21] Brendan Malone and Ross Searle,  “Updating the australian digital soil texture mapping (part 1): re-calibration of field soil texture class centroids and description of a field soil texture conversion algorithm,”  Soil Research, vol. 59, pp. 419–434, 2021.
  * [22] Brendan Malone and Ross Searle,  “Updating the australian digital soil texture mapping (part 2): spatial modelling of merged field and lab measurements,”  Soil Research, vol. 59, pp. 435–451, 2021.
  * [23] Xiaolin Zhu, Jin Chen, Feng Gao, Xuehong Chen, and Jeffrey G. Masek,  “An enhanced spatial and temporal adaptive reflectance fusion model for complex heterogeneous regions,”  Remote Sensing of Environment, vol. 114, no. 11, pp. 2610–2623, 2010.
  * [24] Irina V. Emelyanova, Tim R. McVicar, Thomas G. Van Niel, Ling Tao Li, and Albert I.J.M. van Dijk,  “Assessing the accuracy of blending landsat–modis surface reflectances in two landscapes with contrasting spatial and temporal dynamics: A framework for algorithm selection,”  Remote Sensing of Environment, vol. 133, pp. 193–209, 2013.
  * [25] Yi Yu, Luigi J. Renzullo, Tim R. McVicar, Brendan P. Malone, and Siyuan Tian,  “Generating daily 100 m resolution land surface temperature estimates continentally using an unbiased spatiotemporal fusion approach,”  Remote Sensing of Environment, vol. 297, pp. 113784, 2023.
  * [26] Jerome H. Friedman,  “Greedy function approximation: A gradient boosting machine,”  The Annals of Statistics, vol. 29, no. 5, pp. 1189–1232, 2001.
  * [27] Jungho Im, Seonyoung Park, Jinyoung Rhee, Jongjin Baik, and Minha Choi,  “Downscaling of amsr-e soil moisture with modis products using machine learning approaches,”  Environmental Earth Sciences, vol. 75, pp. 1–19, 2016.
  * [28] Ahmed Samir Abowarda, Liangliang Bai, Caijin Zhang, Di Long, Xueying Li, Qi Huang, and Zhangli Sun,  “Generating surface soil moisture at 30 m spatial resolution using both data fusion and machine learning toward better water resources management at the field scale,”  Remote Sensing of Environment, vol. 255, pp. 112301, 2021.
  * [29] Yi Yu, Luigi J. Renzullo, and Siyuan Tian,  “Continental scale downscaling of awra-l analysed soil moisture using random forest regression,”  MODSIM2021, 24th International Congress on Modelling and Simulation, Sydney, Australia, 5-10 December, 2021.
  * [30] Daniel Clewley, Jane B. Whitcomb, Ruzbeh Akbar, Agnelo R. Silva, Aaron Berg, Justin R. Adams, Todd Caldwell, Dara Entekhabi, and Mahta Moghaddam,  “A method for upscaling in situ soil moisture measurements to satellite footprint scale using random forests,”  IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 10, no. 6, pp. 2663–2673, 2017.
  * [31] Dongying Zhang, Wen Zhang, Wei Huang, Zhiming Hong, and Lingkui Meng,  “Upscaling of surface soil moisture using a deep learning model with viirs rdr,”  ISPRS International Journal of Geo-Information, vol. 6, no. 5, 2017.


#  A Tsallis-like effective exponential delay discounting model and its
implications

Trambak Bhattacharyya<EMAIL_ADDRESS>(corresponding author)
Bogoliubov Laboratory of Theoretical Physics, Joint Institute for Nuclear
Research, Dubna, 141980,
Moscow Region, Russia  Shanu Shukla<EMAIL_ADDRESS>Indian Institute
of Management Indore, Rau-Pithampur Road, Indore, Madhya Pradesh, India,
453556  Ranu Pandey<EMAIL_ADDRESS>JSS Academy of Technical
Education, Noida, India, 201301

###### Abstract

We propose an effective exponential model of delay discounting considering
fluctuation in impulsivity. This model is seen to be dual to the two-parameter
Tsallis model of delay discounting proposed by Takahashi in 2007 (Ref.
takahashi1 ). We demonstrate that the parameters of the model can be related
to the mean and the variance of impulsivity. Our findings provide an intuitive
way to explore the origin of the Tsallis distribution in a social system.

Keywords: Tsallis statistics, delay discounting model, impulsivity fluctuation

##  I Introduction

Delay discounting is associated with the evolution of the value of a reward
with delay to its receipt. For example, money of an amount $x$ received
immediately is more ‘valuable’ than the amount received after a delay time
$t$. However, when some features about a situation remain unspecified,
uncertainties in decision arise and making a choice becomes more difficult.
One such example is choosing between USD 50 right now and USD 70 in a year. In
this case, one has to choose between smaller-sooner rewards and larger-later
rewards. One example of choosing a smaller-sooner reward is a student
unprepared for an exam going to a party the night before the test. The
immediate reward brought out of this action appears to be more desirable than
a delayed reward (which may be of greater value in terms of career, earning
and so on). Choosing the smaller-sooner reward is termed as an impulsive
action by the behavioral analysts ikuk . This implies recognizing more distant
rewards as of lesser value. For a variety of species, reward types and sample,
the value of a reward is discounted over time mazur ; bdm ; lempert ; west ;
cajueiro . However, there is more to the story.

In behavioral studies, the rationality of agents has long been a matter of
discussion and debate kahneman . In the studies of discounting, rational
behavior of an agent is characeterized by inter-temporal consistency which may
be understood from the following example. Let us assume that two sets of a
pair of options are given to participants. Options in set 1 are A1) USD 100
available after 365 days, and B1) USD 120 available after 372 days. Options in
set 2 are A2) USD 100 available right now, and B2) USD 120 available after 7
days. Choices A1 and A2 are more impulsive than B1 and B2 respectively because
these choices manifest a preference for smaller-sooner rewards. However,
things become interesting when a participant is asked to choose one option
from both sets during a single experiment. Many people may opt for the
smaller-sooner reward in set 1, whereas for set 2 they settle in for a larger-
later reward (despite the waiting time in both cases are 7 days). This
inconsistency in degree of impulsivity is known as ‘irrationality’. However,
if in both cases one takes impulsive decisions, that is still rational, and
the same is true for those who consistently opt for larger-later rewards in
both cases.

In the case of rational agents, the exponential model of delay discounting is
seen to be appropriate. The exponential delay discounting model considers an
exponential dependence of the subjective value of a reward on delay given by,

$V(D)=V(0)\exp({-\kappa D}),$ (1)

where $V(D)$ is the discounted value of a reward after a delay $D$,
$V(0)>V(D)$ is the undiscounted value of the reward, and $\kappa$ is the
impulsivity parameter.

Exponential model is related to the situation when the discount rate is
constant over time. However, things change when this is not the case. In that
case, the exponential model becomes unable to explain inter-temporal
inconsistency described in a previous paragraph. Several other models like the
hyperbolic model rachilinhyp  have been proposed to consider this phenomenon
in the behavioral studies. But, the hyperbolic model cannot differentiate
between inconsistency and impulsivity. This is the reason why Takahashi
takahashi1  proposed a delay discounting model based on the Tsallis
$q$-exponential function which was also used by Cajueiro cajueiro  in a study
on inter-temporal choices.

Unlike the exponential model, the Tsallis model is described by two
parameters, $q$ (also called the Tsallis $q$ parameter) and impulsivity
$\kappa$. A fundamental approach involving the mathematical techniques of
complex Hilbert space used in quantum mechanics has been taken by Yukalov and
Sornette yukalov . By considering a time-dependent decay rate in the evolution
of a ‘prospect set’, they obtain a power-law time discounting function which
is similar to the Tsallis function. Hence the Tsallis function comes both as a
necessity to describe the observed phenomenon and as a consequence of a very
general mathematical ritual. Because of these reasons, the Tsallis
$q$-exponential function (henceforth to be called the Tsallis function) has
been used in many behavioral studies over the years takahashi2 ; takahashiqdt
; lempert ; fpsy1 ; fph1 ; ejmbe .

As indicated in one of them takahashi2 , the Tsallis function has its origin
in physics. This function can be obtained from the Tsallis entropy, which is a
generalization of the Boltzmann-Gibbs entropy. While maximization of the
Boltzmann-Gibbs entropy gives rise to the exponential function, maximization
of the Tsallis entropy yields a Tsallis $q$-exponential function. This
approach of getting the Tsallis function may be termed as the ‘thermodynamic
approach’ which is a standard procedure of obtaining the equilibrium (maximum
entropy state, to be understood as the state with the maximum disorderliness)
thermodynamic quantities in systems considered in physics.

However, in the present article, we do not consider the approach involving
maximization in a social system (see for example cajueiro ). Here we rather
try to examine why the Tsallis function may appear in a social system from a
perspective that is different from the ones presented in takahashi2 ; yukalov
. The approach taken by us is known as the so-called ‘super-statistics
approach’ beck ; wilkprl  which helps one build another intuitive
understanding of the model. This approach is based on the fact that social and
physical systems contain fluctuations. Let us take the example of the
temperature of a room. Even though we say that the ‘room temperature’ is
30${}^{\text{o}}$C (say), this may not reflect the reality that the
temperature in one corner of the room may be slightly below (say
29.8${}^{\text{o}}$C), and the other corner may be slightly above (say
31.2${}^{\text{o}}$C) than 30${}^{\text{o}}$. So, where does this number 30
come from? This may be the average temperature of the room. If, for example,
we have $n_{1}$ different points in the room having temperature $T_{1}$,
$n_{2}$ different points in the room having temperature $T_{2}$ and so on, the
average temperature of the room is defined by
$\left(n_{1}T_{1}+n_{2}T_{2}+....\right)/(n_{1}+n_{2}+...)$. In social systems
too we encounter fluctuating quantities that follow some distributions, and
the quantities like mean, variance become of our interest.

The novelty of the present work lies in the fact that from this simple
consideration of a fluctuating system, the parameters used in the Tsallis
model may be given some physical meaning. Although we do not deny that the
Tsallis parameter (to be denoted as $q$) is already identified with the
‘measure of inconsistency’ cajueiro ; takahashi1 , our approach provides
another intuitive way of looking at the model parameters. With this general
consideration, we propose a model called the effective exponential model which
contains a Tsallis function and has some correspondence with the existing
Tsallis models.

The paper is organized in the following way. In the next section, we discuss
the effective exponential model in detail. Section III contains a detailed
discussion of our findings. Next, we summarize and conclude in Sec. IV.

##  II Description of the model and its application

###  II.1 Fluctuating impulsivity: how is it distributed?

We consider a social system in which impulsivity $\kappa$ is a positive,
random variable. But we assume that we do not know anything about how they are
distributed. In a typical statistical problem, we have a random variable whose
distribution is unknown. In that case, the gamma distribution serves as the
least-informative option hogg . Let us assume that the impulsivity
distribution of a sample is given by the normalized (integrates to 1) gamma
distribution,

$f(\kappa)=\frac{1}{\Gamma\left(\frac{1}{q-1}\right)}\left\\{\frac{1}{(q-1)\kappa_{0}}\right\\}\left\\{\frac{\kappa}{(q-1)\kappa_{0}}\right\\}^{\frac{1}{q-1}-1}\exp\left(-\frac{\kappa}{(q-1)\kappa_{0}}\right),$
(2)

where $q>1$ and $\kappa_{0}>0$ (so that $f(\kappa)$ is a positive real number)
are two parameters, and the $\Gamma$ function is represented by the following
integral when $q>1$,

$\Gamma\left(\frac{1}{q-1}\right)=\int_{0}^{\infty}\exp(-x)x^{\frac{1}{q-1}-1}dx.$
(3)

For different $q$ and $\kappa_{0}$ values the gamma distribution in Eq. (2)
give rise to different shapes as can be seen in Figs. 2 and 2.

Figure 1: Gamma distributions for a fixed $\kappa_{0}$ and different $q$
values.

Figure 2: Gamma distributions for a fixed $q$ and different $\kappa_{0}$
values.

###  II.2 Impulsivity fluctuation in the exponential model and an effective
factor

We consider the exponential model of delay discounting given by Eq. (1) and
consider that impulsivity $\kappa$ is a fluctuating quantity. Considering that
$\kappa$ in Eq. (1) is distributed according to Eq. (2), we can define an
‘effective exponential factor’ (essentially an average) defined by,

$\displaystyle\epsilon_{\text{eff}}$ $\displaystyle=$
$\displaystyle\frac{\int_{0}^{\infty}d\kappa\exp({-\kappa
D})f(\kappa)}{\int_{0}^{\infty}d\kappa f(\kappa)}$ (4) $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\kappa\exp({-\kappa D})f(\kappa)$

in which each $\kappa$ in the exponential is now weighted by the normalized
distribution $f(\kappa)$.

The factor $\epsilon_{\text{eff}}$ accommodates fluctuating $\kappa$ values
(which may range from very small to very large values) already present in a
social system. In Eq. (1), we replace the exponential function with
$\epsilon_{\text{eff}}$, and propose the effective exponential model given by,

$V(D)=V(0)\epsilon_{\text{eff}}$ (5)

Defining an effective exponential factor is the reason why the proposed model
is named the ‘effective exponential model’. This approach is used in physical
systems wilkprl  where temperature is a fluctuating quantity and an effective
Boltzmann factor is defined (Boltzmann is exponential). We show that this
effective exponential factor is nothing but the Tsallis function. An
analytical closed form of $\epsilon_{\text{eff}}$ is obtained as described
below.

From Eq. (2) we have,

$\displaystyle\epsilon_{\text{eff}}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\kappa\exp({-\kappa D})f(\kappa)$ (6)
$\displaystyle=$
$\displaystyle\frac{1}{\Gamma\left(\frac{1}{q-1}\right)}\left\\{\frac{1}{(q-1)\kappa_{0}}\right\\}\int_{0}^{\infty}d\kappa\exp({-\kappa
D})\left\\{\frac{\kappa}{(q-1)\kappa_{0}}\right\\}^{\frac{1}{q-1}-1}\exp\left(-\frac{\kappa}{(q-1)\kappa_{0}}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{\Gamma\left(\frac{1}{q-1}\right)}\left\\{\frac{1}{(q-1)\kappa_{0}}\right\\}^{\frac{1}{q-1}}\int_{0}^{\infty}d\kappa\exp\left[{-\kappa\left(D+\frac{1}{(q-1)\kappa_{0}}\right)}\right]\kappa^{\frac{1}{q-1}-1}$

Now, defining
$\kappa\left(D+\frac{1}{(q-1)\kappa_{0}}\right)=\kappa^{\prime}$, the above
integral can be written as,

$\displaystyle\epsilon_{\text{eff}}$ $\displaystyle=$
$\displaystyle\frac{1}{\Gamma\left(\frac{1}{q-1}\right)}\left\\{\frac{1}{(q-1)\kappa_{0}}\right\\}^{\frac{1}{q-1}}\left(D+\frac{1}{(q-1)\kappa_{0}}\right)^{-\frac{1}{q-1}}\int_{0}^{\infty}d\kappa^{\prime}\exp(-\kappa^{\prime})\kappa^{\prime\frac{1}{q-1}-1}.$
(7)

From Eq. (3), the integral in the above equation is nothing but
$\Gamma\left(\frac{1}{q-1}\right)$ and hence,

$\epsilon_{\text{eff}}=\left(1+(q-1)\kappa_{0}D\right)^{\frac{1}{1-q}}=\exp_{q}(-\kappa_{0}D),$
(8)

where $\exp_{q}$ is the $q$-exponential function Tsal88  which becomes the
conventional exponential function in the limit $q\rightarrow 1$. Hence, a
modified form of Eq. (1) is given by,

$\displaystyle V(D)$ $\displaystyle=$ $\displaystyle
V(0)\exp_{q}(-\kappa_{0}D)$ (9) $\displaystyle=$
$\displaystyle\frac{V(0)}{\left(1+(q-1)\kappa_{0}D\right)^{\frac{1}{q-1}}}.~{}~{}~{}~{}~{}~{}\text{(--EEM--)}$

This is the main result of our paper which describes an ‘effective exponential
model (EEM)’ of delay discounting. Takahashi takahashi1  has also proposed a
similar Tsallis model. However, there are still some differences among the EEM
and the other Tsallis model. We will discuss that in Section III.

###  II.3 Significance of the parameters $q$ and $\kappa_{0}$

Now, let us investigate what the parameters $\kappa_{0}$ and $q$ signify. As
shown next, parameter $\kappa_{0}$ is the mean impulsivity of the system, and
$q$ is related to the variance of the impulsivity distribution. From Eq. (2),

$\displaystyle\langle\kappa\rangle$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\kappa~{}\kappa f(\kappa)$ (10)
$\displaystyle=$
$\displaystyle\frac{1}{\Gamma\left(\frac{1}{q-1}\right)}\left\\{\frac{1}{(q-1)\kappa_{0}}\right\\}^{\frac{1}{q-1}}\int_{0}^{\infty}d\kappa~{}\kappa^{\frac{1}{q-1}}\exp\left(-\frac{\kappa}{(q-1)\kappa_{0}}\right)$

By defining $k/(q-1)\kappa_{0}=\kappa^{\prime\prime}$ and using Eq. (3) we
obtain,

$\displaystyle\langle\kappa\rangle$ $\displaystyle=$
$\displaystyle\frac{\Gamma\left(\frac{1}{q-1}+1\right)}{\Gamma\left(\frac{1}{q-1}\right)}\times(q-1)\kappa_{0}$
(11) $\displaystyle=$
$\displaystyle\frac{\Gamma\left(\frac{1}{q-1}\right)}{(q-1)\Gamma\left(\frac{1}{q-1}\right)}\times(q-1)\kappa_{0}$
$\displaystyle=$ $\displaystyle\kappa_{0}$

Now, to find out the variance, we need to find out the mean square
impulsivity. Again, using Eq. (2) and redefining the integration variable
appropriately we obtain,

$\displaystyle\langle\kappa^{2}\rangle$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\kappa~{}\kappa^{2}f(\kappa)$ (12)
$\displaystyle=$
$\displaystyle\frac{1}{\Gamma\left(\frac{1}{q-1}\right)}\left\\{\frac{1}{(q-1)\kappa_{0}}\right\\}^{\frac{1}{q-1}}\int_{0}^{\infty}d\kappa~{}\kappa^{\frac{1}{q-1}+1}\exp\left(-\frac{\kappa}{(q-1)\kappa_{0}}\right)$
$\displaystyle=$
$\displaystyle\frac{\Gamma\left(\frac{1}{q-1}+2\right)}{\Gamma\left(\frac{1}{q-1}\right)}\times(q-1)^{2}\kappa_{0}^{2}$
$\displaystyle=$ $\displaystyle q\kappa_{0}^{2}$

Hence,

$\displaystyle\frac{\langle\kappa^{2}\rangle-\langle\kappa\rangle^{2}}{\langle\kappa\rangle^{2}}=q-1;~{}~{}\Rightarrow
q=1+\frac{\langle\kappa^{2}\rangle-\langle\kappa\rangle^{2}}{\langle\kappa\rangle^{2}}$
(13)

Eqs. (11) and (13) show that the $\kappa_{0}$ parameter is the mean
impulsivity of the sample, and the $q$ parameter is related to the relative
variance of impulsivity of the sample. From Eq. (13) it is evident that $q>1$.
Hence, the $q\rightarrow 1$ limit signifies a very narrow impulsivity
distribution with negligible fluctuation. A similar relation for physical
systems can be found in Ref. wilkprl  which involves temperature, and Ref.
wilkprc  which involves particle number.

##  III Discussion

###  III.1 Facets of impulsivity

In this article, we consider fluctuations present in a system and we argue
that the gamma distribution given by $f(\kappa)$ is the most general
distribution of a fluctuating quantity. In that sense, our approach is based
on a condition (fluctuation) that is ubiquitous in any system. In a social
system, fluctuating impulsivity is a common feature. However, researchers
always distinguish between two kinds of impulsivity, trait and state, also
called ‘facets’, denoted by $\kappa_{\text{T}}$ and $\kappa_{\text{S}}$
respectively. A trait characteristic is inherent, and a person (or a
representative of any other species) is born with it. On the other hand, the
state characteristic develops over time and varies. Hence, when one talks
about fluctuations, it may seem that she implies the state characteristics
only. However, in a study traitstate  it has been found that that university
soccer players with higher trait anxiety scores may experience increased state
anxiety under pressure. Though it discusses a connection between trait and
state anxiety, yet it is possible that the trait and state impulsivity are
also connected. Hence we consider impulsivity obtained from delay discounting
data as a combination of both. That is because a more or less unaltered
quantity combined with an alternating quantity gives rise to another
alternating quantity.

###  III.2 Drift and diffusion of impulsivity

It is interesting to note in Figs. 2 and 2 that when we vary $\kappa_{0}$ and
$q$ of the distribution, the peak and the width of the distribution changes.
This shift in the peak of average impulsivity and width is analogous to the
velocity-space drift and diffusion of a collection of particles traversing
inside a medium. In this case, the average velocity and the width change
because of the interaction with the medium molecules. We can similarly think
about the change in the impulsivity distribution in terms of drift and
diffusion in the impulsivity space. If we study the time evolution of the
impulsivity distribution $f(\kappa)$ of a test sample, it can schematically be
written in terms of two terms representing the drift ($\gamma_{\kappa}>0$) and
diffusion ($\delta_{\kappa}>0$) in the following way,

$\displaystyle\dot{f}(\kappa)=-\text{change~{}due~{}to~{}}\gamma_{\kappa}+\text{change~{}due~{}to~{}}\delta_{\kappa},$
(14)

where the dot on $f(\kappa)$ signifies total time derivative. Time evolution
of the impulsivity distribution of a sample may be a consequence of many
factors like changing social status, social interaction and so on.

###  III.3 Comparison of the EEM with another Tsallis-like model

Another Tsallis-like model proposed by Takahashi in Ref. takahashi1  has the
following form according to our notations,

$\displaystyle V(D)$ $\displaystyle=$
$\displaystyle\frac{V(0)}{\exp_{q}\left({\kappa_{0}D}\right)}$ (15)
$\displaystyle=$
$\displaystyle\frac{V(0)}{(1+(1-q)\kappa_{0}D)^{\frac{1}{1-q}}}.~{}~{}~{}~{}~{}~{}\text{(--TM--)}$

Eq. (15) is to be referred to as the Takahashi model (TM) from now on. If $q$
is replaced by $2-q$ in the TM, one obtains the EEM (Eq. 9), and this is the
reason why the TM and EEM may be termed as ‘dual to each other’. This duality
reminds one of the same correspondence existing between the $q$-exponential
function and its inverse eqinv . The EEM relates $q$ to the relative variance
of the impulsivity distribution, and it is always positive. However, the $q$
values in the TM model can assume negative values. In spite of the difference
in their appearances, it is seen that numerically these two models are
equivalent. However, the parameter values in two approaches have a
$q\leftrightarrow 2-q$ correspondence because of the dual nature.

There is another subtle difference between the two Tsallis-inspired models
which we discuss now. This discussion is based on the fact that we demand $V$
to be a real-valued quantity and not a complex number. From Eq. (9) it is
clear that if $1+(q-1)\kappa_{0}D$ is negative, one may get a complex valued
quantity because the $q$ parameter (also called the Tsallis parameter or the
entropic index) is a real number. For example, for $q=0.1$, $\kappa_{0}=0.2$
and $D=10$, $1+(q-1)\kappa_{0}D=-0.8$ and
$\left(1+(q-1)\kappa_{0}D\right)^{1/(1-q)}$ is a complex number. However,
since we have $q>1,~{}\kappa_{0}>0$ and $D>0$ in our model (Eq. 9), we do not
encounter complex values. The same constraint of positivity of the argument
leads to an upper-bound of the $q$ values in the TM,

$q<1+\frac{1}{\kappa_{0}D}.$ (16)

###  III.4 Different limits of the EEM

Now, let us discuss different limits of the effective exponential model. In
the limit $q\rightarrow 1$ (vanishing fluctuation), Eq. (9) gives back the
exponential model of delay discounting as shown below,

$\displaystyle\lim_{q\rightarrow 1}V(D)$ $\displaystyle=$
$\displaystyle\lim_{q\rightarrow
1}V(0)\left(1+(q-1)\kappa_{0}D\right)^{-\frac{1}{q-1}}$ (17) $\displaystyle=$
$\displaystyle\lim_{q\rightarrow
1}V(0)\exp\left({-\frac{1}{q-1}\ln\left(1+(q-1)\kappa_{0}D\right)}\right)$
$\displaystyle=$ $\displaystyle\lim_{q\rightarrow
1}V(0)\exp\left({-\frac{1}{q-1}\left\\{(q-1)\kappa_{0}D-\mathcal{O}(q-1)^{2}+...\right\\}}\right)$
$\displaystyle=$ $\displaystyle V(0)\exp\left(-\kappa_{0}D\right)$

When we take the $\kappa_{0}D\rightarrow 0$ limit, Eq. (17) gives back the
hyperbolic equation given by,

$\lim_{\kappa_{0}D\rightarrow 0}V(D)\approx\frac{V(0)}{1+\kappa_{0}D}.$ (18)

Also, when $q=2$, Eq. (9) gives back the hyperbolic model. Hence the EEM
generalizes the existing models of delay discounting which may be obtained by
taking appropriate limits of the EEM.

##  IV Summary, conclusions, and outlook

To summarize, considering impulsivity to be a fluctuating quantity we describe
an effective exponential model of delay discounting phenomenon. This model
described in the article is found to be dual to the Tsallis model of delay
discounting proposed by Takahashi takahashi1 . In that sense, our findings
support the existing Tsallis models which efficiently describe different
aspects of delay discounting (like the inter-temporal inconsistency) which are
important in deciding policies in many sectors like finance, medical etc.

The novelty of our article lies in the approach it takes. We consider a very
general characteristic of a system (fluctuation) and propose a modification of
the existing exponential model to arrive at another intuitive understanding of
the Tsallis model. With the help of this approach we explain the significance
of the model parameters. We find that

  * •

the Tsallis parameter $q$ is related to the relative variance (variance
divided by average square) of the impulsivity distribution and is greater than
1.

  * •

the impulsivity parameter $\kappa_{0}$ is actually the mean of the impulsivity
distribution.

It will be interesting to apply this model for participants from varied age
groups and social backgrounds to study how impulsivity changes. It will also
be interesting to see how the personality of an individual affects this
behavior. From the observations, one may parameterize the evolution the $q$
and $\kappa_{0}$ parameters with age to predict the impulsivity behavior of
people from certain age group from a certain socio-cultural background. These
findings will be important in deciding policies in relevant areas.

Our study can throw light upon some other interesting aspects of social
science in general. We try to elaborate on this in the next few sentences. As
mentioned in Section III, we consider impulsivity obtained from delay
discounting data to be a combination of trait and state facets. However, how
to combine these two facets to get the combination is another interesting
question which future researchers may revisit. In many studies it is assumed
to be a sum score given by the addition of facets. Since we consider
impulsivity to be bi-faceted, the total impulsivity may be given by the sum
score of $\kappa_{\text{T}}$ and $\kappa_{\text{S}}$. In that case, one
immediately finds that addition is only one of the options among a large
number of operators which combine these two numbers. There can also be some
other ways to combine the facets. One may propose a gedankenexperiment
gedanken  in which the trait and state impulsivity of a group of participants
are quantified from questionnaires, and the combined impulsivity of the same
participants are measured by fitting the data obtained from the delay-
discounting experiments using the Tsallis model. From these three sets of
numbers, one can try to find out a rule of combining them. This question
becomes really important when one needs to score competitors for providing
grants, or for providing some aid and in many other situations in which a
questionable method of evaluation may be disappointing (for the marginalized
section, for example).

We would like to point out that the superstatistics approach may also be taken
in describing financial market data using the $q$-gaussian distribution
tsfinnewjphys ; tsfinphysrept  that is obtained by replacing the exponential
function in Eq. (4) by a gaussian. In Refs. tsfinnewjphys ; tsfinphysrept ,
authors observe that $q$ approaches unity when the time lag increases for the
financial return distribution. This can be explained if $q$ can be identified
as an internal variable that follows a relaxation-like equation qint ,

$\frac{dq}{d\delta t}=-\frac{q-q_{\text{eq}}}{\tau},$ (19)

where $q_{\text{eq}}=1$, $\delta t$ is the time-lag, and $\tau$ is a
characteristic time scale. The solution of the above equation yields $q(\delta
t)=1+(q_{\text{eq}}-1)e^{-\frac{\delta t}{\tau}}$, and supports the
observation.

To conclude, the effective exponential model lends support to the Tsallis
model of delay discounting. However, its implications may go way beyond that
as discussed in the article. We have provided some insight into the model, and
have raised some questions that seemed to be relevant for a deeper
understanding. With this input, we hope that future studies will shine more
light on this very interesting phenomenon of delay discounting.

## References

  * (1) A.L. Odum, Delay Discounting: I’m a k, You’re a k, J. Exp. Anal. Behav. 96 (3) (2011) 427-439
  * (2) J.E. Mazur, An adjusting procedure for studying delayed reinforcement, in: M.L. Commons, J.E. Mazur, J.A. Nevin, H. Rachlin (Eds.), Quantitative Analysis of Behavior: Vol. 5. The Effect of Delay and of Intervening Events on Reinforcement Value, Hillsdale, NJ: Erlbaum, 1987, pp. 55-73
  * (3) H. Rachlin and B. A. Jones, Social Discounting and Delay Discounting, J. Behav. Dec. Making 21 (2008) 29-43
  * (4) K.M. Lempert, and D.A. Pizzagalli, Delay discounting and future-directed thinking in anhedonic individuals, J. Behav. Ther. & Exp. Psychiat. 41 (2010) 258-264
  * (5) B.J. West and P. Grigolini, A psychophysical model of decision making,Physica A 389 (2010) 3580-3587
  * (6) D.O. Cajueiro, A note on the relevance of the q-exponential function in the context of intertemporal choices, Physica A 364 (2006) 385-388
  * (7) A. Tversky and D. Kahneman, Rational Choice and the Framing of Decisions, The Journal of Business Vol. 59, No. 4, Part 2: The Behavioral Foundations of Economic Theory (Oct., 1986), pp. S251-S278 (28 pages)
  * (8) H. Rachlin, Notes on discounting. Journal of The Experimental Analysis of Behavior, 85 (2006) 425-435.
  * (9) T. Takahashi, A comparison of intertemporal choices for oneself versus someone else based on Tsallis’ statistics, Physica A 385 (2007) 637-644
  * (10) V.I. Yukalov and D. Sornette, Physics of uncertainty in Quantum Decision Making, The European Physical Journal B, 71 No. 4 (2009) pp. 533-548.
  * (11) T. Takahashi, A social discounting model based on Tsallis’ statistics, Physica A 389 (2010) 3600-3603
  * (12) T. Takahashi, H. Nishinaka, T. Makino, R. Han and H. Fukui, An Experimental Comparison of Quantum Decision Theoretical Models of Intertemporal Choice for Gain and Loss, Journal of Quantum Information Science 2 (2012) 119-122
  * (13) K. Ishii and C. Eisen, Cultural Similarities and Differences in Social Discounting: The Mediating Role of Harmony-Seeking, Frontiers in Psychology 9 (2018) 1426
  * (14) S.E. Stegall, T. Collette, T. Kinjo, T. Takahashi and P. Romanowich, Quantitative Cross-Cultural Similarities and Differences in Social Discounting for Gains and Losses, Frontiers in Public Health 7 (2019) 297
  * (15) I.M.P. Oller, S.C. Rambaud and M.d.C.V. Martinez, Discount models in intertemporal choice: an empirical analysis, European Journal of Management and Business Economics, 30 No. 1 (2021) pp. 72-91
  * (16) C.Beck and E.G.D.Cohen, Superstatistics, Physica A: Statistical Mechanics and its Applications, Volume 322 (2003) 267-275
  * (17) G. Wilk, and Z. Włodarczyk, Interpretation of the nonextensivity parameter q in some applications of tsallis statistics and lèvy distributions, Phys. Rev. Lett. 84 (2000) 2770-2773
  * (18) R.V. Hogg, J.W. McKean, A.T. Craig, Introduction to Mathematical Statistics, Pearson Education, Boston (2019)
  * (19) C. Tsallis, Possible generalization of Boltzmann-Gibbs statistics, J. Stat. Phys. 52 (1988) 479-487
  * (20) G. Wilk, and Z. Włodarczyk, Multiplicity fluctuations due to the temperature fluctuations in high-energy nuclear collisions, Phys. Rev. C 79 (2009) 054903
  * (21) M. Horikawa, A. Yagi, The Relationships among Trait Anxiety, State Anxiety and the Goal Performance of Penalty Shoot-Out by University Soccer Players, PLOS ONE, Volume 7 Issue 4 (2012) e35727
  * (22) $e_{q}(x)e_{2-q}(-x)=1$ as can be seen from Eq. (8).
  * (23) https://www.britannica.com/science/Gedankenexperiment
  * (24) Stanisław Drożdż, Jarosław Kwapień, Paweł Oświȩcimka and Rafał Rak, The foreign exchange market: return distributions, multifractality, anomalous multifractality and the Epps effect, New J. Phys. 12 (2010) 105003 (23pp).
  * (25) Jarosław Kwapień and Stanisław Drożdż, Physical approach to complex systems, Phys. Rep. 515 (2012) 115-226.
  * (26) V. N. Pokrovski, ISRN Thermodynamics, ID 906136 (2013).