# Subjectivity in Unsupervised Machine Learning Model Selection

Wanyi Chen1, Mary L. Cummings 2

###### Abstract

Model selection is a necessary step in unsupervised machine learning. Despite
numerous criteria and metrics, model selection remains subjective. A high
degree of subjectivity may lead to questions about repeatability and
reproducibility of various machine learning studies and doubts about the
robustness of models deployed in the real world. Yet, the impact of modelers’
preferences on model selection outcomes remains largely unexplored. This study
uses the Hidden Markov Model as an example to investigate the subjectivity
involved in model selection. We asked 33 participants and three Large Language
Models (LLMs) to make model selections in three scenarios. Results revealed
variability and inconsistencies in both the participants’ and the LLMs’
choices, especially when different criteria and metrics disagree. Sources of
subjectivity include varying opinions on the importance of different criteria
and metrics, differing views on how parsimonious a model should be, and how
the size of a dataset should influence model selection. The results underscore
the importance of developing a more standardized way to document subjective
choices made in model selection processes.

##  1 Introduction

In a world of abundant data, unsupervised machine learning (ML) is popular for
discovering patterns and structures in data without needing labels. Selecting
the best model is a necessary step, as even slightly different models can lead
to different interpretations and decisions. For instance, psychologists can
use unsupervised ML to uncover patterns in human learning behaviors (Visser,
Raijmakers, and Molenaar 2002). Different interpretations of such models would
likely lead to different training designs, with some less optimal than others.

There are a few objectives in selecting the best model in unsupervised ML. On
the one hand, a model should accurately describe the data. On the other, a
parsimonious model is often more desirable because it is more interpretable
and less likely to overfit. However, trade-offs exist between accuracy and
parsimony. A model with more variables is more likely accurate but less
parsimonious (Dziak et al. 2020).

Several criteria guide model selection, including information criteria (IC)
such as Akaike Information Criterion (AIC) (Akaike 1974) and Bayesian
Information Criterion (BIC) (Schwarz 1978). Cross-validation metrics such as
the most consistent, best worst case, and best average may provide additional
guidance.

Despite numerous seemingly objective criteria and metrics, model selection
remains subjective. Different ICs place varying emphasis on parsimony versus
predictive accuracy (Dziak et al. 2020). ML researchers and practitioners hold
various beliefs about the importance of certain criteria or metrics, so the
definition of the “best” model becomes subjective. While it is well known that
biases in datasets cause biases in machine learning models (Buolamwini and
Gebru 2018), little is known about how the modelers’ preferences influence
model selection results. The subjective decisions made in the model selection
process can be considered the modeler’s “degree of freedom.”

The modeler’s degree of freedom yields insights into the repeatability and
reproducibility of results, essential in scientific research (Plesser 2018).
If results are not repeatable or reproducible, the research community cannot
critically assess the correctness of scientific claims made and conclusions
drawn from the results (Plesser 2018). Consequently, we cannot have confidence
in the robustness of various real-world ML models based on these research.

This study investigates the subjectivity involved in model selection using the
Hidden Markov Model (HMM) as an example of unsupervised ML. In the following
sections, we first review the background on HMM, information criterion, model
selection, and researcher degrees of freedom. Then, we describe an experiment
where 33 participants and three Large Language Models (LLMs) were asked to
select models in three scenarios.

##  2 Background and related work

An HMM is a stochastic model with two layers: a lower layer of observations
and a higher layer of hidden states that is not directly observable.
Mathematically, an HMM consists of a transition matrix that represents the
transition probabilities among $n$ different hidden states, an emission matrix
that represents the emission probabilities from $n$ hidden states to $m$
different observations, and a vector that represents the initial probability
distributions among the $n$ hidden states. An HMM’s number of free parameters
$n_{p}$ is calculated by counting the freely-estimated entries in the
transition and emission matrices. HMMs with more hidden states are more
complex and have larger $n_{p}$ values.

The Baum-Welch algorithm is typically used to train HMMs (Rabiner 1989). The
number of hidden states is unknown and must be chosen before training. Given
the same dataset, the algorithm can train different models with different
numbers of hidden states.

###  2.1 Information criterion

IC | Penalty term | Emphasis
---|---|---
AIC | $2n_{p}$ | Predictive accuracy
AIC3 | $3n_{p}$ | Predictive accuracy
AICc | $(2Tn_{p})/(T-n_{p}-1)$ | Predictive accuracy
AICu | $(2Tn_{p})/(T-n_{p}-1)+T\log(T/(T-n_{p}-1))$ | Predictive accuracy
CAIC | $n_{p}(\log T+1)$ | Parsimony
BIC | $n_{p}\log T$ | Parsimony
ABIC | $(n_{p}\log(T+2))/24$ | Depends on $T$
Table 1: Summary of ICs

Information criteria, calculated for each candidate model, are often used in
model selection. In theory, the candidate model with the minimum IC value
should be selected. Different ICs have different theoretical and practical
motivations, but they share a common form:

$IC=-2(\log L)+\texttt{penalty term}$

where $\log L$ stands for the log likelihood of the dataset given the model.
The first term $-2(\log L)$ characterizes the model’s accuracy: lower values
indicate that the model has a higher likelihood. In HMM’s case, a model’s
likelihood can be calculated by the forward algorithm (Rabiner 1989).

The model is likely to become more accurate with more parameters, so to keep
the model parsimonious, the ICs add a penalty term based on $n_{p}$. The value
of the penalty term increases as $n_{p}$ increases. Some ICs’ penalty terms
also involve $T$, the dataset size. ICs with smaller penalty terms emphasize
predictive accuracy, and ICs with larger penalty terms emphasize parsimony.
Table 1 summarizes seven common ICs’ penalty terms and emphasis. Appendix A
includes a more detailed review of the theoretical and practical motivation of
these ICs.

Among these seven ICs, AIC, AIC3, AICc, and AICu have lighter penalty weights,
making them more likely to choose more complex models and risk overfitting.
CAIC and BIC have heavier penalty weights, making them more likely to choose
simpler models and risk underfitting. ABIC’s penalty weight is lighter than
that of BIC’s and may be lighter or heavier than that of AIC’s depending on
$T$ (Dziak et al. 2020). ABIC’s penalty weight becomes heavier with larger
$T$. Because of the different penalty terms, these information criteria often
suggest different best models when given the same set of candidate models.

###  2.2 Model selection

Recognizing that different ICs often lead to different models, some
researchers have compared the performance of different ICs, with mixed
results. Costa et al. (2010) performed a Monte Carlo experiment with HMMs and
found that the dataset size, the conditional state-dependent probabilities,
and the latent transition matrix are the main factors influencing the IC test
results. Bacci et al. (2014) conducted another Monte Carlo simulation study
with HMMs. They also found that the capability of the different ICs in
detecting the true number of hidden states varies case by case. Particularly,
lower values of persistence probabilities in the same state and/or a greater
uncertainty in the allocation of the observations to the latent states may
lead to worse performance of some ICs (Bacci, Pandolfi, and Pennoni 2014).

These Monte Carlo experiments assume some ground truth exists. However, in
many real-world applications, it is impossible to know the true latent
transition matrix or many other salient factors. Therefore, it still remains
uncertain whether choosing some ICs over others will yield better model
selection results for a particular real-world scenario.

Other researchers propose abandoning the use of ICs altogether. Celeux et al.
(2008) proposed to use cross-validation to instead assess the number of hidden
states in HMMs, i.e., partitioning the dataset into different folds and
iteratively assigning different folds as training and testing sets. The model
that maximizes the cross-validated likelihood on testing sets is selected.
This method selects models with the best predictive performances (Celeux and
Durand 2008). However, this method is often computationally intensive (Pohle
et al. 2017), and evidence shows that it tends not to outperform BIC (Costa
and De Angelis 2010).

Some researchers have developed pragmatic guides for practitioners engaged in
model selection. Dziak et al. (2020) advised practitioners to choose AIC if
underfitting is a concern and choose the more parsimonious BIC if overfitting
is more of a concern. Pohle et al. (2017) advised that in addition to using
the ICs, practitioners should also use their expert knowledge to narrow the
range of candidate models and use model validation methods to ensure models
accurately represent the data. These pragmatic guides attempt to provide
selection structure, but little is known about human preferences in this
process.

Beyond HMMs, model section occurs in many other fields and contexts. For
example, statisticians are interested in analyzing data clusters by means of
mixture distributions, and the ICs can help them select the best mixture
models (Bozdogan 1994). The ICs are also used for determining the number of
significant components in Principle Component Analyses, which has a wide range
of applications (Bai, Choi, and Fujikoshi 2018). For the purposes of this
study, we use HMM as an example of unsupervised ML to investigate how human
preferences influence model selection.

###  2.3 Researcher degrees of freedom

Researchers’ subjective choices in data collection, processing, and analysis
cause studies to have different results and sometimes support opposite
hypotheses (Simmons, Nelson, and Simonsohn 2011; Huntington-Klein et al.
2021). Simmons et al. (2011) first coined the term “researcher degrees of
freedom” when they discovered that a seemingly small change in researchers’
subjective choices could significantly increase the false-positive rate in
psychology studies. Similarly, Huntington-Klein et al. (2021) found that
unreported researcher variation could decrease the validity of study results.

ML researchers recognize that biases in datasets cause biases in resulting
models (Buolamwini and Gebru 2018), yet practitioner-induced biases are less
discussed. More recently, Cummings et al. (2021) demonstrated that subjective
judgments made in the development and interpretation of ML models could lead
to inconsistent results and errors of commission and omission. Still, none
have explored the impact of researcher degrees of freedom in unsupervised
learning model selection, which is the aim of this study.

##  3 Method

We created three model selection scenarios by training HMMs using publicly
available replay files from Star-Craft II (Blizzard 2021), a popular online
multiplayer Real-Time Strategy game. We focused on Zerg (a role in the game)
and grouped its more than 100 different actions into 20 categories. The
sequences of player actions, called replays, are our observation sequences.

We created three datasets corresponding to three model selection scenarios: a
small dataset (scenario 1) comprising 30 replay files, a mid dataset (scenario
2) comprising 490 replay files, and a large dataset (scenario 3) comprising
2980 replay files. All files were randomly sampled for all datasets. For each
dataset, we trained 19 candidate HMMs with the number of hidden states ranging
from 2 to 20. Then, we calculated the IC values for each candidate model
(Appendix B). Table 2 summarizes the best model chosen by each IC in each
scenario. The models are denoted by their number of hidden states. For
example, model 9 refers to the candidate model with 9 hidden states.

IC | Scenario 1 | Scenario 2 | Scenario 3
---|---|---|---
AIC | model 9 | model 19 | model 20
AIC3 | _model 6_ | model 17 | model 20
AICc | model 9 | model 17 | model 20
AICu | _model 6_ | model 13 | model 20
CAIC | model 2 | model 4 | _model 13_
BIC | model 2 | _model 7_ | model 14
ABIC | model 4 | _model 7_ | model 14
Table 2: Best model selected by each IC in each scenario. Bold = best average
model, italics = most consistent model, underline = best worst case model.

Next, for the candidate models chosen by at least one IC, we ran additional
cross-validation tests. For example, in scenario 1, candidate models 2, 4, 6,
and 9 appeared at least once in the second column of Table 2, so we further
investigated these models. We evaluated their negative loglikelihood terms
$-2(\log L)$ using 5-fold cross-validation. For scenarios 2 and 3, the
datasets are larger, so we ran 10-fold cross-validation instead of 5-fold. We
report the full results in Appendix C.

Table 2 also reflects the “best average” metric selects the model with the
lowest $-2(\log L)$ term average across the 5 or 10 folds. Lower $-2(\log L)$
term values suggest that the model has higher likelihood. We also calculated
the $-2(\log L)$ term’s standard deviation (SD) across the folds. The “most
consistent” metric selects the model with the smallest SD. Additionally, for
each model, the largest $-2(\log L)$ term value among the folds is considered
the “worst case.” The “best worst case” metric selects the model with the
smallest worst-case value.

Given these models, we designed a self-paced survey to measure modelers’
preferences in the three model selection scenarios (Appendix D). First,
participants completed an informed consent form and a demographics survey
(Appendix E). The demographics survey also included a seven-question risk
propensity scale survey that measures an individual’s general risk-taking
tendencies (Meertens and Lion 2008). We hypothesized that one’s risk
propensity might correlate to model selection choices. Then, participants were
given background information about HMMs and the ICs.

Next, participants completed an untimed model selection survey, which gave the
participants the three model selection scenarios in randomized order.
Participants did not know that the models were trained on StarCraft II data.
For each scenario, the survey presents the best models chosen by various ICs,
the cross-validation results, and the summaries. The participants also have
optional access to the tables in Appendix B by clicking a “click to expand”
option. The participants are then asked to select only one “best” model based
on available information, indicate which factors influenced their decisions
the most, and briefly explain their reasoning (see Appendix F for the full
survey).

Most participants completed the survey in 20 to 40 minutes. The participants
were recruited as volunteers but were each given a $20 Amazon gift card as an
appreciation. This study was approved by the Institutional Review Board.

In addition to the human participants, we also incorporated three LLMs into
our study: ChatGPT 3.5, ChatGPT 4.0, and Google Bard. Given that LLMs are
trained on a vast corpus of human texts, encompassing numerous research papers
on model selection, we hypothesized that the LLMs might emulate the most
typical human responses to model selection tasks. Given that model selection
is an onerous task, we also wanted to explore whether the LLMs could be
helpful tools for assisting this task. For each model selection scenario, the
LLMs were presented with identical information as our human participants,
including the dataset size, the best models chosen by various ICs, and a
summary of the cross-validation results (Appendix I). The next section
presents the results.

##  4 Results and discussion

We recruited a total of 33 participants, all with some ML experience. Ages
ranged from 20-45 yrs, with an average of 27 yrs (SD = 5 yrs). The group
included 8 women and 25 men. 27 were undergraduate or graduate students
majoring in computer science or related fields, and 6 were professionals in
technical roles. Appendix G provides more details on the demographics.

The participants’ risk propensity scores range from 2.14 to 6.17, with a mean
of 3.77 (SD = 1.04). Higher scores indicate more risk-taking. Typical risk
propensity scale scores range from 2.00 to 7.00, with a mean of 4.63 (SD =
1.23) (Meertens and Lion 2008). Our mean scores indicate that our group of
participants is considerably more risk-averse than the original group,
although we cannot establish statistical significance without access to the
original data.

###  4.1 Participant model selection preferences

In two out of three scenarios, participants’ model selection choices varied.
Figure 1 presents participants’ choices in each scenario. When two of the
three cross-validation metrics (best average, most consistent, best worst
case) point to the same model, most participants chose that model. For
example, in scenario 2, model 7 is both the most consistent model and the
model with the best worst case. 31 out of 33 participants chose model 7 as the
best model. A few participants based their choices on evidence that more than
one cross-validation metric agreed on model 7. “7 seemed like the obvious
choice because it was both the most consistent and the best worst case,” wrote
one participant. Another participant also noted, “My choice was based
primarily on consistency results and best worst case.”

Figure 1: Model selection results

However, in the other two scenarios, all three cross-validation metrics point
to different models, resulting in significantly more disagreements among
participants. Some participants had almost exact opposite reasons for
selecting models. For instance, in scenario 1, the simplest HMM possible, a
model with 2 hidden states, is one of the candidate models. 13 participants
chose this model, with one explaining, “When it comes to model selection, I
often prefer parsimony over accuracy.” Similarly, another participant simply
wrote “the fewest number of states” when asked to explain his or her
reasoning.

On the contrary, some participants refrained from choosing such extreme
values. One participant wrote, “I just don’t like 2 because it is the smallest
possible number. Explaining everything with just two variables sounds
dangerous.” Similarly, another participant noted, “Given 19 observation
categories, only 2 latent states (as suggested by BIC) might be too
reductionist.” A third participant stated, “I usually do not select models at
the boundary of the range, e.g., I don’t choose 2 in range 2-20.” These
different attitudes towards extreme simplicity contributed to disagreements
about the best model in this scenario.

Participants also held divergent views on how dataset size should influence
model selection. The survey did not explicitly ask how the participants
reasoned about dataset sizes, and only six participants (18.2%) mentioned
dataset size as a factor they considered. One participant believed that the
“most consistent” metric was more important for larger dataset sizes, stating,
“Given the large size of the dataset, results on consistency are likely to
hold and be a good reflection of the strength of the model.” However, another
participant held the opposite view, noting, “If the data size is small, I may
prefer a more consistent model.”

There was also disagreement about the importance of accuracy for models
trained on small datasets. One person thought accurate predictions were
difficult with a small sample size, stating, “The sample size is really small
compared with the previous two datasets. I think it’s difficult to get an
accurate prediction.” Conversely, another participant prioritized accuracy
over parsimony for the same scenario, explaining, “There was not much data, so
I selected accuracy over parsimony for information criterion.”

Information criteria is another factor that affected participants’ selections.
Different participants placed importance on different ICs. Most participants
indicated that more than one IC was important. Table 3 summarizes their
responses. The second column is the number of participants who considered the
corresponding IC listed in the first column. Across all three scenarios, BIC
is considered most frequently. Other ICs have varying degrees of popularity
and are not consistent across scenarios. Among the seven ICs, BIC (Schwarz
1978) is also the most cited one, reaching over 55,000 citations to date. It
is possible that the participants considered BIC more frequently because it is
more well-known.

Additionally, ten participants (30.3%) commented that they would like to know
more context since it might influence their model selection preferences. Our
survey did not give the participants context about the dataset, how it was
processed, and the purpose of the model. One participant commented, “The
metrics above are useful in a general sense, but domain knowledge of the
problem and how and where the models will be used or deployed also play an
important role in model selection.”

IC considered | Count
---|---
BIC | 25
CAIC | 19
AIC3 | 17
AICu | 15
AIC | 13
ABIC | 11
AICc | 6
(a)

IC considered | Count
---|---
BIC | 30
ABIC | 26
CAIC | 14
AIC | 10
AIC3 | 7
AICc | 6
AICu | 6
(b)

IC considered | Count
---|---
BIC | 25
AIC | 19
CAIC | 18
ABIC | 18
AICu | 15
AICc | 13
AIC3 | 12
(c)

Table 3: IC considered for each scenario Figure 2: Individual consistency
scores

| Scenario 1 | Scenario 2 | Scenario 3
---|---|---|---
a. One or more information criteria | 4 | 4 | 5
b. Best average performance | 4 | 5 | 4
c. Most consistent performance | 4 | 4 | 4
d. Best worst-case performance | 2 | 2 | 2
(a)

Pairwise comparison | Scenario 1 | Scenario 2 | Scenario 3 | Pairwise inconsistency score
---|---|---|---|---
$a$-$b$ | $a=b$ | $a<b$ | $a>b$ | 2
$a$-$c$ | $a=c$ | $a=c$ | $a>c$ | 1
$a$-$d$ | $a>d$ | $a>d$ | $a>d$ | 0
$b$-$c$ | $b=c$ | $b>c$ | $b=c$ | 1
$b$-$d$ | $b>d$ | $b>d$ | $b>d$ | 0
$c$-$d$ | $c>d$ | $c>d$ | $c>d$ | 0
(b)

Table 4: Example consistency score calculation

###  4.2 Participant inconsistencies

We analyzed participants’ choice consistency across scenarios in terms of
whether their criteria preferences changed between scenarios. The survey asked
each participant to rate the importance of four factors (first column of Table
4(a)) in each scenario. Each factor was rated on a Likert scale of 1 to 5,
with 1 being the least important and 5 being the most important. Table 4(a)
shows one participant’s responses as an example.

Then, we compared participants’ ratings across scenarios. For instance, did
participants always think that “$a.$ one or more information criteria” is more
important than “$b.$ best average performance in cross-validation”? If they
rated choice $a$ the same as choice $b$ in one scenario, we denote that as
$a=b$. Similarly, $a>b$ means that a participant rated choice $a$ higher than
choice $b$, and $a<b$ means that the participant rated choice $b$ higher than
choice $a$. If the comparison symbol is consistent across all three scenarios,
we denote a pairwise inconsistency score of 0 (the $a$-$d$ row in Table 4(b)
offers an example). Similarly, a pairwise inconsistency score of 1 means the
comparison symbol in one scenario differs from the other two (the $a$-$c$ row
offers an example). A pairwise inconsistency score of 2 means the comparison
symbol differs across all three scenarios (the $a$-$b$ row offers an example).

We added all six pairwise inconsistency scores to obtain an overall
consistency score for each participant. In the Table 4 example, the
participant’s consistency score is 4. Higher scores indicate more
inconsistency. The overall mean score is 3.70 (SD = 2.53). Figure 2 shows the
score distribution. We used k-means clustering to define three clusters. We
labeled the six participants (18.2%) whose scores ranged from 0 to 1 as
“consistent”, the fourteen participants (42.4%) whose scores ranged from 2 to
4 as “slightly inconsistent,” and the thirteen participants (39.4%) whose
scores ranged from 5 to 8 as “inconsistent.”

By this definition, nearly 40% of participants were inconsistent in their
model selection preferences. The cause of the inconsistencies is not clear.
Only two of thirteen inconsistent participants explained that they favored
different ICs and cross-validation metrics because of different dataset sizes.
The others did not explicitly explain why their preferences varied across
scenarios, which is inline with behavioral economics research, which has shown
that individuals can exhibit inconsistent preferences when presented with the
same choice in different forms (Kahneman and Tversky 2013). While no two
scenarios are identical in our survey, participants faced similar choices
across scenarios, relying on information presented to them in similar forms.
The high degree of inconsistency observed in our study suggests that
predicting an individual’s model selection preferences in a new scenario is
challenging, unless many criteria and metrics agree and the choice relatively
straightforward.

###  4.3 Participant correlations

We investigated whether consistency scores correlated with gender, risk
attitude, and years of ML experience. For risk attitude, we labeled the six
participants (18.2%) whose risk propensity scale scores were one standard
deviation below the mean as “risk avoiders”, the five participants (15.2%)
whose scores were one standard deviation above the mean as “risk takers,” and
the remaining as “neutral.” For years of ML experience, the survey asked
whether participants had less than one year, one to three years, or more than
three years of experience. Since only four participants indicated “less than
one year”, we merged the first two experience levels. We clustered the 23
participants (69.7%) with $\leq 3$ years of experience into one group and the
10 participants (30.3%) with $>3$ years of experience into another group.

Then, the Fisher’s exact test determined whether nonrandom associations exist.
Fisher’s exact test is an alternative to the $\chi^{2}$-test and is
recommended when the sample size is small (Freeman and Halton 1951). We found
no correlations between consistency scores and gender (p = 0.88), risk
attitudes (p = 0.78), and years of ML experience (p = 0.34) at the
significance level of 0.05.

We also ran Fisher’s exact test to investigate whether one’s gender, risk
attitude, years of ML experience, and individual consistency correlate with
one’s preferences in model selection. Scenario 2 is omitted for this analysis
because 31 out of 33 participants chose the same model. We report the results
and p-values in Appendix H. We found no correlations at the significance level
of 0.05.

###  4.4 LLM model selection preferences

LLM |  $1^{st}$ try |  $2^{nd}$ try |  $3^{rd}$ try
---|---|---|---
ChatGPT 3.5 | model 2 | model 6 | model 6
ChatGPT 4.0 | model 2 | model 6 | model 4
Google Bard | model 2 | model 2 | model 6
(a)

LLM |  $1^{st}$ try |  $2^{nd}$ try |  $3^{rd}$ try
---|---|---|---
ChatGPT 3.5 | model 7 | model 4 | model 7
ChatGPT 4.0 | model 7 | model 7 | model 7
Google Bard | model 7 | model 7 | model 7
(b)

LLM |  $1^{st}$ try |  $2^{nd}$ try |  $3^{rd}$ try
---|---|---|---
ChatGPT 3.5 | model 14 | model 20 | model 20
ChatGPT 4.0 | model 14 | model 14 | model 14
Google Bard | model 14 | model 14 | model 14
(c)

Table 5: LLM model selection results

Unlike the participants who only took the survey once, LLMs have the capacity
to answer the same question with multiple responses. This can be achieved by
either selecting “regenerate response” or by starting a new conversation with
the identical prompt. For each scenario, we repeated the prompt three times.
Table 5 summarizes the “best” models selected by each LLM in each scenario.

In two out of three scenarios, the LLMs’ response patterns closely resembled
those of the human participants. In scenario 1, the LLMs’ responses varied the
most. Responses were evenly split between models 2 and 6, with one mention of
model 4. Human preferences paralleled this, with roughly equal preferences for
models 2 and 6 and a lesser preference for model 4. In scenario 2, both LLMs
and humans had the highest amount of agreements. The vast majority agreed that
model 7 is the best model, as this is a relatively easy case. However, in
scenario 3, while LLMs gravitated towards model 14, human responses had more
variance. The reason for this difference is unknown.

When asked to state model selection reasons, different LLMs also showed
varying preferences for parsimony. For instance, in scenario 3, ChatGPT 3.5
stated, “The preference for simpler models (14 hidden states) make it a
reasonable choice. The preference for simpler models is rooted in the
principle of Occam’s razor, which states that among competing hypotheses or
models, the simplest one should be preferred until evidence suggests
otherwise.” In comparison, instead of simply favoring simpler models, ChatGPT
4.0 stated, “It’s also important to consider the trade-off between model
simplicity and the model’s ability to capture complex structures in the data.”
These responses mirrored the humans’ responses, where different participants
showed different attitudes towards parsimony.

Moreover, similar to 30.3% of the participants, all LLMs recognized that model
selection may be context-dependent. For instance, ChatGPT 4.0 commented, “This
[model selection] also largely depends on the specific context and what you
prioritize in model performance. You may consider factors like
interpretability, computation cost, and the consequences of misclassification
when making a final decision.” Similarly, Google Bard said, “I would also
consider the following factors when selecting the best model: the domain
knowledge of the problem, the computational resources available, and the time
constraints.”

Like their human counterparts, all LLMs showed inconsistencies in their
answers. Each provided different answers to the same prompt in at least one
scenario. We prompted them to rate the importance of the four factors shown in
Table 4(a), in the same way we asked the human participants to rate these
factors (see Appendix I for an example prompt). Using the method detailed in
section 4.2, we derived a consistency score for each LLM. Both ChatGPT 3.5 and
Google Bard scored 3 (slightly inconsistent), while ChatGPT 4.0 scored 7
(inconsistent).

In summary, the LLMs mirrored human behavior in model selection. Both the LLMs
and the humans expressed varying preferences for parsimony, recognized that
model selection may be context-dependent, and showed notable inconsistencies
in their decision processes.

##  5 Conclusion

The results highlight that subjectivity can greatly influence the seemingly
objective model selection process, especially when different criteria and
metrics disagree. Furthermore, a large percentage of human participants and
all the LLMs are inconsistent in their model selection preferences, with no
trends seen in experience, background, or LLM versions.

Researchers and practitioners must recognize their inherent subjectivity and
potential inconsistencies. High levels of subjectivity and inconsistencies can
cast doubts on the assertions made in ML research and compromise the
robustness of ML models. When leveraging tools like LLMs to aid model
selection, it’s equally important to be cognizant of the limitations and
inconsistencies inherent in these tools.

To improve the repeatability and reproducibility of ML studies, the ML
community should standardize the reporting of subjective choices in model
selection. Our study revealed several sources of subjectivity. First, both the
participants and the LLMs had very different views on how parsimonious models
should be. The preferred degree of parsimony should be reported and justified
in application contexts. This preference should also guide which ICs are
considered more important because different ICs quantify different trade-offs
between parsimony and accuracy.

Another source of subjectivity is how different dataset sizes should influence
model selection. Multiple studies agree that AIC (and other ICs that emphasize
accuracy) is preferred for small datasets because underfitting is more likely
a concern, while BIC (and other ICs that emphasize parsimony) is preferred for
large datasets because overfitting is more likely a concern (Costa and De
Angelis 2010; Dziak et al. 2020). Any model selection preferences that violate
this recommendation should be reported and justified.

If researchers and practitioners run cross-validation tests in addition to
information criteria, then cross-validation metrics should also be reported,
since it is a major source of subjectivity. Furthermore, if many metrics
disagree, researchers and practitioners should consider whether the model
should be trained in the first place. For instance, if many different ICs and
cross-validation metrics suggest HMMs with widely-varying hidden states,
researchers should consider whether an HMM best describes the data.

###  5.1 Limitations and future work

One limitation of our study is that it is unknown whether multiple ICs and
cross-validation metrics are used in practice. Some studies use BIC alone to
select models without discussing alternative ICs (Le Strat and Carrat 1999;
Wang et al. 2016).

Another limitation is that our survey did not provide any context about the
dataset. As ten participants and all three LLMs commented, context might
influence their model selection preferences. Future studies should examine how
different contexts affect model selection. However, providing more context
might introduce more confounding variables. For example, had participants been
informed that the models were trained on StarCraft II data, their prior
knowledge of StarCraft II might have influenced their preferences. Hence,
future studies will need careful experimental designs to effectively
disentangle these confounding influences.

This study revealed that many participants were either unfamiliar with model
selection guidelines in the literature or chose to ignore them, as evidenced
by varying opinions regarding dataset sizes. Future studies should examine
whether or when researchers and practitioners are inclined to follow given
technical guidelines and why. Another direction for future studies is whether
a data visualization tool could make the model selection process more
transparent. Five participants (15.2%) mentioned the information presented in
the survey seemed overwhelming, and some visualizations of the data might
help.

Ultimately, a better understanding of the researcher degrees of freedom in
model selection will inform a more standardized model selection process, which
will help improve the repeatability and reproducibility of various ML studies
and ultimately lead to more confidence in various ML models deployed in the
real world.

## Acknowledgements

This research was sponsored, in part, by Charles River Analytics and the
Department of the Navy.

## References

  * Akaike (1974) Akaike, H. 1974.  A new look at the statistical model identification.  _IEEE transactions on automatic control_ , 19(6): 716–723.
  * Bacci, Pandolfi, and Pennoni (2014) Bacci, S.; Pandolfi, S.; and Pennoni, F. 2014.  A comparison of some criteria for states selection in the latent Markov model for longitudinal data.  _Advances in Data Analysis and Classification_ , 8: 125–145.
  * Bai, Choi, and Fujikoshi (2018) Bai, Z.; Choi, K. P.; and Fujikoshi, Y. 2018.  Consistency of AIC and BIC in estimating the number of significant components in high-dimensional principal component analysis.  _The Annals of Statistics_ , 46(3): 1050–1076.
  * Blizzard (2021) Blizzard. 2021.  s2client-proto.  https://github.com/Blizzard/s2client-proto.  Accessed: 2021-09-01.
  * Bozdogan (1994) Bozdogan, H. 1994.  Mixture-model cluster analysis using model selection criteria and a new informational measure of complexity.  In _Proceedings of the first US/Japan conference on the frontiers of statistical modeling: An informational approach_ , 69–113. Springer.
  * Buolamwini and Gebru (2018) Buolamwini, J.; and Gebru, T. 2018.  Gender shades: Intersectional accuracy disparities in commercial gender classification.  In _Conference on fairness, accountability and transparency_ , 77–91. PMLR.
  * Celeux and Durand (2008) Celeux, G.; and Durand, J.-B. 2008.  Selecting hidden Markov model state number with cross-validated likelihood.  _Computational Statistics_ , 23(4): 541–564.
  * Costa and De Angelis (2010) Costa, M.; and De Angelis, L. 2010.  Model selection in hidden Markov models: a simulation study.  Technical report, Department of Statistics, University of Bologna.
  * Cummings and Li (2021) Cummings, M. L.; and Li, S. 2021.  Subjectivity in the creation of machine learning models.  _ACM Journal of Data and Information Quality_ , 13(2): 1–19.
  * Dziak et al. (2020) Dziak, J. J.; Coffman, D. L.; Lanza, S. T.; Li, R.; and Jermiin, L. S. 2020.  Sensitivity and specificity of information criteria.  _Briefings in bioinformatics_ , 21(2): 553–565.
  * Freeman and Halton (1951) Freeman, G.; and Halton, J. H. 1951.  Note on an exact treatment of contingency, goodness of fit and other problems of significance.  _Biometrika_ , 38(1/2): 141–149.
  * Huntington-Klein et al. (2021) Huntington-Klein, N.; Arenas, A.; Beam, E.; Bertoni, M.; Bloem, J. R.; Burli, P.; Chen, N.; Grieco, P.; Ekpe, G.; Pugatch, T.; et al. 2021.  The influence of hidden researcher decisions in applied microeconomics.  _Economic Inquiry_ , 59(3): 944–960.
  * Kahneman and Tversky (2013) Kahneman, D.; and Tversky, A. 2013.  Prospect theory: An analysis of decision under risk.  In _Handbook of the fundamentals of financial decision making: Part I_ , 99–127. World Scientific.
  * Le Strat and Carrat (1999) Le Strat, Y.; and Carrat, F. 1999.  Monitoring epidemiologic surveillance data using hidden Markov models.  _Statistics in medicine_ , 18(24): 3463–3478.
  * Meertens and Lion (2008) Meertens, R. M.; and Lion, R. 2008.  Measuring an individual’s tendency to take risks: the risk propensity scale 1.  _Journal of applied social psychology_ , 38(6): 1506–1520.
  * Plesser (2018) Plesser, H. E. 2018.  Reproducibility vs. replicability: a brief history of a confused terminology.  _Frontiers in neuroinformatics_ , 11: 76.
  * Pohle et al. (2017) Pohle, J.; Langrock, R.; van Beest, F. M.; and Schmidt, N. M. 2017.  Selecting the number of states in hidden Markov models: pragmatic solutions illustrated using animal movement.  _Journal of Agricultural, Biological and Environmental Statistics_ , 22(3): 270–293.
  * Rabiner (1989) Rabiner, L. R. 1989.  A tutorial on hidden Markov models and selected applications in speech recognition.  _Proceedings of the IEEE_ , 77(2): 257–286.
  * Schwarz (1978) Schwarz, G. 1978.  Estimating the dimension of a model.  _The annals of statistics_ , 461–464.
  * Simmons, Nelson, and Simonsohn (2011) Simmons, J. P.; Nelson, L. D.; and Simonsohn, U. 2011.  False-positive psychology: Undisclosed flexibility in data collection and analysis allows presenting anything as significant.  _Psychological science_ , 22(11): 1359–1366.
  * Visser, Raijmakers, and Molenaar (2002) Visser, I.; Raijmakers, M. E.; and Molenaar, P. 2002.  Fitting hidden Markov models to psychological data.  _Scientific Programming_ , 10(3): 185–199.
  * Wang et al. (2016) Wang, F.; Tan, S.; Yang, Y.; and Shi, H. 2016.  Hidden Markov model-based fault detection approach for a multimode process.  _Industrial & Engineering Chemistry Research_, 55(16): 4613–4621.


# Bridging Binarization: Causal Inference with Dichotomized Continuous
Treatments

Kaitlyn J. Lee
Division of Biostatistics
UC Berkeley  Alan Hubbard
Division of Biostatistics
UC Berkeley  Alejandro Schuler
Division of Biostatistics
UC Berkeley

###### Abstract

The average treatment effect (ATE) is a common parameter estimated in causal
inference literature, but it is only defined for binary treatments. Thus,
despite concerns raised by some researchers, many studies seeking to estimate
the causal effect of a continuous treatment create a new binary treatment
variable by dichotomizing the continuous values into two categories. In this
paper, we affirm binarization as a statistically valid method for answering
causal questions about continuous treatments by showing the equivalence
between the binarized ATE and the difference in the average outcomes of two
specific modified treatment policies. These policies impose cut-offs
corresponding to the binarized treatment variable and assume preservation of
relative self-selection. Relative self-selection is the ratio of the
probability density of an individual having an exposure equal to one value of
the continuous treatment variable versus another. The policies assume that,
for any two values of the treatment variable with non-zero probability density
after the cut-off, this ratio will remain unchanged. Through this equivalence,
we clarify the assumptions underlying binarization and discuss how to properly
interpret the resulting estimator. Additionally, we introduce a new target
parameter that can be computed after binarization that considers the status-
quo world. We argue that this parameter addresses more relevant causal
questions than the traditional binarized ATE parameter. Finally, we present a
simulation study to illustrate the implications of these assumptions when
analyzing data and to demonstrate how to correctly implement estimators of the
parameters discussed.

##  1 Introduction

Causal inference methodology historically has considered settings with binary
treatment variables. For example, Neyman focused on estimating the average
treatment effect (ATE), which compares the average outcome under treatment to
the average outcome under no treatment [1]. The ATE continues to be a popular
target parameter in causal literature focused on binary treatments. However,
across a variety of disciplines, there is growing interest in quantifying the
effects of continuous exposures. For example, recent research has explored the
causal effects of income inequality [2], alcohol consumption [3], and
pesticide use [4].

The ATE is not defined for a continuous exposure; however, one popular method
of quantifying the causal effect of a continuous exposure is by binarizing or
dichotimizing the treatment variable, which enables researchers to target a
parameter somewhat akin to the ATE. This method involves creating a new dummy
variable that equals 1 when the exposure of interest is within some range of
values and 0 otherwise. Then, researchers calculate the ATE using this new
dummy variable as the exposure. We term this estimand the“binarized ATE.”

What does this look like in practice? A 2015 paper by Hu et al. investigated
the relationship between pesticide exposure and health effects in Chinese
farmers [4]. Individuals were divided into two groups: those with relatively
high pesticide exposure levels (more than 50 pesticide applications in
2009-2011) and those with relatively low pesticide exposure levels (less than
50 pesticide applications in the same three year period). The researchers
implemented a regression estimator including baseline covariates. They found
evidence suggestive of some causal relationship between “high” vs. “low”
levels of pesticide use and abnormal health test results (including blood
tests, nerve conduction studies, and neurological exams).

There are numerous examples of researchers binarizing continuous treatments in
public health contexts. For example, Hu et al. investigated the relationship
between in a “low-risk” group and risk of type 2 diabetes in women [5]. In
another example, Wu et al. calculated the effect of air pollution and
mortality [6]. In addition, researchers employ the binarized ATE across the
social sciences, like in political science and economics. Kurtz and Lauretig
investigated the effect of economic liberatlization on political protest
acitivty [7]. Carneiro, Lokshin, and Umapathi looked at the effect of
increased schooling on income in Indonesia [8].

In this paper, we show that the binarized ATE is the same as the difference in
average outcome under two related modified treatment policies (MTP). MTPs are
an example of a population-level intervention in which one considers a
counterfactual world in which changes to underlying distributions of different
variables are considered [9]. In particular, MTPs are policies that change the
underlying distribution of the treatment variable in a population such that
the new distribution depends on covariates as well as an individual’s natural
value of treatment, as defined by Hanuese and Rotnitzky [10]. An individual’s
natural value of treatment is the value of treatment an individual would
receive in the status-quo world, prior to any intervention, which could differ
across individuals. By considering shifts in the underlying population
distribution of the treatment, MTPs allow for stochastic interventions rather
than deterministic interventions. These types of interventions can therefore
represent counterfactual worlds that are more realistic in modeling the type
of behaviors that researchers may expect in a target population.

Utilizing this framework, we show that the binarized ATE estimates the
difference in expected outcome under two policies that impose cut-offs
corresponding to the cut-offs used to generate the binarized treatment
variable and assume preservation of relative self-selection. Relative self-
selection is the ratio of the probability density of an individual having an
exposure equal to one value of the continuous treatment variable versus
another. The policies assume that, for any two values of the treatment
variable with non- zero probability density after the cut-off, this ratio will
remain unchanged.

This equivalence addresses many of the common critiques put forward by methods
researchers have regarding binarization in causal inference. Even textbooks
written on epidemiological methods heavily dissuade researchers from
categorizing/binarizing treatment variables [11]. Some researchers contend
that binarization violates the Stable Unit Treatment Values Assumption
(SUTVA). SUTVA is one of the underlying assumptions that researchers make when
prescribing a causal meaning to their statistical estimate of the ATE [12].
VanderWeele and co-authors argue that binarization violates this assumption as
encoding a continuous treatment into a binary one implies that there are
several values of treatment embedded in one value of the binary treatment
variable [12] [13]. Rothman, Greenland, and Lash make a similar argument,
stating that such “open-ended categories” such as cut-offs are harmful because
one category can contain many different exposure values or confounder effects
[11]. However, our proposed MTP circumvents these issues by defining the
precise mechanism by which the value of the treatment changes. Specifically,
we can re-write the average outcome under two seemingly deterministic
binarized treatment levels in terms of the average outcomes under two policies
that change the underlying distribution of treatment in a fully specified way.
In addition to this critique, Bennette and Vickers contend that binarization
assumes“homogeneity of risk within categories” [14]. However, our application
of MTPs in this paper does not imply that the effect is the same for everyone;
rather, in thinking about the effect of a modified treatment policy at a
population level, we are considering a target parameter that is an average of
the treatment effect across the distribution of the treatment of interest in
the target population.

This type of dichotmization was discussed by Stitleman et al. [15]. These
interventions were also mentioned by Díaz and van der Laan as one example of
an intervention which imposes some type of cut-off on the distribution of
treatment [16]. McCoy et al. introduced the “average regional effect,” which
considers how the average outcome changes under modified treatment policies
which preserve relative self-selection [17]. Rose and Shem-Tov consider
binarization in the context of noncompliance [18]. Additionally, the binarized
ATE was considered in the context of non-random missingness by van der Laan et
al. [19]. The estimand in [19] reduces to the BATE when assuming missingness
at random. In this paper, we aim to provide a thorough and detailed treatment
of the assumptions underlying the BATE, as well as examples of when the
assumptions may be more or less plausible, to provide practioners with more
clarity regarding whether the BATE is the answer to their particular research
question. We also present a novel estimand, the policy effect of binarization,
which extends the additive risk estimand considered in [15]. We believe that
this estimand answers a more relevant causal question than the BATE. In
addition, we provide details for regression, inverse probability weighted, and
augmented inverse probability weighted estimators that, under certain
conditions outlined, will consistently estimate the BATE and the policy effect
of binarization. To our knowledge, these estimators and their underlying
assumptions have not been discussed in the context of binarization.

The rest of the paper proceeds as follows. In Section 2, we present our proof
of equivalence between the binarized ATE and the effect of modified treatment
polices that restrict the value of the treatment variable will preserving
relative self-selection. We also introduce the policy effect of binarization
and discuss its relevance. Finally, we introduce regression adjustment-based,
inverse probability weighting-based, and doubly-robust estimators for the two
target parameters. Section 3 follows with a simulation study in which we
implement both estimators for the binarized ATE and our new target parameter.
Lastly, in Section 4, we provide a discussion of our results and potential
directions for future work.

##  2 Methods

###  2.1 Setup and Notation

Consider an observational study with baseline covariates $W\in\mathbb{R}^{p}$,
treatment (possibly continuous) $A\in\mathbb{R}$, and a single outcome
$Y\in\mathbb{R}$. Let the observed data be $O=(W,A,Y)$. Let supp$(A)$ and
supp$(W)$ denote the support of treatments and covariates, respectively.

We denote the data generating distribution as $P_{0}$ for the observed data
$O$. Assume that $O_{1},...,O_{n}$ are $n$ i.i.d. draws of $O\sim P_{0}$. If
$p_{Y,A,W}(y,a,w)$ is the joint distribution of the data, then we can write
the decomposition $p_{Y,A,W}(y,a,w)=p_{Y|A,W}(y,a,w)p_{A|W}(a,w)p_{W}(w)$. We
can define our statistical model $\mathcal{M}$ as the set of all possible
distributions. We make no assumptions about these densities, and thus our
statistical model is nonparametric.

Let our statistical target parameter $\Psi:\mathcal{M}\rightarrow\mathbb{R}$
be a mapping from our statistical model to the parameter space (in this case,
real line). Then, our target parameter can be thought of as a function of
different distributions where $\Psi(P_{0})$ gives us our true estimand of
interest.

Each set of observed values $\\{O_{1},...,O_{n}\\}$ can be fully described
using an empirical probability distribution $P_{n}$ that puts point mass
$\frac{1}{n}$ at each observation $O_{i}$. Each draw from $P_{0}$ generates a
new $P_{n}$; therefore, $P_{n}$ is also a random distribution. We can now
define our statistical estimator $\hat{\Psi}(P_{n})$ as a function of our
observed data that maps to the parameter space $\mathbb{R}$. Since $P_{n}$ is
random, $\hat{\Psi}(P_{n})$ is also random, meaning that each set of observed
data $\\{O_{1},...,O_{n}\\}$ can correspond to a different value of the
estimator.

Let $\mathcal{A}$ be a subset of supp$(A)$ and
$\mathcal{A}^{\mathsf{c}}=\text{supp}(A)/\mathcal{A}$ be the complement of
$\mathcal{A}$. Consider a binary treatment variable $T$. An individual with a
natural value of treatment $A=a$ is assigned a value of $T=1$ when
$a\in\mathcal{A}$ and assigned $T=0$ when $a\in\mathcal{A}^{\mathsf{c}}$. In
the traditional binarization setting, $\mathcal{A}$ is some interval on
$\mathbb{R}$ (for example, $T=1$ when $A\geq a$, and $T=0$ when $A<a$, for
some $a\in\mathbb{R}$).

As is standard practice in causal inference literature, we will adopt the
potential outcome framework [20]. The framework assumes that, for each
individual and each permissible value of treatment $A=a$, there exists a
counterfactual outcome $Y_{a}$ corresponding to the value that the outcome
would have taken in a world in which the treatment were deterministically set
to $a$. The observed outcome is $Y=Y_{A}$.

###  2.2 Equivalence of Binarized Average Treatment Effect and Average
Treatment Effect Under Preservation of Relative Self-Selection

Presume that we are analyzing the binarized data $Y,T,W$. A common practice is
to regress $Y$ onto $T$ and $W$ to obtain an estimate of the regression
function $E[Y|T,W]$ and then report the coefficient on $T$. This is a specific
case of a more general plug-in estimation strategy where the estimand
(binarized ATE; BATE) is given by

$\psi_{\text{BATE}}=E[E[Y|T=1,W]]-E[E[Y|T=0,W]]$ (1)

In other words, we use the regression model to impute or predict both
“potential outcomes” for all subjects, then average across them all and take
the difference.

This estimand is represented in terms of the binarized treatment $T$. Our
question is: what causal contrast (if any) does this estimand represent in
terms of the original treatment $A$?

Consider the estimand $E[E[Y|T=1,W]]$. Using total expectation,

$\displaystyle E[E[Y|T=1,W]$ $\displaystyle=E[E[Y|A\in\mathcal{A},W]]$
$\displaystyle=E_{W}[E_{A|W,A\in\mathcal{A}}[E[Y|A,A\in\mathcal{A},W]]]$
$\displaystyle=E_{W}[E_{A|W,A\in\mathcal{A}}[E[Y|A,W]]]$
$\displaystyle=E[\mu(\tilde{A}_{1},W)]$

where we have defined $\mu(a,w)=E[Y|A=a,W=w]$ and a variable $\tilde{A}_{1}$
such that $\tilde{A}_{1}|W\sim A|W,A\in\mathcal{A}$ for every $W$. In terms of
densities, the distribution of $\tilde{A}_{1}|W$ follows directly from the
conditioning: for $a\in\mathcal{A}$, we have
$p_{A|W,A\in\mathcal{A}}(a,w)=p_{A|W}(a,w)/\pi_{\mathcal{A}}(w)$ where we have
defined $\pi_{\mathcal{A}}(w)=Pr(A\in\mathcal{A}|W=w)$. Thus, extending to the
support of $A$, the variable $\tilde{A}_{1}|W$ has density

$\displaystyle p_{\tilde{A}_{1}|W}(a,w)$
$\displaystyle=\frac{1_{\mathcal{A}}(a)p_{A|W}(a,w)}{\pi_{\mathcal{A}}(w)}.$

The distribution of $\tilde{A}_{1}$ preserves what we term relative self-
selection of treatment values within strata $W$ of the population (conditional
on the treatments being within the desired region). Specifically,
$\frac{p_{\tilde{A}_{1}|W}(a,w)}{p_{\tilde{A}_{1}|W}(a^{\prime},w)}=\frac{p_{A|W}(a,w)}{p_{A|W}(a^{\prime},w)}$
for any two treatment values $a,a^{\prime}\in\mathcal{A}$. We can therefore
think of $\tilde{A}_{1}$ as the unique treatment assignment mechanism
$\tilde{A}$ for which $P(\tilde{A}\in\mathcal{A})=1$ that also preserves the
relative self-selection preferences expressed by
$\frac{p_{A|W}(a,w)}{p_{A|W}(a^{\prime},w)}$.

This bears a little more explanation: consider a simple case with three
treatments $\\{1,2,3\\}$ and a population where all people have an equal
chance $1/3$ of receiving any of them. If we “outlaw” treatment 1 (i.e.
threshold treatment so $\mathcal{A}=\\{2,3\\}$), then under treatment policy
$\tilde{A}_{1}$, we would now have an equal $1/2$ chance of seeing treatments
2 or 3 (preserving the relative probabilities between them within strata of
the population).

Preservation of relative self-selection is akin to Luce’s definition of
independence of irrelevant alternatives (IIA) [21]. In the definition of IIA,
actors are typically considering choosing between only two options before the
introduction of irrelevant alternatives. However, in the BATE case, we can
compare the ratio of probability density of selecting one value of treatment
$a$ relative to another $a^{\prime}$ conditional on covariates given the
choices present in the smaller set $\mathcal{A}$ or the choices present in the
full support $\text{supp}(A)$. Under IIA, adding additional choices for
treatment should not change the relative self-selection ratio. Therefore,
interventions which act on individuals who follow IIA can be modeled by this
MTP.

Given this exposition, we now define $\tilde{A}_{0}$ such that
$\tilde{A}_{0}|W\sim A|A\in\mathcal{A}^{\mathsf{c}},W$. By identical arguments
as above, $E[E[Y|T=1,W]=E[\mu(\tilde{A}_{0},W)]$. Therefore,

$\psi_{\text{BATE}}=E[\mu(\tilde{A}_{1},W)]-E[\mu(\tilde{A}_{0},W)].$ (2)

###  2.3 Identification and Causal Assumptions

The reason we have cast the binarized ATE in the form given above is because
it has a very natural causal interpretation under standard causal assumptions.
Specifically, assume:

  1. 1.

Positivity: $0<Pr(A\in\mathcal{A}|W=w)<1$ for all $w\in\text{supp}(W)$.

  2. 2.

Conditional independence: $Y_{a}\perp\\!\\!\\!\\!\perp A\,|\,W$ for all $a\in$
supp($A$).

Consider an arbitrary modified treatment $\tilde{A}=\tilde{A}_{1}$ or
$\tilde{A}=\tilde{A}_{0}$. Under these assumptions we can write:

$\displaystyle E[\mu(\tilde{A},W)]$
$\displaystyle=E_{\tilde{A},W}[E[Y|A=a,W=w]]$
$\displaystyle=E_{\tilde{A},W}[E[Y_{A}|A=a,W=w]]\quad\text{(definition of
$Y$)}$
$\displaystyle=E_{\tilde{A},W}[E[Y_{a}|A=a,W=w]]\quad\text{($Y_{a}=Y_{A}$ when
$A=a$)}$
$\displaystyle=E_{\tilde{A},W}[E[Y_{a}|\tilde{A}=a,W=w]]\quad\text{(conditional
independence)}$ $\displaystyle=E_{\tilde{A},W}[E[Y_{\tilde{A}}|W=w]]$
$\displaystyle=E[Y_{\tilde{A}}]$

Note that the positivity assumption allows us to say that the expectations
$E[Y|\tilde{A}_{1},W]$ and $E[Y|\tilde{A}_{0},W]$ will be well-defined.
Moreover the conditional independence between $Y(a)$ and $A$ implies the same
for $\tilde{A}$.

Therefore, finally:

$\psi_{\text{BATE}}=E[Y_{\tilde{A}_{1}}]-E[Y_{\tilde{A}_{0}}]$ (3)

Under these assumptions, $\Psi_{BATE}$ is equivalent to the difference in
expected outcome under two modified treatment policies that affect the
underlying distribution of treatment in the population. Specifically, we can
consider setting $T=1$ equivalent to enacting a policy that restricts
$A\in\mathcal{A}$ while preserving relative self-selection, and setting $T=0$
equivalent to enacting a policy that restricts $A\in\mathcal{A}^{\mathsf{c}}$
while preserving relative self-selection. This interpretation is completely
transparent and there is actually no ambiguity in terms “different versions of
treatment”. Given the (standard) identifying assumptions, this parameter is a
perfectly reasonable causal contrast to use when estimating the effects of
policies under which subjects would have treatment probabilities in the
thresholded region proportional to what they would have been without the
thresholding. In the discussion we give examples of policies where this
assumption is sensible and when it is not.

###  2.4 New Parameter

In the previous section, we showed that the binarized ATE is equivalent to the
difference between two average counterfactual outcomes. However, when
considering the impact of a policy, it seems most natural to consider the
difference between a counterfactual outcome and the outcome under the status
quo with no policy, as suggested by Hubbard and van der Laan [9]. For example,
the estimand proposed by Haneuse and Rotnitzky is a difference that
“quantifies the average effect in the studied population of a policy that
switches each received dose…to the modified dose” [10].

To this end, we propose a new estimand of interest, the policy effect of
binarizing. It requires no data beyond what is used to calculate the binary
ATE. The policy effect of binarizing is the difference between the expected
outcome under one of the two policies and the expected outcome under the
status quo:

$\Psi_{PEB}=E[Y_{\tilde{A}_{t}}]-E[Y]\,\,\,\text{for }t\in\\{0,1\\}.$

Estimands of interest should correspond with the question researchers are
interested in exploring. Researchers are usually interested in asking causal
questions to help guide decisions when faced with a choice. For example, often
in medicine, researchers are interested in knowing whether a particular
treatment is effective in decreasing the probability of a bad outcome relative
to no treatment. Therefore, the ATE is an appropriate choice, as it compares a
counterfactual world where everyone in the target population receives the
treatment versus a counterfactual world in which no one in the target
population receives the treatment. In the case of testing a brand-new
treatment, the second counterfactual world essentially represents the status-
quo world, which currently exists without the treatment. In the case of
considering a binarized treatment variable, we argue that the relevant choices
are not between the two counterfactual worlds under two MTPs, but rather
between one counterfactual world and the status-quo world. Put in other words,
the relevant question is, “On average, how would the outcome change from the
status-quo if we were to force the treatment values into a particular interval
of allowable values?” $\Psi_{PEB}$ is the answer to this question.

In addition to answering a more relevant research question, $\Psi_{PEB}$
should be convenient for researchers to implement as identification of the
estimand requires the same set of assumptions to identify $\Psi_{BATE}$.
Therefore, data generated by processes that satisfied conditions allowing for
$\Psi_{BATE}$ to be estimated are also data in which $\Psi_{PEB}$ may be
estimated. In the following section, we describe two methods by which
researchers may implement estimators for $\Psi_{PEB}$, including details for
calculating standard errors.

###  2.5 Estimators

We consider three estimators of $\Psi_{BATE}$ and $\Psi_{PEB}$. The first is a
regression-based adjustment estimator [22]. The second is a standard inverse
probability weighting (IPW) estimator [23]. Both of these estimators are
commonly taught in introductory causal inference classes and represent two
classes of estimators that researchers might use to estimate causal effects
from observational data while adjusting for covariates. We finally consider
efficient, doubly-robust estimators. The augmented IPW (A-IPW) estimator is
discussed in the main body of the text, and a discussion of the targeted
maximum likelihood estiamtor (TMLE) is presented in Appendix A.

Under certain conditions, these estimators are consistent, meaning that as the
sample $n$ grows to infinity, the estimators will converge on the correct
population level values of the estimands. Specifically, the regression
estimator will be consistent if the outcome regression $E[Y|T=t,W=w]$ for
$t\in\\{0,1\\},w\in\text{supp}(W)$ is correctly specified. The IPW estimator
will be consistent if researchers have a correctly specified model for the
propensity score $\pi_{\mathcal{A}}(w)=Pr(T=1|W=w)$ for $w\in\text{supp}(W)$.
As previously mentioned, the A-IPW and TMLE estimators are doubly-robust. This
means that the estimators will be consistent estimators if researchers are
able to either consistently estimate the outcome regression $E[Y|T=t,W=w]$ or
the propensity score $\pi_{\mathcal{A}}(w)$ for
$t\in\\{0,1\\},w\in\text{supp}(W)$.

All four estimators are asymptotically normal, meaning that researchers can
calculate normal-based confidence intervals using the estimated standard
errors. We include details for calculating the standard errors for each
estimator below. In addition, if both the outcome regression $E[Y|T=t,W=w]$
and the propensity score $\pi_{\mathcal{A}}(w)$ are consistently estimated for
$T\in\\{0,1\\},w\in\text{supp}(W)$, then the A-IPW and TMLE estimators are
efficient, meaning that they have the smallest asymptotic variance.

As noted above, the four estimators rely upon estimators of conditional
expectations or conditional probabilities. Given that the MTPs defined above
result in treatment distributions that depend on the status-quo distribution
of $A$, one may assume that one must first estimate the density of $A$.
However, our equivalence result shows that the only relevant piece of the
distribution of $A$ that researchers need to estimate $\Psi_{BATE}$ and
$\Psi_{PEB}$ is the distribution of $T$, the binarized treatment. More
concretely, if researchers only had access to the data $(W,T,Y)$ for a sample
and no other information about the actual value of treatment $A$, they would
still be able to consistently estimate $\Psi_{BATE}$ and $\Psi_{PEB}$. Our
result also shows that, to estimate $\Psi_{BATE}$, one can use any standard
ATE estimation technique plugging in $T$ as the treatment variable. However,
the estimators for $\Psi_{PEB}$ take a different form, which we detail below.

####  2.5.1 Regression Estimator

The regression estimator $\widehat{\Psi}^{reg}$ first requires researchers to
calculate the following OLS coefficient estimates with Y as the outcome

$\hat{Y}=\hat{\beta_{0}}+T\hat{\beta}_{T}+\widecheck{W}^{\prime}\hat{\beta}_{W}+T\widecheck{W}^{\prime}\hat{\beta}_{interact}$

where $\widecheck{W}=W-\frac{1}{n}\sum_{i=1}^{n}W_{i}$ is the vector of
demeaned covariates. Let $\overline{T}=\frac{1}{n}\sum_{i=1}^{n}T_{i}$ and
$\overline{T\widecheck{W}^{\prime}}=\frac{1}{n}\sum_{i=1}^{n}T_{i}\widecheck{W}^{\prime}$
be the mean of $T$ and $T\widecheck{W}^{\prime}$ in the sample, respectively.
Then, the regression estimates for $\Psi_{BATE}$ and $\Psi_{PEB}$ are given by

$\displaystyle\widehat{\Psi}_{BATE}^{reg}$ $\displaystyle=\hat{\beta}_{T}$ (4)
$\displaystyle\widehat{\Psi}_{PEB}^{reg}$
$\displaystyle=\begin{cases}(1-\overline{T})\hat{\beta}_{T}-\overline{T\widecheck{W}^{\prime}}\hat{\beta}_{interact}\text{
for }t=1\\\
(\overline{T})\hat{\beta}_{T}-\overline{T\widecheck{W}^{\prime}}\hat{\beta}_{interact}\text{
for }t=0\\\ \end{cases}$ (5)

respectively.

In order to calculate the standard errors, we first employ M-estimation [24].
For $\widehat{\Psi}_{BATE}^{reg}$, one could also use a robust standard error
estimate for $\hat{\beta}_{T}$, such as the Eicker-Huber-White sandwich
variance estimator.

Let
$r_{i}=Y_{i}-(\hat{\beta_{0}}+T_{i}\hat{\beta}_{T}+\widecheck{W_{i}}^{\prime}\hat{\beta}_{W}+T_{i}\widecheck{W}_{i}^{\prime}\hat{\beta}_{interact})$
be the residual from the regression. First, note that
$\widehat{\theta}=(\overline{T},\overline{T\widecheck{W}},\hat{\beta_{0}},\hat{\beta}_{T},\hat{\beta}_{W},\hat{\beta}_{interact})^{\prime}$
is the vector of estimates that solves the moment condition

$\displaystyle\widehat{\mathbb{V}}(\hat{\beta})$
$\displaystyle=A_{n}(\hat{\beta})^{-1}B_{n}(\hat{\beta})(A_{n}(\hat{\beta})^{-1})^{\prime},$

where

$\displaystyle C_{i}$
$\displaystyle=\begin{bmatrix}1&T_{i}&\widecheck{W_{i}}^{\prime}&T_{i}\widecheck{W_{i}}^{\prime}\end{bmatrix}\begin{bmatrix}1\\\
T_{i}\\\ \widecheck{W_{i}}\\\ T_{i}\widecheck{W_{i}}\\\ \end{bmatrix}$
$\displaystyle A_{n}(\hat{\beta})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\begin{bmatrix}1&0&0\\\ 0&1&0\\\
0&0&C_{i}\\\ \end{bmatrix}$ $\displaystyle B_{n}(\hat{\beta})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\psi(O_{i},\widehat{\theta})\psi(O_{i},\widehat{\theta})^{\prime}$

Then, by delta method, we have that
$\sqrt{n}(\widehat{\Psi}_{BATE}^{reg}-\Psi_{BATE})$ and
$\sqrt{n}(\widehat{\Psi}_{PEB}^{reg}-\Psi_{PEB})$ are both asymptotically
normal with estimated variances given by

$\displaystyle\widehat{\mathbb{V}}\left(\widehat{\Psi}_{BATE}^{reg}\right)$
$\displaystyle=\begin{bmatrix}0&0&0&1&0&0\\\
\end{bmatrix}\mathbb{V}(\hat{\beta})\begin{bmatrix}0\\\ 0\\\ 0\\\ 1\\\ 0\\\
0\\\ \end{bmatrix}$ (7)
$\displaystyle\widehat{\mathbb{V}}\left(\widehat{\Psi}_{PEB}^{reg}\right)$
$\displaystyle=\begin{cases}\begin{bmatrix}-\hat{\beta_{T}}&-\hat{\beta}_{interact}&0&(1-\overline{T})&0&-\overline{T\widecheck{W}^{\prime}}\end{bmatrix}\mathbb{V}(\hat{\beta})\begin{bmatrix}-\hat{\beta_{T}}\\\
-\hat{\beta}_{interact}\\\ 0\\\ (1-\overline{T})\\\ 0\\\
-\overline{T\widecheck{W}}\end{bmatrix}\text{ for }t=1\\\
\begin{bmatrix}-\hat{\beta_{T}}&-\hat{\beta}_{interact}&0&-\overline{T}&0&-\overline{T\widecheck{W}^{\prime}}\end{bmatrix}\mathbb{V}(\hat{\beta})\begin{bmatrix}-\hat{\beta_{T}}\\\
-\hat{\beta}_{interact}\\\ 0\\\ -\overline{T}\\\ 0\\\
-\overline{T\widecheck{W}}\end{bmatrix}\text{ for }t=0\end{cases}$ (8)

####  2.5.2 Inverse Probability Weighting (IPW) Estimator

The IPW esitmator $\widehat{\Psi}^{IPW}$ requires an estimator of the
propensity score conditional on covariates $\pi_{\mathcal{A}}(w)$. Denote this
estimate $\widehat{\pi_{\mathcal{A}}}(w)$. Then, the IPW estimates for
$\Psi_{BATE}$ and $\Psi_{PEB}$ are given by

$\displaystyle\widehat{\Psi}_{BATE}^{IPW}$
$\displaystyle=\frac{1}{n}\left(\sum_{i=1}^{n}\frac{T_{i}Y_{i}}{\widehat{\pi_{\mathcal{A}}}(w)}-\sum_{i=1}^{n}\frac{(1-T_{i})Y_{i}}{1-\widehat{\pi_{\mathcal{A}}}(w)}\right)$
(9) $\displaystyle\widehat{\Psi}_{PEB}^{IPW}$
$\displaystyle=\left\\{\begin{array}[]{@{}l@{\quad}l@{}}\cfrac{1}{n}\left({\sum}_{i=1}^{n}\cfrac{T_{i}Y_{i}}{\widehat{\pi_{\mathcal{A}}}(w)}-{\sum}_{i=1}^{n}Y_{i}\right)\text{
for }t=1\\\
\cfrac{1}{n}\left({\sum}_{i=1}^{n}\cfrac{(1-T_{i})Y_{i}}{1-\widehat{\pi_{\mathcal{A}}}(w)}-{\sum}_{i=1}^{n}Y_{i}\right)\text{
for }t=0\\\ \end{array}\right.\kern-1.2pt$ (12)

respectively. The IPW estimator is asymptotically linear and therefore
converges in distribution to a normal distribution [25].

The IPW estimator has known influence function; however, calculating the
standard error by estimating the standard deviation of the influence function
results in overly-conservative inference [26]. Therefore, we suggest that
researchers calculate the standard errors for the IPW estimators described
here via the nonparametric bootstrap.

####  2.5.3 Augmented IPW Estimator

The A-IPW estimator $\widehat{\Psi}^{AIPW}$ requires both an estimator of the
conditional mean $m(t,w)\equiv E[Y|t,w]$ and an estimator of the propensity
score conditional on covariates $\pi_{\mathcal{A}}(w)$. Denote these estimates
$\widehat{m}(t,w)$ and $\widehat{\pi_{\mathcal{A}}}(w)$, respectively.
Researchers can use a variety of machine learning (ML) algorithms to create
these estimates. Super learning can also be used to determine the optimal
predictor (or ensemble of predictors) from a class of many candidate
estimators based on the data [27]. Researchers should employ sample splitting
and cross fitting to reduce over-fitting and ensure that the estimates
converge to the truth sufficiently quickly [28].

The efficient influence curves (EIFs) for $\Psi_{BATE}$ and $\Psi_{PEB}$ are
given by

$\displaystyle\phi(O)_{BATE}$
$\displaystyle=\left(\frac{T}{{\pi_{\mathcal{A}}}(w)}-\frac{(1-T)}{1-{\pi_{\mathcal{A}}}(w)}\right)(Y-{m}(T,W))+{m}(1,W)-{m}(0,W)-\Psi_{BATE}$
$\displaystyle\phi(O)_{PEB}$
$\displaystyle=\left\\{\begin{array}[]{@{}l@{\quad}l@{}}\cfrac{T}{{\pi_{\mathcal{A}}}(w)}(Y-{m}(1,W))+{m}(1,W)-Y-\Psi_{PEB}\text{
for }t=1\\\
\cfrac{(1-T)}{1-{\pi_{\mathcal{A}}}(w)}(Y-{m}(0,W))+{m}(0,W)-Y-\Psi_{PEB}\text{
for }t=0\\\ \end{array}\right.\kern-1.2pt$

respectively.

These EIFs motivate the A-IPW estimators, which for $\Psi_{BATE}$ and
$\Psi_{PEB}$ are given by

$\displaystyle\widehat{\Psi}_{BATE}^{AIPW}$
$\displaystyle=\frac{1}{n}\left[\sum_{i=1}^{n}\left(\frac{T_{i}}{\widehat{\pi_{\mathcal{A}}}(w)}-\frac{(1-T_{i})}{1-\widehat{\pi_{\mathcal{A}}}(w)}\right)(Y_{i}-\widehat{m}(T_{i},W_{i}))+\widehat{m}(1,W_{i})-\widehat{m}(0,W_{i})\right]$
(13) $\displaystyle\widehat{\Psi}_{PEB}^{AIPW}$
$\displaystyle=\left\\{\begin{array}[]{@{}l@{\quad}l@{}}\cfrac{1}{n}\left[{\sum}_{i=1}^{n}\cfrac{T_{i}}{\widehat{\pi_{\mathcal{A}}}(w)}(Y_{i}-\widehat{m}(1,W_{i}))+\widehat{m}(1,W)-Y_{i}\right]\text{
for }t=1\\\
\cfrac{1}{n}\left[{\sum}_{i=1}^{n}\cfrac{(1-T_{i})}{1-\widehat{\pi_{\mathcal{A}}}(w)}(Y_{i}-\widehat{m}(0,W_{i}))+\widehat{m}(0,W_{i})-Y_{i}\right]\text{
for }t=0\\\ \end{array}\right.\kern-1.2pt$ (16)

respectively.

The variance of the two estimators can be calculated either via estimation of
the influence curve or via nonparametric bootstrap.

The influence curve evaluated at a given observation $O_{i}$ for
$\widehat{\Psi}_{BATE}^{AIPW}$ and $\widehat{\Psi}_{PEB}^{AIPW}$ can be
estimated as

$\displaystyle\widehat{\phi}(O_{i})_{BATE}$
$\displaystyle=\left(\frac{T_{i}}{\widehat{\pi_{\mathcal{A}}}(w)}-\frac{(1-T_{i})}{1-\widehat{\pi_{\mathcal{A}}}(w)}\right)(Y-\widehat{m}(T_{i},W_{i}))+\widehat{m}(1,W_{i})-\widehat{m}(0,W_{i})-\widehat{\Psi}_{BATE}$
$\displaystyle\widehat{\phi}(O_{i})_{PEB}$
$\displaystyle=\left\\{\begin{array}[]{@{}l@{\quad}l@{}}\cfrac{T_{i}}{\widehat{\pi_{\mathcal{A}}}(w)}(Y_{i}-\widehat{m}(1,W_{i}))+\widehat{m}(1,W)-Y_{i}-\widehat{\Psi}_{PEB}^{AIPW}\text{
for }t=1\\\
\cfrac{(1-T_{i})}{1-\widehat{\pi_{\mathcal{A}}}(w)}(Y_{i}-\widehat{m}(0,W_{i}))+\widehat{m}(0,W_{i})-Y_{i}-\widehat{\Psi}_{PEB}^{AIPW}\text{
for }t=0.\\\ \end{array}\right.\kern-1.2pt$

Then, the variance estimates can be calculated as

$\displaystyle\widehat{\mathbb{V}}\left(\widehat{\Psi}_{BATE}^{AIPW}\right)$
$\displaystyle=\frac{1}{n}{\sum}_{i=1}^{n}\widehat{\phi}(O_{i})_{BATE}^{2}$
$\displaystyle\widehat{\mathbb{V}}\left(\widehat{\Psi}_{PEB}^{AIPW}\right)$
$\displaystyle=\frac{1}{n}{\sum}_{i=1}^{n}\widehat{\phi}(O_{i})_{PEB}^{2}$

Note that, while the estimator is doubly-robust, the estimator for the
variance is not - the influence curve based inference is only consistent when
both $\widehat{m}(t,w)$ and $\widehat{\pi_{\mathcal{A}}}(a)$ are consistent
estimators [25].

##  3 Simulation Study

###  3.1 Data Generating Processes

To illustrate the equivalence between the binarized ATE and difference in
expected outcomes under the proposed MTPs, we performed a simulation study.

$\displaystyle W$ $\displaystyle\sim\text{ Bern}(p=0.5)$ $\displaystyle A|W$
$\displaystyle\sim\mathcal{N}(5+2*W,1)$ $\displaystyle Y|A,W$
$\displaystyle\sim A^{3}+\text{sin}(A)+100*W+\mathcal{N}(0,1)$ $\displaystyle
T|A$ $\displaystyle=\begin{cases}1\text{ when }A\geq 6\\\ 0\text{ when }A<6\\\
\end{cases}$

The estimands of interest are the binarized ATE and the policy effect of
binarization when implementing a cut off at a treatment value of $A\geq 6$.
Here, we define the policy effect of binarization as
$\Psi_{PEB}=E[Y_{\tilde{A}_{1}}]-E[Y]$, which contrasts the average outcome
under the MTP which restricts values $A\geq 6$ and preserves relative self-
selection with the average outcome in the status-quo world. Under the data
generating process, the true parameter values are $\Psi_{BATE}=201.806$ and
$\Psi_{PEB}=90.872$.

Over 2000 simulations, we drew samples of size $n=150,300,500$ from the data
generating process and calculated the estimates of $\Psi_{BATE}$ and
$\Psi_{PEB}$ using the two different estimators. The regression estimator is
consistent for this DGP because $E[Y|T,W]$ is linear in $T$ and $W$, as $W$ is
binary. For a more general discussion of when the regression estimator will be
consistent, see Appendix B. The estimate $\hat{\pi}_{\mathcal{A}}(W)$
necessary for implementing the IPW estimator was calculated using logistic
regression with $W$ as the predictor and $T$ as the binary outcome, which is
consistent for the propensity score for this DGP. Therefore, both estimators
are consistent. We also estimated the standard errors of the estimates using
the methods described in 2.5. Finally, we present the percent bias of the
estimates as well as the simulated standard error to evaluate the performance
of the estimators.

###  3.2 Visualizing the MTPs

1 illustrates the resulting densities after implementing the cut-off for the
two values of covariates. The light blue line represents the cut-off value of
6. The green line is the status-quo density $A$, the red dotted line is the
density of $\tilde{A}_{1}$, and the blue dotted line is the density of
$\tilde{A}_{0}$. As illustrated by this figure, the densities $\tilde{A}_{1}$
and $\tilde{A}_{0}$ differ depending on the value of the covariate $W$. This
is generally true of our result; the resulting distributions of treatment
after binarizing can depend both on the status-quo distribution of treatment
as well as the pre-treatment covariates.

Figure 1: Density of treatment $A$ by different values of $W$ in the status-
quo world and after implementing the proposed modified treatment policies.

Both $\tilde{A}_{1}$ and $\tilde{A}_{0}$ have density equal to 0 for values of
treatment that are not within the allowed range ($A<6$ for $\tilde{A}_{1}$ and
$A\geq 6$ for $\tilde{A}_{0}$). It is also evident that both preserve relative
self-selection for the values of treatment within the allowed range. The non-
zero portions of the MTP distributions of treatment are simply scaled versions
of the density of the status-quo $A$.

###  3.3 Results

| Regression Estimator | IPW Estimator
---|---|---
n | Estimate | Bias | Est SE | Sim SE | Estimate | Bias | Est SE | Sim SE
150 | 201.840 | 0.515 | 14.702 | 14.807 | 202.479 | 1.155 | 14.775 | 14.629
300 | 202.243 | 0.918 | 10.072 | 10.160 | 202.325 | 1.000 | 10.096 | 10.207
500 | 202.580 | 1.256 | 7.805 | 7.957 | 202.274 | 0.950 | 7.842 | 7.734
Table 1: Simulation results for estimating $\Psi_{BATE}$. The true value of the estimand is 201.8065. | Regression Estimator | IPW Estimator
---|---|---
n | Estimate | Bias | Est SE | Sim SE | Estimate | Bias | Est SE | Sim SE
150 | 89.580 | -1.362 | 12.331 | 11.447 | 90.299 | -0.573 | 11.561 | 11.698
300 | 89.927 | -0.945 | 8.719 | 7.912 | 89.940 | -0.932 | 7.903 | 7.847
500 | 90.141 | -0.731 | 6.830 | 6.283 | 90.026 | -0.846 | 6.146 | 6.216
Table 2: Simulation results for estimating $\Psi_{PEB}$. The true value of the
estimand 90.872.

Using the methods described in 2.5, we calculated estimates for the target
parameters and the standard errors of the estimates. What bears repeating is
that the only information the estimators require about the treatment $A$ is
the value of the binarized treatment $T$. One may assume that, in order to
estimate the expected outcomes under the MTPs, one must first estimate the
full densities of $\tilde{A}_{1}$ and $\tilde{A}_{0}$. However, through the
equivalence with the binarized ATE, the relevant quantity one must estimate is
instead the conditional probability $Pr(T=1|W=w)=Pr(A\geq 6|W=w)$ for
$w\in\\{0,1\\}$, which is easier to estimate.

As the results in tables 1 and 2 show, both the regression and the IPW
estimators for $\Psi_{PEB}$ and $\Psi_{PATE}$ have almost negligible bias,
even for the smaller sample size of $n=150$. When estimating both $\Psi_{PEB}$
and $\Psi_{PATE}$, the regression and IPW estimators have similar standard
errors, and thus would result in relatively similar normal-based confidence
intervals.

##  4 Discussion

In this paper, we show the equivalence between the binarized ATE and the
difference between the average counterfactual outcomes under two MTPs that
impose the relevant cut-off defined by the binarization and preserve relative
self-selection. Relative self-selection is the ratio of the probability
density of one value of treatment to another. Through this equivalence, we
affirm that the binarized ATE is a valid causal estimand and is identified
under standard causal identification assumptions, addressing many of the
common concerns methods researchers have positted about such coarsening of a
continuous treatment variable.

We also introduce a new estimand, the policy effect of binarization, which
contrasts the expected outcome under one of the MTPs associated with
binarization with the expected outcome in the status-quo world. This estimand
requires the same identification assumption as the binarized ATE and therefore
can be estimated from the same data researchers may use to estimate the
binarized ATE. We argue that this estimand answers a more relevant causal
question as it compares the counterfactual world under one MTP (the “treated”
counterfactual) to the status-quo world (which is akin to a “control”
counterfactual). As this estimand cannot be calculated using standard methods
developed for the ATE, we also provide details for how to implement three
different estimators of the policy effect of binarization as well as how to
estimate their standard errors.

Finally, we conducted a simulation study to show analytically that our results
hold using the estimators that we proposed. Specifically, using only the
binarized variable $T$ based on the cut-off and no other information about the
status-quo treatment distribution of $A$, we were able to provide consistent
estimates for both $\Psi_{BATE}$ and $\Psi_{PEB}$. This result may be
surprising, as typically estimating the effects of MTPs requires estimating
the density of $A$. For example, this is the case in the shift intervention
considered by Díaz and van der Laan [29]. Therefore, the two MTPs discussed
which impose a cut-off and preserve relative self-selection lend themselves to
much easier estimation from a practioner’s point of view, which perhaps
explains part of the reason why the method is so prevalent in the literature
across disciplines. In fact, it has been shown that for all MTPs which impose
a cut-off and restrict treatments to a certain region of the support, the
causal effects for the unique MTP which preserves relative self-selection are
indeed the easiest to estimate from a theoretical view point [19].

However, although the binarized ATE and the policy effect of binarization may
be easy to estimate, they may not be relevant causal quantities depending on
the specific research context. Ultimately, whether the binarized ATE or the
policy effect of binarization should be the estimand of interest depends on
the expected behavior of individuals after implementing the intervention of
interest. For example, suppose researchers were interested in estimating the
effect of a policy that limits the amount of energy households can use by
fining households that consume more energy than a certain threshold. In this
case, the assumption of preservation of relative self-selection is not very
plausible, as one could reason that households above the threshold would
change their energy usage to just barely meet the new limits, as there is no
incentive to reduce energy consumption any further. Conversely, households
that already use less energy than the threshold will likely make no change to
their energy use. Therefore, the resulting distribution of household energy
usage would have a “piling up” of households at the imposed threshold with the
distribution below the threshold looking the same as status-quo. Therefore, in
this case, it would likely not make sense to binarize household energy usage
based on the proposed cut-off.

On the other hand, suppose researchers were interested in estimating the
effect of the distance an individual lives from a manufacturing plant.
Specifically, the researchers are interested in estimating the effect of a
proposed policy which would give individuals who currently live within a cut-
off distance a housing voucher which they can use to relocate throughout the
city. In this case, it might be more reasonable that preservation of relative
self-selection would hold after implementing the policy. Deciding where to
live is a multifacted decision based on individual needs that may depend on
different factors such as access to public transportation, quality of local
schools, and proximity to grocery stores, among other factors. The individuals
who chose to live close to the factory may have done so based on different
needs; therefore, it seems plausible that these individuals would relocate
evenly across the other populated areas of the city if given the economic
assistance to do so.

These are simply two examples of policies researchers may be interested in.
Ultimately, the decision to binarize should depend on domain-specific
knowledge and predictions regarding the effect of the relevant policy on the
distribution of the treatment variable in the population. Unfortunately,
before the policy is implemented, this is a completely untestable assumption.
However, researchers can still make plausible educated guesses and choose a
causal estimand of interest accordingly.

## References

  * [1] Neyman’s Repeated Sampling Approach to Completely Randomized Experiments.  In Donald B. Rubin and Guido W. Imbens, editors, Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction, pages 83–112. Cambridge University Press, Cambridge, 2015.
  * [2] Michal Shimonovich, Anna Pearce, Hilary Thomson, Gerry McCartney, and Srinivasa Vittal Katikireddi.  Assessing the causal relationship between income inequality and mortality and self-rated health: protocol for systematic review and meta-analysis.  Systematic Reviews, 11(1):20, February 2022.
  * [3] Hang Zhou and Vasilis Vasiliou.  Alcohol Use and Use Disorder and Cancer Risk: Perspective on Causal Inference.  Complex Psychiatry, 8(1-2):9–12, August 2022.
  * [4] Ruifa Hu, Xusheng Huang, Jikun Huang, Yifan Li, Chao Zhang, Yanhong Yin, Zhaohui Chen, Yanhong Jin, Jinyang Cai, and Fang Cui.  Long- and Short-Term Health Effects of Pesticide Exposure: A Cohort Study from China.  PLoS ONE, 10(6):e0128766, June 2015.
  * [5] Frank B. Hu, JoAnn E. Manson, Meir J. Stampfer, Graham Colditz, Simin Liu, Caren G. Solomon, and Walter C. Willett.  Diet, Lifestyle, and the Risk of Type 2 Diabetes Mellitus in Women.  New England Journal of Medicine, 345(11):790–797, September 2001\.  Publisher: Massachusetts Medical Society.
  * [6] Xiao Wu, Danielle Braun, Marianthi-Anna Kioumourtzoglou, Christine Choirat, Qian Di, and Francesca Dominici.  Causal inference in the context of an error prone exposure: Air pollution and mortality.  The Annals of Applied Statistics, 13(1):520–547, March 2019.  Publisher: Institute of Mathematical Statistics.
  * [7] Marcus J. Kurtz and Adam Lauretig.  Does Free-Market Reform Induce Protest? Selection, Post-Treatment Bias, and Depoliticization.  British Journal of Political Science, 52(2):968–976, April 2022\.  Publisher: Cambridge University Press.
  * [8] Pedro Carneiro, Michael Lokshin, and Nithin Umapathi.  Average and Marginal Returns to Upper Secondary Schooling in Indonesia.  Journal of Applied Econometrics, 32(1):16–36, 2017.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/jae.2523.
  * [9] Alan E. Hubbard and Mark J. van der Laan.  Population intervention models in causal inference.  Biometrika, 95(1):35–47, 2008.
  * [10] Sebastien Haneuse and Andrea Rotnitzky.  Estimation of the effect of interventions that modify the received treatment.  Statistics in Medicine, 32(30):5260–5277, December 2013.
  * [11] Kenneth J. Rothman, Sander Greenland, and Timothy L. Lash.  Modern Epidemiology.  Wolters Kluwer, Philadelphia, United States, 2008.
  * [12] Tyler J. VanderWeele and Miguel A. Hernán.  Causal Inference Under Multiple Versions of Treatment.  Journal of causal inference, 1(1):1–20, May 2013.
  * [13] Tyler J. VanderWeele, Yu Chen, and Habibul Ahsan.  Inference for causal interactions for continuous exposures under dichotomization.  Biometrics, 67(4):1414–1421, December 2011.
  * [14] Caroline Bennette and Andrew Vickers.  Against quantiles: Categorization of continuous variables in epidemiologic research, and its discontents.  BMC Medical Research Methodology, 12(1):21, February 2012.
  * [15] Ori M Stitelman, Alan E Hubbard, and Nicholas P Jewell.  The Impact Of Coarsening The Explanatory Variable Of Interest In Making Causal Inferences: Implicit Assumptions Behind Dichotomizing Variables.  2010\.
  * [16] Iván Díaz Muñoz and Mark van der Laan.  Population Intervention Causal Effects Based on Stochastic Interventions.  Biometrics, 68(2):541–549, 2012.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1541-0420.2011.01685.x.
  * [17] David McCoy, Alan Hubbard, Alejandro Schuler, and Mark van der Laan.  Cross-validated decision trees with targeted maximum likelihood estimation for nonparametric causal mixtures analysis, 2023.
  * [18] Evan K. Rose and Yotam Shem-Tov.  On Recoding Ordered Treatments as Binary Indicators, July 2023\.  arXiv:2111.12258 [econ, stat].
  * [19] Lars van der Laan, Wenbo Zhang, and Peter B. Gilbert.  Nonparametric estimation of the causal effect of a stochastic threshold-based intervention.  Biometrics, 79(2):1014–1028, June 2023.
  * [20] Donald Rubin.  Estimating Causal Effects of Treatments in Experimental and Observational Studies.  ETS Research Bulletin Series, 1972(2):i–31, 1972.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/j.2333-8504.1972.tb00631.x.
  * [21] Paramesh Ray.  Independence of Irrelevant Alternatives.  Econometrica, 41(5):987–991, 1973.  Publisher: [Wiley, Econometric Society].
  * [22] Peng Ding.  A First Course in Causal Inference, October 2023.  arXiv:2305.18793 [stat].
  * [23] James M. Robins, Andrea Rotnitzky, and Lue Ping Zhao.  Estimation of Regression Coefficients When Some Regressors are not Always Observed.  Journal of the American Statistical Association, 89(427):846–866, September 1994.
  * [24] Leonard A. Stefanski and Dennis D. Boos.  The Calculus of M-Estimation.  The American Statistician, 56(1):29–38, 2002.  Publisher: [American Statistical Association, Taylor & Francis, Ltd.].
  * [25] Mark J. van der Laan and Jamies M. Robins.  Unified Methods for Censored Longitudinal Data and Causality.  Springer Series in Statistics. Springer, 2003.
  * [26] Jared K. Lunceford and Marie Davidian.  Stratification and weighting via the propensity score in estimation of causal treatment effects: a comparative study.  Statistics in Medicine, 23(19):2937–2960, 2004.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/sim.1903.
  * [27] Rachael V Phillips, Mark J van der Laan, Hana Lee, and Susan Gruber.  Practical considerations for specifying a super learner.  International Journal of Epidemiology, 52(4):1276–1285, August 2023\.
  * [28] 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, February 2018.
  * [29] Iván Díaz Muñoz and Mark van der Laan.  Population Intervention Causal Effects Based on Stochastic Interventions.  Biometrics, 68(2):541–549, 2012.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1541-0420.2011.01685.x.
  * [30] Megan S. Schuler and Sherri Rose.  Targeted Maximum Likelihood Estimation for Causal Inference in Observational Studies.  American Journal of Epidemiology, 185(1):65–73, January 2017.
  * [31] Susan Gruber and Mark J. van der Laan.  tmle: An R package for targeted maximum likelihood estimation.  Journal of Statistical Software, 51(13):1–35, 2012.  doi:10.18637/jss.v051.i13.

##  Appendix A Targeted Maxmimum Likelihood Estimator (TMLE)

The TMLE estimator $\widehat{\Psi}^{TMLE}$ is an efficient, doubly-robust
estimator. We provide a brief overview of the steps required below. For a more
comprehensive introduction to TMLE for applied work in the binary treatment
setting, we refer researchers to the 2017 paper “Targeted Maximum Likelihood
Estimation for Causal Inference in Observational Studies” by Schuler and Rose
[30]. The TMLE presented for the BATE is equivalent to the binTMLE presented
in [19] and is efficient under the assumption of missingness at random.

  1. 1.

Estimate the conditional expectation function $E[Y|T=t,W=w]$ for
$t\in\\{0,1\\},w\in\text{supp}(W)$. Use this estimator to generate initial
estimates $\widehat{Y}_{\tilde{A}_{1}}$ and $\widehat{Y}_{\tilde{A}_{0}}$ for
each individual in the sample, which are the predicted expectations when $T=1$
and $T=0$, respectively. Let
$\hat{Y}=1(T=1)\hat{Y}_{\tilde{A}_{1}}+1(T=0)\hat{Y}_{\tilde{A}_{0}}$ be the
prediction for the observed outcome.

  2. 2.

Estimate the propensity score function $Pr(T=1|W=w)$ for $w\in\text{supp}(W)$.
Denote this estimator $g_{0}(w)$.

  3. 3.

Calculate the “clever covariate” $H(T=t,W=w)$, which differs for the two
estimators, for each individual in the sample.

     * •

If estimating $\Psi_{BATE}$,
$H=\frac{1(T=1)}{g_{0}(W)}-\frac{1(T=0)}{1-g_{0}(W)}$

     * •

If estimating $\Psi_{PEB}$, $H=1-\frac{1(T=1)}{g_{0}(W)}$ when considering
$E[Y_{\tilde{A}_{1}}]$ and $H=1-\frac{1(T=0)}{1-g_{0}(W)}$ when considering
$E[Y_{\tilde{A}_{0}}]$.

Then, use logistic regression to estimate the coefficient $\beta$ regressing H
on the observered outcome $Y$ using $\hat{Y}$ as a fixed intercept:

$\text{logit}(Y)=\text{logit}(\hat{Y})+\beta H(T,W)$

Use this estimated coefficient $\hat{\beta}$ to calculate new, “updated”
predictions for the potential outcomes $\hat{Y}_{\tilde{A}_{1}}^{*}$ and
$\hat{Y}_{\tilde{A}_{0}}^{*}$ for each individual in the sample:

$\displaystyle\text{logit}(\hat{Y}_{\tilde{A}_{1}}^{*})$
$\displaystyle=\text{logit}(\widehat{Y}_{\tilde{A}_{1}})+\hat{\beta}\frac{1}{g_{0}(W)}$
$\displaystyle\text{logit}(\hat{Y}_{\tilde{A}_{0}}^{*})$
$\displaystyle=\text{logit}(\widehat{Y}_{\tilde{A}_{0}})+\hat{\beta}\frac{1}{1-g_{0}(W)}$

  4. 4.

Let
$\hat{Y}^{*}=1(T=1)\hat{Y}_{\tilde{A}_{1}}^{*}+1(T=0)\hat{Y}_{\tilde{A}_{0}}^{*}$
be the new prediction for the observed outcome. Now, we can calculate

$\displaystyle\hat{\Psi}_{BATE}^{TMLE}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\hat{Y}_{\tilde{A}_{1}}^{*}-\hat{Y}_{\tilde{A}_{0}}^{*}$
(17) $\displaystyle\hat{\Psi}_{PEB}^{TMLE}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\hat{Y}^{*}-\hat{Y}_{\tilde{A}_{t}}^{*}$
(18)

To estimate the variance, researchers can either calculate influence curve
based estimates or use the nonparametric bootstrap. The influence curve based
estimates are the same as the AIPW estimator as in Section 2.5.3 and have the
same consistently conditions.

If estimating $\Psi_{BATE}$, researchers can use the R package tmle, using $T$
as the treatment variable [31]. The package will calculate the estimate as
well as the standard errors based on the influence curve.

##  Appendix B When is the regression estimator consistent?

Suppose the expectation of $Y$ given $A$ and $W$ is linear in $W$ and $W$
times some function of $A$. In order for the regression estimator to be
consistent, the expectation of $Y$ given $T$ and $W$ must be linear. When will
this be the case? Let $g_{A}(A)$ and $g_{int}(A)$ be two (potentially non-
linear) functions of $A$.

Consider a general DGP where the expectation of $Y$ given $A$ and $W$ is a
linear function in $W$, $g_{A}(A)$, and in an interaction term $W*g_{int}(A)$,
i.e.

$E[Y|A,W]=\beta_{0}+\beta_{W}W+\beta_{A}g_{A}(A)+\beta_{int}W*g_{int}(A).$

Consider a cut-off at $A=a_{0}$ such that $T=1$ when $A\geq a_{0}$ and $T=0$
when $A<a_{0}$. Then, we can consider the expectation of $Y$ given $T=1$ and
$W$.

$\displaystyle E[Y|T=1,W]$
$\displaystyle=\frac{1}{\pi_{\mathcal{A}}(w)}\int_{a_{0}}^{\infty}E[Y|A,W]da$
$\displaystyle=\frac{1}{\pi_{\mathcal{A}}(w)}\int_{a_{0}}^{\infty}\beta_{0}+\beta_{W}W+\beta_{A}g_{A}(a)+\beta_{int}W*g_{int}(a)da$
$\displaystyle=\frac{1}{\pi_{\mathcal{A}}(w)}\left(\beta_{0}+\beta_{W}W\right)+\frac{\beta_{A}}{\pi_{\mathcal{A}}(w)}\int_{a_{0}}^{\infty}g_{A}(a)da+\beta_{int}*W\frac{1}{\pi_{\mathcal{A}}(w)}\int_{a_{0}}^{\infty}g_{int}(a)da$
$\displaystyle=\frac{1}{\pi_{\mathcal{A}}(w)}\left(\beta_{0}+\beta_{W}W\right)+\beta_{A}E[g_{A}(\tilde{A}_{1})|W]+\beta_{int}WE[g_{int}(\tilde{A}_{1})|W]$

Similarly for $E[Y|T=0,W]$, we have

$E[Y|T=0,W]=\frac{1}{1-\pi_{\mathcal{A}}(w)}\left(\beta_{0}+\beta_{W}W\right)+\beta_{A}E[g_{A}(\tilde{A}_{0})|W]+\beta_{int}WE[g_{int}(\tilde{A}_{0})|W]$

Therefore, the regression estimator for $\Psi_{BATE}$ will be consistent when
$E[g_{A}(\tilde{A}_{t})|W]$ and $E[g_{int}(\tilde{A}_{t})|W]$ are linear in
$W$ for $t\in\\{0,1\\}$. The regression estimator for $\Psi_{PEB}$ when
considering $T=t$ will be consistent when $E[g_{A}(\tilde{A}_{t})|W]$ and
$E[g_{int}(\tilde{A}_{t})|W]$ are linear in $W$.