## Methodological Foundations of Modern Causal Inference in Social Science Research

By Guanghui Pan ∗

## I. Notations, and Fundamental Analytical Framework

Inference is the core of statistical analysis. Compared to descriptive analysis of samples, statisticians are more interested in observing a small part of the individuals in a population. Through appropriate inductive reasoning, we can gain knowledge about the population. Therefore, they propose hypotheses and design mathematical models, whether simple or complex, aimed at inferring the characteristics of the population.

Mathematically, the population is abstracted as a closed set. Some attributes of the population are measurable . We define measurable sets Z as a set of non-empty, closed (under complement, countable unions, and countable intersections) subsets of the population; 1 .If we have a set of functions P measuring (assign numbers to the characteristics of) the measurable set, or mathematically, mapping the measurable set onto the real numbers ( P : Z → R ), we call the functions as measures , or, data-generating process (DGP) .In other words, suppose in a measurable set Z we have two measurable attributes, Z = ( X,Y ); the data-generating process gives out the measurable distributions X and Y .

However, the DGP is, at most times, hard to describe. Consider X to represent the heights and Y to represent the weights of individuals within the British population. It is hard to decipher the data-generating distribution for heights and weights, even if we could possibly obtain them from the census data. An intuitive way to describe it is to set specific statistical target parameters to describe it, for instance, the mean and variance of the height and weight of the British people. The function to map the DGP to the target parameter is called estimandor ψ , and the number of the target parameter ψ ( P ) is usually called the (statistical) estimand (i.e., ψ ( P ) = E P ( Z )).

A more common condition is that we could never know the true DGP P ; we

∗ Department of Sociology, University of Oxford. 41 Park End Street, Oxford, the United Kingdom, OX1 1HH. Email: guanghui.pan@sociology.ox.ac.uk. Manuscript in progress, do not cite or redistribute without permission.

1 which is called σ -algebra. Indeed, suppose the population is D and the measurable sets Z , we define a measurable space with the pair ( D , Z ). Further, with measure P , we define the triple ( D , Z, P ) as the measure space .

could only observe the samples Z i = ( X , Y i i ). Statistical inference is the process by which we use the samples to speculate the characteristics of the population. In this thesis, we assume the samples we used in statistical inference are randomly drawn from the population through the same process, or in other words 2 :

ASSUMPTION I.1 (Independent and Identically Distributed): The samples are independent and identically distributed . 3

With IID samples, we have the empirical measure (distribution) based on the observations P n ( n denotes the sample size). Correspondingly, we have an estimator denoted as ˆ and the empirical ψ point estimate ψ ( P n ). 4 If ψ ( P n ) contains the sample mean, or any linear combination regarding the sample means, under the IID assumption, we have the central limit theorem (CLT):

THEOREM I.1: Central Limit Theorem (CLT): suppose µ and σ 2 separately denote the expectation and the variance for the IID observations. If σ 2 is finite, as the sample size n approaches infinity, difference between the average sample mean ˆ ψ and the expectation µ approaches a normal distribution with mean 0 and variance σ 2 , at the rate of 1 / √ n :

$$\sqrt { n } ( \psi ( \mathbb { P } _ { n } ) - \mu ) \to \mathcal { N } ( 0, \sigma ^ { 2 } )$$

If µ is an unbiased estimator for ψ ( P ), Equation I.1 describes the asymptotic relationship between ˆ and ψ ψ , and the difference between ˆ and ψ ψ is called statistical error . I will elaborate more on the analysis of statistical error in Section III.

Indeed, in some statistical analyses, our ultimate goal is not to estimate the statistical estimand. Instead, we use statistical estimand to approach the 'real' estimand in our problems. Consider the following common scenarios in sociological or demographic research:

- 1) Suppose we are interested in the causal relations between a treatment and the outcome. Assume that our treatment is binary with only two values:

2 Rigorously, the DGP for random samples Z i are called, in most textbooks, the random variables. Therefore, the random variable is indeed a function and usually is denoted as Z = Z i .

3 For instance, using convenience sampling does not draw IID samples. We cannot identify each sample's probability of being sampled from the whole population, and the samples are not connected (consider an extreme scenario in which we sample the heights from a British professional basketball league that can never represent the heights of the British population). Nonetheless, samples from nonIID settings may also perform statistical inference, but I will not address the techniques in this thesis.

4 We could understand the estimandor and the estimated separately as the function and the quantity to be estimated, while the estimator and the estimate separately as the function and the quantity we perform estimation based on the observed data. However, in most statistical papers and textbooks, estimator and estimate are interchangeable and rarely rigorously defined, and most researchers assume that the estimator and the estimador have the same functions ψ . In this thesis, if the notation is clear, I also use ψ to refer to the estimand ψ ( P ) and ˆ to refer to the estimator ψ ψ ( P n )).

treatment and control. Also, the treatment should be assigned before the outcomes being observed. We use Y (1) to denote the outcome under treatment, while Y (0) to denote the outcome under control (in some textbooks, they are denoted as Y t and Y c ). Let A denote the treatment variable, and X denote the set of covariates. The causal measurable set can be written as Z ∗ = ( X,A,Y (1) , Y (0)) under the causal DGP P ∗ . Suppose our target of interest is the average treatment effect defined as the difference between the expectation of Y (1) and Y (0): ψ ∗ ( P ∗ ) = E P ∗ [ E Y [ (1)] -E Y [ (0)] ] . In the real world, due to the fundamental problem in causal inference, we can not observe the outcomes under the treatment and control simultaneously; we could only have the measurable set Z = ( X,A,Y ) under the statistical DGP P and try to use a reasonable statistical estimand to approach ψ ∗ , for instance, the conditional expected outcome given the treatment is assigned versus the control is assigned: ψ ( P ) = E P [ E Y A [ | = 1] -E Y A [ | = 0] . ] Only under specific assumptions can we conclude that ψ ( P ) is equivalent to ψ ∗ ( P ∗ ).

- 2) Suppose our parameter of interest is the mean survival time for a (sub)population. A simple measurable data frame is Z ∗ = ( X,T ), in which X denotes the covariates and T denotes the survival time (under specific survival probability), and therefore, the estimand is ψ ∗ ( P ∗ ) = E T [ ] 5 . However, we may encounter the problem that the survival time is censored, meaning that we cannot directly know the survival time T but know the censoring time T C . Therefore, the data structure for us is Z = ( X,T,δ,T C ), where δ signals if we occur censoring and T C denotes the censoring time. A possible (but uncommon) statistical estimand with the data structure is ψ ( P ) = E T δ [ | = 0], which ignores the censoring data. Further assumptions are required to make our statistical estimand approach the survival estimand. The detailed techniques for survival inference with truncation and censoring are the main part of Section 2 of this thesis.
- 3) Suppose we are interested in how a mediator interferes with the causal relationship between the treatment and the outcome, for instance, how much the causal effect directly goes from the treatment to the outcome (the direct effect) and how much it goes through the mediator onto the outcome (the indirect effect) (VanderWeele 2015, Chapter 2). We still assume that our treatment is a binary one A = 1 denotes the treatment, and A = 0 denotes the control. The mediator takes the value of m (1) when A = 1 and m (0) when A = 0. Thus, there are four combinations of

5 Of course, in real life analysis, we usually have a more complex data frame like Z = ( X,S t ( 1 ) , S ( t 2 ) , ..., S ( t n )), where S t ( i ) denotes the survival probability at t i . Thus, we first calculate the survival time under specific survival probability (for instance, half-life survival time, in which S t ( m ) = 0 5) with . S -1 , and average the survival time to get E T [ ].

the potential outcome: Y (1 , m (1)) , Y (1 , m (0)) , Y (0 , m (1)) , and Y (0 , m (0)). Our estimands are the two components of the total average treatment effect E Y [ (1 , m (1))] -E Y [ (0 , m (0))]: the direct effect, defined as the average effect of treatment in the absence of the mediator E Y [ (1 , m (0))] -E Y [ (0 , m (0))]; and the indirect effect, defined as the difference between the causal effects with and without the mediator E Y [ (1 , m (1))] -E Y [ (1 , m (0))]. Indeed, the most important parts, as shown in the decomposition, are the conditional response function : E Y [ ( a, m )] and we let the target mediation estimator ψ ∗ ( P ∗ ) to be that. From our statistical model, in which we have the data formed in a tuple: Z = ( X,A,M,Y ), we could yield the conditional expectation term, ψ ( P ) = E Y A [ | = a, M = m ]. Similarly, we need specific assumptions that allow us to equalize ψ ∗ ( P ∗ ) and ψ ( P ). The details of the technique are discussed in Section 3 of this thesis.

Since the content of causal inference (the first scenario above) goes throughout the thesis, we specifically analyze it in this chapter. Causal inference , broadly speaking, is a process that determines an independent effect of a particular object (the treatment) on another (the outcome), and it is usually contained in a larger system. For instance, when Galileo experimented on the Leaning Tower of Pisa, he isolated the independent effect of the weights of the two balls, observing that the heavier ball and the lighter ball fell on the ground simultaneously, and therefore inferred that the weights (masses) of the balls (the treatment) have nothing to do with the gravity (the outcome) of the two balls. The strategy to isolate the treatment effect is a controlled experiment , in which we are assured that the only difference between the two is the object we deem the treatment. However, in many scientific disciplines, we can not artificially manipulate and completely isolate the treatment 6 . Moreover, due to the requirement of repeatability in modern science, we usually observe the group-level differences between treated and untreated instead of the two individuals (balls). Based on this, experiment designers need techniques such as randomization or blind control for treatment assignment 7 . Experiments with randomized group assignment are commonly called Randomized Control Trials (RCT) .

Experiments are becoming increasingly common in social science studies, but most research still relies on observational data to infer causal relationships. The observational data does not have the RCT design: they do not randomize the samples to the treatment and control groups. They are not designed to isolate the effects of the treatment variable. Therefore, researchers need to use statistical

6 Indeed, in Galileo's experiment, the masses could not be the only difference between the two balls, as either the materials (densities) or the size must be different because the densities and the volumes determine the masses (in Galileo's original experiment he ensured the two balls were irony). Therefore, further experiment designs are needed: controlling the size of the two balls and controlling the type of materials of the two balls.

7 Randomization and blind-control might refer to different techniques. Consider the example of Pavlov's classic conditioning experiment with the dog, which is a blind instead of randomized control.

techniques to transfer the observational cases to approximate the RCT process, and therefore, causal inference with observational data is a pseudo-RCT . Our thesis mainly addresses statistical ideas on causal inference with observational data.

In this thesis, causal inference with observational data is the process using the observable estimator from the random samples ψ ( P n ) to estimate the causal estimand ψ ∗ ( P ∗ ) via the statistical estimand ψ ( P ). In empirical studies, the estimation functions for the estimand (population distribution) (estimandor) and for the observables (sample distribution) are always the same (for instance, the expectation ψ ( P ) = E P [ Z ] and ˆ ( ψ P n ) = E P n [ Z i ] have the same functional form), and the only difference between the two functions here is in the measurement choice. Thus, below, we will use ˆ and ψ ψ ( P n ) interchangeably to represent the estimator from the empirical dataset, we use ψ and ψ ( P ) interchangeably to represent the estimator (estimandor) from the statistical estimand.

Statisticians prefer ˆ as an unbiased estimator on ψ ψ ∗ and regard the divergence between the two terms as error terms. Moreover, we can decompose the error term into the statistical error, which is the divergence between the estimator and the statistical estimand, and the causal error, which is the divergence between the statistical estimand and the causal estimand:

$$\hat { \psi } - \psi ^ { * } = \underbrace { ( \hat { \psi } - \psi ) } _ { \text{statistical error} } + \underbrace { ( \psi - \psi ^ { * } ) } _ { \text{causal error} }$$

For the statistical error, we could further decompose it into the statistical variance and the statistical bias:

$$\hat { \psi } - \psi = \underbrace { ( \hat { \psi } - E [ \hat { \psi } ] ) } _ { \text{statistical variance} } + \underbrace { ( E [ \hat { \psi } ] - \psi ) } _ { \text{statistical bias} }$$

Analyzing the error terms is the core of the causal analysis, as in substantive work, especially from observational data, our estimator at most times could not be satisfied as unbiased yielding the causal estimand (and even in the RCT, further assumptions may be required). In Section II, we discuss the causal error, and in Section III, we discuss the statistical error.

## II. Analyzing the Causal Error

## A. Unbiased estimators and assumptions

The causal error is defined as the divergence between the statistical estimand ψ ( P ) (simply ψ ) and the causal estimand ψ ∗ ( P ∗ ) (simply ψ ∗ ). We cannot directly apply the statistical estimand as the causal one because of the fundamental problems of causal inference that we do not observe the outcome under different treatment conditions simultaneously, as each individual only receives one identifiable treatment. In Neyman-Rubin's (NR) causal framework (Neyman et al. 1935; Rubin 1990), the outcomes under different treatment statuses from the causal DGP ( Y ( a )) are the potential outcomes (Holland 1986). In the statistical data frame, the existing outcome is conditioned on the assigned value of the treatment ( Y A | = a ); nevertheless, we could not directly obtain the outcomes conditioned on the treatment assigned to other values. Therefore, the unobservable outcomes are called counterfactuals .

NR's counterfactual framework is not the only way to understand the causal and statistical estimand relationship. Computer scientist Judea Pearl (2009) created a framework called 'do-calculus' (DoC) . From Pearl's perspective, the outcome is deterministic if the treatment action has been triggered and the mapping rule from the treatment to the outcome is determined: if the treatment is A = a, , then a → Y ( a ); if A = a ′ , then a ′ → Y ( a ′ ). Based on the mapping symbol, Pearl developed the graphical expression for causal analysis called the Direct Acyclic Graphs (DAG) , and it is commonly used in structural equation models and causal mediation analysis.

Indeed, both NR's counterfactual and Pearl's DoC framework require additional assumptions to make ψ ∗ a ( P ∗ ) = E P ∗ [ Y ( a )] and ψ a ( P ) = E P [ Y A | = a ] equivalent. In this whole thesis, we consider the binary treatment. We suppose that the probability of being assigned as treatment and control is a positive number between 0 and 1:

ASSUMPTION II.1 (Positivity Assumption): 8 The probability of being assigned as treatment and control is a positive number between 0 and 1:

$$P ( A = 1 ) \in ( 0, 1 ) ; P ( A = 0 ) \in ( 0, 1 )$$

Where A is the treatment.

And suppose that in the 'omniscient' causal data frame Z ∗ = ( X,A,Y (1) , Y (0), we could observe the following conditional expectations: E Y [ (1) | A = 1], E Y [ (1) | A = 0], E Y [ (0) | A = 1], E Y [ (0) | A = 0], while in the statistical data frame Z = ( X,A,Y ) we could only observe E Y A [ | = 1] and E Y A [ | = 0]. With the following assumption, we could link between E Y A [ | = 1] and E Y [ (1) | A = 1], and between E Y A [ | = 0] and and E Y [ (0) | A = 0]:

ASSUMPTION II.2 (Consistent Assumption): The potential outcome under treatment received is the same as the observed outcome. That is,

$$Y = Y ( A )$$

8 In some literature, this assumption is also called overlap assumption (Heckman et al. 1998).

where Y is the observed outcome, Y ( A ) is the potential outcome, and A is the treatment.

Indeed, in the DoC framework, consistency has been implied since the mapping function from the execution of the treatment to the outcome is defined. Also, it is worth noting that when we use observational data to infer the causal estimand, consistent assumption needs to be held on the individual level: Y i = Y i ( A i ). Furthermore, since we have Assumption I.1 for the IID samples, we could infer that there's no interference among individuals: the treatment assigned to one observational sample does not affect the outcome of the others. Consistency and no interference assumptions on the individual observational level are collectively known as the Stable Treatment Unit Value Assumption , or SUTVA.

If we have established the relationship between the causal and statistical data frames, that E P [ Y A | = a ] = E P ∗ [ Y ( a ) | A = a ]. We suppose the proportion assigned to the treatment group is ρ , and therefore, E Y [ (1)] = ρE Y [ (1) | A = 1] + (1 -ρ E Y ) [ (1) | A = 0] and E Y [ (0)] = ρE Y [ (0) | A = 1] + (1 -ρ E Y ) [ (0) | A = 0] We may calculate the difference between the causal average treatment effect and the statistical treatment effect with a simple calculation:

$$& \underbrace { \left ( E [ Y ( 1 ) | A = 1 ] - E [ Y ( 0 ) | A = 0 ] \right ) } _ { \text{ statistical average treatment effect} } - \underbrace { \left ( E [ Y ( 1 ) ] - E [ Y ( 0 ) ] \right ) } _ { \text{ causal average treatment effect} } \\ & \quad = \left ( E [ Y ( 1 ) | A = 1 ] - E [ Y ( 0 ) | A = 0 ] \right ) \\ & \quad - \left [ \left ( \rho E [ Y ( 1 ) | A = 1 ] + ( 1 - \rho ) E [ Y ( 1 ) | A = 0 ] \right ) \\ & \quad - \left ( \rho E [ Y ( 0 ) | A = 1 ] + ( 1 - \rho ) E [ Y ( 0 ) | A = 0 ] \right ) \right ] \\ & \quad = \underbrace { E [ Y ( 0 ) | A = 1 ] - E [ Y ( 0 ) | A = 1 ] } _ { \text{ difference in baseline} } + \underbrace { \left ( 1 - \rho \right ) \left ( \delta _ { 1 } - \delta _ { 0 } \right ) } _ { \text{ heterogeneous treatment effect} } \\ \text{where}$$

where

$$\delta _ { 1 } = E [ Y ( 1 ) | A = 1 ] - E [ Y ( 0 ) | A = 1 ] \quad \text{and} \quad \delta _ { 0 } = E [ Y ( 1 ) | A = 0 ] - E [ Y ( 0 ) | A = 0 ].$$

## PROOF:

Let E Y [ (1) | A = 1] = α 1 , E Y [ (1) | A = 0] = α 2 , E Y [ (0) | A = 1] = α 3 , and E Y [ (0) | A = 0] = α 4 . Therefore, the left side of the equation is:

$$( \alpha _ { 1 } - \alpha _ { 4 } ) - \left [ \rho \alpha _ { 1 } + ( 1 - \rho ) \alpha _ { 2 } - \rho \alpha _ { 3 } - ( 1 - \rho ) \alpha _ { 4 } \right ].$$

Simplifying this, we have:

$$( \alpha _ { 1 } - \alpha _ { 4 } ) - [ \rho \alpha _ { 1 } + ( 1 - \rho ) \alpha _ { 2 } - \rho \alpha _ { 3 } - ( 1 - \rho ) \alpha _ { 4 } ]$$

$$& = ( \alpha _ { 1 } - \alpha _ { 4 } ) - \rho \alpha _ { 1 } - ( 1 - \rho ) \alpha _ { 2 } + \rho \alpha _ { 3 } + ( 1 - \rho ) \alpha _ { 4 } \\ & = \alpha _ { 1 } - \alpha _ { 4 } - \rho \alpha _ { 1 } - ( 1 - \rho ) \alpha _ { 2 } + \rho \alpha _ { 3 } + ( 1 - \rho ) \alpha _ { 4 } \\ & = \underbrace { ( \alpha _ { 3 } - \alpha _ { 4 } ) } _ { \text{difference in baseline} } + ( 1 - \rho ) \underbrace { ( \alpha _ { 1 } - \alpha _ { 3 } ) - ( \alpha _ { 2 } - \alpha _ { 4 } ) } _ { \text{heterogeneous treatment effect} } \.$$

Rewriting this back in terms of the original expectations, we have the right side of the equation.

Equation II.1 reveals the two origins of bias in causal inference when using the statistical estimator to infer the causal estimand: the baseline difference, or the selection bias, which is the pre-treatment divergence when grouping individuals to the treatment and control groups; and the heterogeneous treatment effect between the treatment and control group, which is the post-treatment divergence between the treatment and the control group. For example, consider evaluating the impact of a training program on workers' productivity. Initially, we measure their productivity levels before the training. Next, we divide the workers into two groups: a treatment group that receives the training and a control group that does not. After a certain period, we measure the change in productivity in both groups to assess the effect of the training program. The bias in this measurement comes from two sources: firstly, a pre-training bias, where workers in the treatment group might have different initial productivity levels compared to those in the control group; and secondly, a post-training bias, where workers in the treatment group might experience a greater improvement in productivity than those in the control group, even if they had all received the training.

Therefore, to eliminate the potential pre and post-treatment bias, we need further assumptions for identification. Since the pre-treatment selection and posttreatment heterogeneity can be attributed to the non-randomization in the treatment assignment, we have the ignorability/ unconfoundedness assumption:

ASSUMPTION II.3 (Ignorability/Unconfoundedness Assumption): The treatment assignment A is independent to the potential outcomes Y(1) and Y(0):

$$Y ( 1 ), Y ( 0 ) \perp A ^ { 9 }$$

Since the potential outcome is independent of the assignment of the treatment, we can infer that E Y [ (1) | A = 1] = E Y [ (1) | A = 0] and E Y [ (0) | A = 1] = E Y [ (0) | A = 0]. Therefore, under the unconfoundedness assumption, the pre-treatment baseline difference is 0. Meanwhile, the gap between E Y [ (1) | A = 1] and E Y [ (0) | A = 1] is the same as the gap between E Y [ (0) | A = 1] and E Y [ (0) | A = 0], eliminating the post-treatment heterogeneity 10 . In this sense, with Assumptions II.1,II.2, and

9 The symbol ⊥ ⊥ is the independent symbol, it means the two objects (vectors or matrices) are uncorrelated: a ⊥ ⊥ b ⇐⇒ cov a, b ( ) = 0 ⇐⇒ E a [ | b ] = 0. The three expressions are exchangeable in our use.

10 The elimination of post-treatment heterogeneity with the unconfoundedness assumption does not

II.3, we could finally conclude that the statistical average treatment effect is an unbiased estimand on the causal average treatment effect: ψ ( P ) = ψ ∗ ( P ∗ ).

In most circumstances, indeed, we may find a set of covariates X in the statistical model are correlated with both A and Y , violating the ignorability/unconfoundedness assumption. Thus, we may randomize the treatment assignment conditioned on X , which is a pseudo-randomization . We update Assumptions II.1 and II.3 to make them include the conditions of the covariates (the consistency hypothesis remains unchanged):

ASSUMPTION II.4 (Causal Inference Assumptions): Suppose a statistical DGP Z = ( X,A,Y ) , in which X denotes the covariates, A denotes the treatment, and Y denotes the outcome. To make the statistical estimand ψ ( P ) = E E [ X [ Y A | = 1 , X ]] -E E [ X [ Y A | = 1 , X ]] equivalent to the causal estimand ψ ∗ ( P ∗ ) = E Y [ (1)] -E Y [ (0)] from the causal DGP Z ∗ = ( X,A,Y (1) , Y (0)) (where Y (1) , Y (0) denote the potential outcomes under treatment and control, respectively), we need the following hypotheses:

- 1) Positivity: the probability to be assigned to treatment and control group conditioned on the covariates, is a positive number between 0 and 1:

$$P ( A = 1 | X ) \in ( 0, 1 ) ; P ( A = 0 |$$

- 2) Consistency: the potential outcome under the treatment received is the same as the observed outcome:

$$Y = Y ( A )$$

- 3) Unconfoundedness: conditional on a set of observed covariates X , the potential outcomes Y (1) and Y (0) are independent of the treatment assignment A :

$$\{ Y ( 1 ), Y ( 0 ) \} \perp A | X$$

Assumption II.4 is the sufficient and necessary condition for the statistical estimand on the average treatment effect to be equivalent to the causal estimand on the average treatment effect. To simplify, consider our target parameter is the potential outcome ψ ∗ a = E Y [ ( a )], and our statistical estimand is ψ a = E E [ X [ Y A | =

eliminate what econometrists called the heterogeneous treatment effect (HTE) in causal inference. First, although the average post-treatment heterogeneity between the treatment and control groups is eliminated, the individual treatment effect Y i (0) | A i = 1 and Y i (1) | A i = 0 still exists. Secondly, the HTE is indeed a conditional expectation E Y [ (1) | X ] -E Y [ (0) | X ] instead of what the unconfoundedness assumption controls E Y [ (1)] -E Y [ (0)].

a, X ]] (as E X [ Y ( a ) | X ] = ∫ yf Y ( a ) | X ( y X dy | ) ), we have:

$$\begin{array} {$$

Meanwhile, Equation II.2 can be also written as:

$$\text{Meanwhile, Equation II.2 can be also written as: } \\ \psi _ { a } ^ { a } & = E [ Y ( a ) ] \\ & = E [ E [ Y ( a ) | X ] ] \text{ (conditional expectation)} \\ & = E \left [ E [ Y ( a ) | X ] \frac { E [ 1 _ { A } | X ] } { E [ 1 _ { A } | X ] } \right ] \text{ (positivity; 1 _ { A } = 1 if A = a; 0 otherwise)} \\ & = E \left [ \frac { E [ Y ( a ) 1 _ { A } | X ] } { E [ 1 _ { A } | X ] } \right ] \text{ (unconfoundedness)} \\ & = E \left [ \frac { E [ Y 1 _ { A } | X ] } { P [ A = a | X ] } \right ] \text{ (consistency)} \\ & = E \left [ \frac { E [ Y 1 _ { A } | X ] } { \pi _ { a } ( X ) } \right ] \text{ (define $\pi_{a} ( X) = P ( A = a | X))} \\ & = E \left [ \frac { Y 1 _ { A } } { \pi _ { a } ( X ) } \right ] \text{ (reverse conditional expectation)} \\ \text{Equations II.2 and II.3 illustrate that we could infer the potential outcome $\psi_{a}$ or the propensity score function for} \\ \text{with either the conditional expectation $\psi_{a}$ or the propensity score function for} \\ \text{the treatment $\pi(a) as the unbiased estimators.  The results are the foundation}$$

Equations II.2 and II.3 illustrate that we could infer the potential outcome ψ ∗ a with either the conditional expectation ψ a or the propensity score function for the treatment π a ( ) as the unbiased estimators. The results are the foundation of what we refer to as the doubly robust/debiased estimation later in this chapter.

## B. Violations on causal assumptions

In most social science research scenarios with observational (survey) data, finding an unbiased causal estimator is challenging, as the three conditions in Assumption II.4 are not always satisfied, especially the positivity and the unconfoundedness assumptions.

## Violations on positivity

In social science, violating the positivity assumption is common if our treatment involves policy/reform/law enforcement that affects all our research objects. For instance, suppose our target is to measure how the rules on sports gambling may affect suicide risks for the residents of a state. Since the law affects everyone in the state, it violates the positivity assumption as P A ( ) = 1 (everyone in the state is grouped as treated). Violating the positivity assumption will make the casual estimand unidentifiable , as potential outcomes under counterfactual scenarios are nonexistent. Therefore, in the situations discussed above, further assumptions are required for causal inference.

## · Difference in Differences 11

A common technique to address the causal inference if the positivity assumption is violated is the difference in difference (DID) method (especially with longitudinal data), with additional assumptions and control group settings. For instance, in the example above, we may artificially choose the comparable control group as a neighboring state with similar socioeconomic factors but without changes to sports gambling laws. In this sense, we have implicitly assumed that the consistency and unconfoundedness assumptions to be true, since the treatment state and the control state will not affect the outcome of each other, and there are no other covariates affecting the change in suicide rates in both states during the observational time. Additionally, we need to assume that the suicide rate trend in the control state does not change before and after the change of the gambling law. In other words, it is called the parallel assumption (Roth et al. 2023):

ASSUMPTION II.5 ( Parallel Assumption for Difference in Differences): In the absence of the treatment, the difference in the average outcomes between the treated and control groups remains constant over time:

$$E [ Y _ { t _ { 1 } } ( 0 ) \ | \ A = 1 ] - E [ Y _ { t _ { 0 } } ( 0 ) \ | \ A = 1 ] = E [ Y _ { t _ { 1 } } ( 0 ) \ | \ A = 0 ] - E [ Y _ { t _ { 0 } } ( 0 ) \ | \ A = 0 ]$$

where Y t 1 and Y t 0 separately denotes the outcome before and after the treatment.

With Assumption II.5, we indeed transfer the causal effect of the policy change into the causal effect of time. With the introduction of a 'comparable' control group, the probability of receiving treatment becomes a real number between 0 and 1. Meanwhile, since E Y [ t 1 (0) | A = 1] is unobservable from the observational data, the assumption is not testifiable. The causal estimand, in this regard, should be the difference between the observable expected post-treatment outcome and the counterfactual post-treatment outcome (assuming the treatment was not received):

$$E [ Y _ { t _ { 1 } } ( 1 ) | A = 1 ] - E [ Y _ { t _ { 1 } } ( 0 ) | A = 1 ]$$

11 It is worth noting that methods mentioned here may be applied to scenarios which violate the other assumptions, or under which no assumption is violated.

Although the assumption is not testifiable, we can still justify the effectiveness of using the DID method. For instance, we can use statistical or visualization methods for the palacbo test , which examines whether the trends for the treatment group and the control group before the treatment were paralleled. Moreover, we could also apply a falsification test , which sets an outcome that is impossible to be affected by the treatment, to see if the trend for the outcome changes in the treatment group.

DID framework can be extended to the Difference in difference in differences (DDD) . In the DID framework, we assume that the unconfoundedness assumption is not violated, as we adjust for the time-invariant differences between groups and common trends affecting all units. Meanwhile, the unconfoundedness assumption may not hold, as some confounding variables might vary across groups and over time and further affect the parallel assumption. For instance, assuming we would like to evaluate the effectiveness of a specific educational policy on students' academic performance and choose a province for the experiment to further promote it to the whole nation, our treatment effect should be inferrable. Suppose the experimental province has only implemented it in its cities but not in the rural areas; if we only use the DID method to compare the effectiveness with the neighbor province while not noticing that we mixed the heterogeneous effects between urban and rural areas, our estimation may somehow underestimate the effectiveness. Therefore, we need to differentiate the DID results on the layer of urban or rural areas and approximate the causal effects more precisely.

## · Regression Discontinuity

Like DID, Regression Discontinuity (RD) is another method that, with additional assumptions of continuity at the breakpoint, addresses violations of positivity by assuming that the probability of treatment after the breakpoint is not always 1. For instance, suppose our target is to evaluate how tax cuts stimulate consumption for those with an annual income exceeding $ 100,000. Imagine the fiscal policy increases the tax rate for individuals earning beyond $ 100,000 from 25% to 30% (while remaining constant for people earning below $ 100,000 at 25%). Without further assumptions, we cannot identify the policy's effect because all individuals earning beyond $ 100,000 are assigned to the treatment group, violating the positivity assumption. As demonstrated in the DID section, we could assume the parallel trend to compare the treated and untreated groups (who earn below $ 100,000) before and after the treatment and identify the causal effect of the tax reform on consumption.

Moreover, we could have a continuity assumption under the RD setting, assuming the distribution is continuity at the policy's breakpoint. In the example above, we assume that there's a slight (no) difference between the consumption behaviors for those earning just below and just above $ 100,000 (how much the just is the radius of the threshold).

ASSUMPTION II.6 (Continuity Assumption for Regression Discontinuity): The expected potential outcomes E Y [ (0) | X = x ] and E Y [ (1) | X = x ] are continuous at the threshold c :

$$\lim _ { x \to c ^ { - } } E [ Y ( 0 ) | X = x ] = \lim _ { x \to c ^ { + } } E [ Y ( 0 ) | X = x ]$$

$$\lim _ { x \to c ^ { - } } E [ Y ( 1 ) | X = x ] = \lim _ { x \to c ^ { + } } E [ Y ( 1 ) | X = x ]$$

Where the potential outcomes are denoted as Y (1) if treated, and Y (0) if untreated. X i denotes the running variable.

Based on Assumption II.6, we made up the potential outcome where 0 &lt; lim x → c -P Y ( (1) | X = x ) &lt; 1 and 0 &lt; lim x → c + P Y ( (1) | X = x ) &lt; 1 to get over the positivity assumption. Therefore, due to Assumption II.4, we could assign the statistical estimand below the threshold as the causal estimand for the untreated group, whereas the statistical estimand above the threshold as the causal estimand for the treatment group:

$$\lim _ { x \to c ^ { - } } E [ Y | X = x ] = \lim _ { x \to c ^ { - } } E [ Y ( 0 ) | X = x ]$$

$$\lim _ { x \to c ^ { + } } E [ Y | X = x ] = \lim _ { x \to c ^ { + } } E [ Y ( 1 ) | X = x ]$$

and the difference in the conditional expectations of the observed outcomes on either side of the cutoff can be attributed to the treatment effect:

$$\lim _ { x \to c ^ { + } } E [ Y | X = x ] - \lim _ { x \to c ^ { - } } E [ Y | X = x ] = \lim _ { x \to c ^ { + } } E [ Y ( 1 ) | X = x ] - \lim _ { x \to c ^ { - } } E [ Y ( 0 ) | X = x ]$$

.

In the example above, suppose we take a radius of $ 5,000 and the threshold is $ 100,000. Therefore, we may calculate the average consumption costs for people after the tax reform earnings between $ 95,000 to $ 100,000 and for those earnings between $ 100,000 to $ 105,000. The treatment effect can be captured by the difference between the two expected consumption costs.

Applying the RD method has some restrictions on the running variable X and potential confounding variables. For the running variable, we have to assume that the patterns in the radius of the threshold are consistent 12 . In the above example, if the tax rates originally were different for people

12 In some literature, this is also called 'no manipulation' to ensure the continuity below and above the cutoff point.

earning below and above $ 100,000 (for instance, originally, people earning below $ 100,000 only had a tax rate of 15% and remained constant after the reform), the RD method requires further auxiliaries for identification. Or if the original tax rate cutoff point was at $ 97,500 (people earning below $ 97,500 received 15% tax rate while those earning above $ 97,500 received 25% tax rate), we need to redesign the bandwidth of the radius around the threshold. For the covariates that may interfere with the potential outcomes, the RD method also requires them to cross the threshold smoothly. In this way, the only difference between people under and above the threshold is whether they are assigned to the treatment group (above the threshold). In other words, the continuity assumption does not ensure the unconfoundedness condition is satisfied: as our estimand lim x → c + E Y [ (1) | X = x ] -lim x → c -E Y [ (0) | X = x ] is indeed the statistical estimand E Y [ (1) | A = 1] -E Y [ (0) | A = 0] in Equation II.1 and we still have the bias into baseline difference and heterogeneous treatment effect. For instance, in the previous example, we need to ensure that educational levels for people earning between $ 95,000 to $ 100,000 and $ 100,000 to $ 100,500 are almost the same, excluding the interference from the covariate of educational levels.

The positivity assumption requires that the probability assigned to the treatment and control group is a real number between 0 and 1, which, if violated, will make the causal estimand unidentifiable. DID and RD are, with additional assumptions, appropriate ways to construct the positivity in the probability of assigning cases to the treatment and control groups. Both methods are also helpful when addressing the violation of the unconfoundedness assumption- ignoring the possible interference from the unequal distribution of the covariates between the treatment and the control groups.

## Violations on unconfoundedness

Unlike the violation of the positivity assumption, which makes the causal effect unidentifiable, the violation of the unconfoundedness may only yield bias for estimating the causal effect. As Equation II.1 suggests, when the unconfoundedness assumption does not hold, the bias using statistical estimand to infer the causal estimand can be attributed to the pre-treatment selection bias and the post-treatment heterogeneity. The unconfoundedness assumption states that the outcome ( Y (1) , Y (0)) is independent of the treatment A given a set of covariates X .In econometrics terms, we call treatment A the endogenous treatment, while covariates X are exogenous covariates. Based on the relationship among X,Y, and A , previous researchers developed two intrinsically consistent ways to address the violation of the unconfoundedness assumption: the propensity score function-based matching and weighting method and eliminating exogeneity based on the local treatment average effect method, with instrumental variables and fixed effects.

## · Propensity Function based Matching and Weighting

With covariates affecting group assignment known, matching or weighting is the intuitive choice to eliminate the unconfoundedness bias. If covariates X are all known, as Equation II.3 indicates, with the intermediate estimator π a ( X ) = P A ( = a | X ), we can have the unbiased inference of the causal estimand E Y [ ( a )]. Traditionally, we call the intermediate estimator π a ( X ) as the propensity score function, as it measures the likelihood of treatment group assignment given the covariates X . The expression using the identity function divided by the propensity score function ✶ a π a ( X ) is called the inverse probability weighting (IPW) . Therefore, the IPW method yields the unbiased (and also regular and asymptotically linear, see Section III) estimation for causality if all assumptions hold.

The IPW method achieves pseudo-randomization since the reweighting process reassigns individuals into the pseudo-treatment groups based on their propensity (likelihood) rather than their true assignment status (which lacks the randomization process). The precondition for weighting is that covariates X affect the propensity and cause selection bias if not controlled. Therefore, if we have a longitudinal study, X should be variables ex-ante the treatment assignment. Suppose we have a set of variables M that affect the outcome ex-post of the treatment assignment. In that case, they will not affect our estimation of the causal effect A → Y . Still, they will decompose the total treatment effect into the direct one (the treatment directly affects the outcome) and the indirect (the causal effect goes via M onto Y ) effects, which we will discuss in the mediation section 13 .

Based on the propensity score function, weighting with inverse probability (propensity) provides a way to manipulate the treatment assignment to attain pseudo-randomization. Similarly, a matching method based on the propensity score function 14 can also achieve the effect of pseudo-randomization. This is the classic propensity score matching (PSM) method for causal inference (Rosenbaum and Rubin 1983). Suppose we could match the cases assigned in the treatment group and the control group with exactly the same propensity score and calculate the difference between the matched

̸

̸

13 The unconfoundedness assumption states ( Y (1) , Y (0)) ⊥ ⊥ A | X ⇐⇒ cov Y, A ( | X ) = 0, which does not indicate either cov Y, A ( ) = 0, cov X, A ( ) = 0, or cov Y, X ( ) = 0. Instead, cov Y, A ( ) = 0 and cov X, A ( ) = 0 are the foundations of causal inference (especially for Equation II.3. The relationship between Y and X is the trickiest: cov Y, X ( ) is not necessary 0, but cov Y, A ( | X ) ≡ 0 states that there might be some correlations between Y and X ), but such correlation has to be blocked under treatment A . Under this circumstance, the treatment assignment is randomized given X . In the next section, we will see that X is indeed the instrumental variable.

14 Propensity score function is not a score function that will be introduced in Section III. Therefore, to avoid confusion, in the following part of this chapter, we only call it 'propensity function.'

cases along the propensity score spectrum, and we finally calculate the average; we could yield the causal estimand for the average treatment effect based on the matching method.

The method seems plausible, but when dealing with real observational data, researchers have to make a tradeoff between the quality of the matching algorithm and the selection of cases. The real problem is that we rely on the observational data to generate the estimator ˆ π a to estimate π a . We can imagine that when using the observational data to estimate the propensity function, for the individuals in the treatment group, the distribution is likely to be dense at the end towards 1 (if 1 indicates being assigned to the treatment group) and relatively sparse at the end towards 0, while for the control group individuals tend to distribute denser on the side of 0 and more sparse on the side of 1. Therefore, it is infeasible to have a one-on-one match between the individuals from the treatment group and the control group, with the exact same propensity value, and get everyone matched (see Figure 1 for the illustration). Researchers have to adopt methods either to allow the divergence (caliper) in propensity scores between the matched cases, to drop the unmatched cases, or a method with the combination of the two (for instance, set a threshold for nearest neighborhood matching and drop the cases beyond the threshold).

The problem of PSM, as King and Nielsen (2019) advocate, is not on the process of causal inference (from π a to ψ ∗ a , as Equation II.3 suggests the unbiased process), nor on the process of statistical inference ˆ π k a to π a (where k denotes the k -th method for propensity estimation, and all ˆ π k a can be an unbiased estimator for π a if causal assumptions Assumption II.4 hold), it lies in the choice of ˆ π k a , or the problem they call 'model dependence': we rely on empirical observational data to simulate the DGP for propensity function, and further use the simulated function to predict the propensities for individuals from treatment and control groups (the dots in Figure 1), leaving cases unmatched due to uneven densities for the treatment and control groups. The bias, as King and Nielsen (2019) suggest, is a subjective bias originating from model choices: and the subjective choice of model somehow increases the imbalance 15 , model dependence, and bias for the causal estimation.

Intrinsically, the problem with PSM, if any, still results from the models to specify propensity function cannot satisfy the unconfoundedness assumption. Suppose the unconfoundedness assumption holds, which specifies the independence between the outcome and the treatment under the covariates.

15 Imbalance refers to the derivations to the exact match.

In that case, we can imagine that the expected outcomes for the treated and the control should have a constant distance along the propensity score (the slope does not necessarily have to be zero though). Thus, as long as we have a balanced match along the propensity score, the choices of ˆ π k a would be the same and unbiased estimator for π a (see the illustration of the middle and lower panel of Figure 1). However, since the remaining covariates still influence the distribution of outcomes for treatment and control groups based on the propensity scores (for example, in the lower panel of Figure 1, we have different slopes for treatment and control groups along the x-axis), and different matching criteria result in bias.

In summary, weighing and matching methods achieve causal inference by identifying the propensity function: π a ( X ) = P A [ = a X | ] and expect the estimated ˆ π a ( X ) is an unbiased estimator for π a ( X ). Unlike the local average treatment effect models, which restrict the exogeneity for the covariates (discussed below), we do not have any requirements for X in relation to Y . In other words, the propensity function only matters for X and A 16 , and it is almost unavoidable that we have statistical error between the estimated propensity score and the true value (similar to the omitted variable bias in regression analysis).

Thus, different models that researchers adopt will unavoidably generate different sample estimators. So, how do researchers claim that the causal effect they captured makes sense by adopting a specific weighting or matching method? We believe two things researchers need to claim before they describe their causal findings: one is the preconditions the models rely on: for instance, what covariates they have included and how they contribute to address or reduce the bias from confounding effects; the other is the theoretical guidance for them to choose the preconditions: the choice of the specific causal identification with variables included would be better theory-driven than pure data-driven.

Methods based on propensity functions with weighting and matching have multiple variant forms other than IPW and PSM. For instance, researchers could adopt an evolutionary search algorithm (Diamond and Sekhon 2013) or a Hungarian algorithm (Rosenbaum 1989) to optimize the matching process, or stratify the samples and match within different stratification (Rosenbaum and Rubin 1984), or match cases based on Mahalanobis distance (MDM, see Rubin 1980; King and Nielsen 2019).

16 Strictly speaking, the unconfoundedness assumption in the propensity function based method turns to: the outcome is independent to the treatment conditioned on the propensity score, ( Y (1) , Y (0)) ⊥ ⊥ A | π a ( X ).

## · Local Average Treatment Effect

The unconfoundedness assumption does not require the full model predicting the propensity for the treatment; indeed, we need the conditions 'isolating' the treatment and outcome relationship. We regard the causal relationship between the treatment and the outcome as the endogenous effect (in the isolated system), and the endogenous treatment effect is actually the local average treatment effect (LATE) (Imbens and Angrist 1994). The LATE captures the causal effect by nullifying the exogenous variation. According to Equation II.2, with no exogenous variability, the outcome conditioned on specific treatment status equals the potential outcome (i.e., E Y [ | A = a ] = Y ( a )). Classical instrumental variables (IV) method uses the exogenous instrumental variables, and fixed effect (FE) method focuses on the endogeneity effects within entities (groups, individuals, etc.). We discuss them briefly in this subsection. Meanwhile, as we suggest in the weighting and matching method, to specify the causal effect with the LATE method, we suggest researchers point out what the precondition, and in this case, the exogeneity condition, is for the inference.

## -Instrumental Variables

We start with the assumption of the IV method. A comprehensive review from the biostatistical and clinical perspective of the IV can be found in Baker et al. 2016. Specifically, we illustrate a binary IV scenario for simplifcation. In an empirical study, suppose we have assigned individuals to the treatment and control groups to measure the posttreatment divergence as the treatment effect. We have discussed that the main concern in this setting is the action of 'assignment'- whether it is 'random.' With random assignment, as we have discussed,

$$\hat { \psi } = E [ Y _ { i } | A _ { i } = 1 ] - E [ Y _ { i } | A _ { i } = 0 ] \\ = & \psi = E [ Y | A = 1 ] - E [ Y | A = 0 ] \\ = & \psi ^ { * } = E [ Y ( 1 ) ] - E [ Y ( 0 ) ],$$

as the fundamental of causal inference with statistical estimands.

Now that we consider some 'naughty' individuals do not want to follow the group assignment in the process described above, violating the randomization in the grouping. Let X denote the binary status of whether they are assigned to the treatment group, and A denote whether they receive treatment. We then have four observational groups: assigned treatment, received treatment: ( X i = 1 , A i = 1); assigned treatment, received control: ( X i = 1 , A i = 0); assigned control, received treatment: ( X i = 0 , A i = 1); and assigned control, received treatment:

( X i = 0 , A i = 0). The naughty individuals who choose whether to comply with the treatment assignment are either random or affect the outcome only through treatment 17 . X are the instruments, with the (instrumental) exogeneity assumption 18 :

ASSUMPTION II.7 (Instrumental Exogeneity Assumption): The instrumental variable selection needs to satisfy the following two conditions:

- a) Instrument relevance: he instrumental variable should be correlated with the treatment:

̸

$$c o v ( X, A ) \ne 0$$

- b) Instrument Exogeneity: the instrumental variable should be uncorrelated with the outcome fitted with the treatment:

$$c o v ( X, Y | A ) = 0$$

For each of the four observational groups, we may guess their motivations for complying with or violating the groups they are assigned to. For instance, individuals complying with the treatment group may be always-takers , as they intend to take the treatment no matter which group they are assigned, or compliers , who just comply with their assignment, while individuals violating the treatment group assignment may be defiers , who simply go to the opposite group they are assigned, or never-takers , who intend to go to the treatment group no matter which group they are assigned. Similarly, those who comply with the control group are either compliers or never-takers, and those who are assigned control but go to the treatment group are either defiers or always-takers. In other words, always-takers have individuals with ( X i = 1 , A i = 1) and ( X i = 0 , A i = 1), compliers are individuals with ( X i = 0 , A i = 0) and ( X i = 1 , A i = 1), never-takers are those ( X i = 0 , A i = 0) and ( X i = 1 , A i = 0), and defiers are those ( X i = 0 , A i = 1) and ( X i = 1 , A i = 0).

It is easy to notice that we can not identify the causal effects for alwaystakers and never-takers with the instrument, as the treatment effect is not affected by the instrument. Further, we need to assume that there are no defiers in the experiment, as we need the relationship between

17 This is not to say that the 'naughty decisions' are not affected by other covariates. For instance, some individuals care about their past medical history and decide not to take the treatment. Researchers do not care about the covariates: as long as there's no other reasonable correlation with the outcome except via the treatment.

18 The general exogeneity assumption suggests the explanatory variables should be uncorrelated with any error term: cov ϵ, X ( ) = 0, as ϵ = Y -f ( X ) if f ( X ) is the predicted value for the outcome given covariates are X .

the treatment and the instrument monotonic:

ASSUMPTION II.8 (Monotonic Assumption for Instrumental Variables): There are no defiers in the instrumental settings. In other words, the instruments affect the treatment in one direction (either always increases or always decreases the likelihood of receiving the treatment).

Thus, the treatment effect on the outcome for the compliers can be reliably estimated using the instrument, provided the exogeneity and monotonicity assumptions hold. Since adding the instrument meets all the assumptions necessary for causal inference, particularly unconfoundedness, the outcome is now independent of the treatment given the instrument. The Local Average Treatment Effect (LATE) for the compliers is determined by dividing the covariance between the outcome and the instrument by the covariance between the treatment and the instrument.

$$( \text{II.4} ) \quad \ \ L A T E = \frac { c o v ( X, Y ) } { c o v ( X, A ) } = \frac { E [ Y | X = 1 ] - E [ X = 0 ] } { E [ A | X = 1 ] - E [ A | X = 0 ] }$$

for the binary instrument. We call the generalized form (the division for the two covariances) as the IV estimand 19 . The pseudo-

19 IV method first appears in econometric literature not for identifying the causal effects or the endogeneity. It aimed to address the parameter identification in simultaneous equation models , in which the outcome is the equilibria of some structural relationships (we call these equilibria the structural form ). For instance, consider the outcome as an equilibrium of the price and the quality, say Y = f ( P, Q ). We have the equilibrium: Q s = Q d , where Q s and Q d represent the quality from the supply and the demand. Suppose we have the supply function (in linear form, with intercept excluded), and the demand function:

$$\begin{cases} P _ { s } = f ( Q _ { s }, X _ { s } ) + \epsilon _ { s } = \alpha _ { s } Q _ { s } + \beta _ { s } X _ { s } + \mu _ { s } \\ P _ { d } = f ( Q _ { d }, X _ { d } ) + \epsilon _ { d } = \alpha _ { d } Q _ { d } + \beta _ { d } X _ { d } + \mu _ { d } \end{cases}$$

where X ,X s d are separately the covariates affecting the relationship between price and quality for supply and demand (they are more likely different sets). With equilibrium conditions, Q and P are endogenous variables (since they appear in both structural forms), and we put the endogenous variables on the left side, while the exogenous variables on the right, so we have:

$$\begin{cases} P = \pi _ { 1 } X _ { s } + \pi _ { 2 } X _ { d } + \nu _ { p } \\ Q = \theta _ { 1 } X _ { s } + \theta _ { 2 } X _ { d } + \nu _ { q } \end{cases}$$

which are the reduced forms . With some algebra transformation, we could easily get:

$$\pi _ { 1 } = \frac { - \beta _ { d } } { \alpha _ { d } - \alpha _ { s } } ; \pi _ { 2 } = \frac { \beta _ { s } } { \alpha _ { d } - \alpha _ { s } } ; \theta _ { 1 } = \frac { - \alpha _ { s } \beta _ { d } } { \alpha _ { d } - \alpha _ { s } } ; \theta _ { 2 } = \frac { \alpha _ { d } \beta _ { s } } { \alpha _ { d } - \alpha _ { s } }$$

$$\alpha _ { s } = \frac { \theta _ { 1 } } { \pi _ { 1 } } ; \alpha _ { d } = \frac { \theta _ { 2 } } { \pi _ { 2 } } ; \beta _ { s } = \pi _ { 2 } ( \frac { \theta _ { 1 } } { \pi _ { 1 } } - \frac { \theta _ { 2 } } { \pi _ { 2 } } ) ; \beta _ { d } = - \pi _ { 1 } ( \frac { \theta _ { 1 } } { \pi _ { 1 } } - \frac { \theta _ { 2 } } { \pi _ { 2 } } ) \quad _ {. }$$

so, for example, if our target is the slope in the demand curve between Q and P : α d = θ 2 π 2 , we need the auxiliary variable X d and X s satisfying X d ⊥ ⊥ X s (otherwise, the estimation on θ 2 and π 2 are biased as we omitted the variables in the reduced form) a to establish models separately for Q and P , so that we get the unbiased estimator on the slope. The tricky thing is, in the simultaneous structural model for observational data, X d and X s will overlap; thus, in reality, we cannot differentiate θ 1 with θ 2 , and π 1

randomization is achieved with the assistance of the instrument for the compliers.

A classic example in social science is to measure the causal effect of military service on future earnings (Angrist 1990). Since the endogeneity of voluntary enlistments, individuals who choose to enlist may differ from those who do not, the causal relationship estimated directly from the model between military services and earnings yields bias. To address the endogeneity issue, Angist uses the draft lottery during the Vietnam War as the instrument: the draft lottery numbers are randomly assigned, while the eligibility to be drafted influences the likelihood of military service. Thus, those who are and are not induced to serve due to the draft lottery are the compliers. For the compliers, the study showed that those with military service experiences have lower future earnings than those without, so it verifies the negative causal effect of military service on future earnings.

Does this study have a good research design? First, we need to point out the external validity of the causality: in this study, Agrist focuses on the draft between 1970 and 1972 (so there are only three years for the draft) and men born between 1950 and 1953. So, rigorously speaking, the result can only be applied to the very specific men group in the 1950-53 cohort. Second, since many people resist being drafted, does draft evasion invalidate the results? Not at all! This is because the draft lottery is an exogenous variable, so if the results are restricted to with π 2 , so we will doubtlessly yield a biased estimator.

So the solution is to deploy an instrument and to circumvent direct estimation on the slope for the demand function; rather, we find the exogenous function which identifies Q s using the variable has nothing to do on the demand side, then rely on Q d = Q s to replace the quantity in demand with the quantity identified on the supply side, and finally calculate the slope between Q s and P s . For example, consider we are describing the demand curve between the price of grains and the quantity. We introduce the instrument, the rainfall, which affects only the quantity in supply instead of demand, and with the steps described, we can identify the demand curve slope (Pearl 2015; Haavelmo 1944; Stock and Trebbi 2003). With the exogenous variable uncorrelated with the error term in the output model, we address the consistent and unbiased correlation between the endogenous variables.

The process discussed above reveals the essence of the instrument variable lies in the unidentifiable problem in matrix operation. Let's further simplify the simultaneous structural model as:

$$\begin{cases} Y _ { n \times 1 } = X _ { n \times m } A 1 _ { m \times 1 } + Z _ { n \times t } B 1 _ { t \times 1 } + e 1 _ { n \times 1 } \\ X _ { n \times m } = Y _ { n \times 1 } A 2 _ { 1 \times m } + Z _ { n \times t } B 2 _ { t \times m } + e 2 _ { n \times m }. \end{cases}$$

In the simultaneous structural model, X and Y are endogenous variables, while Z is exogenous. The suffix denotes the dimension of the matrix (we suppose there are n cases, m endogenous covariates, and t exogenous covariates (assume t ≤ m ). To estimate the covariates, our first step is to bring the second equation directly to the right of the first equation, simplifying as Y n × 1 = Y n × 1 A 2 1 × m A 1 m × 1 + Z n × t ( B 2 t × m A 1 m × 1 + B 1 t × 1 )+( 2 e n × m A 1 m × 1 + 1 e n × 1 ). Let A = A A ,C 2 1 = ( B A 2 1+ B 1)( I -A ) -1 , We thus capture the reduced form: Y n × 1 = Z n × t C t × 1 + . Since e C is a t × 1 vector, we could not calculate specifically A ,A ,B ,B 1 2 1 2, unless we make some assumptions that some vectors in the matrices are 0. This is why we need the exogeneity assumption to presume further some uncorrelated relationship between the instrument and the error term to get the estimation in the endogenous relationship.

the compliers. But, this relies on the monotonic assumption, which is to say that those drafted but did not serve in the military are nevertakers: they determined that whatever the results are, they will not serve in the military- in other words, they are not defiers. Under this circumstance, the results are valid for the compliers- the average treatment effect yields between those drafted and did military service, and those not drafted and didn't go to the military.

In today's social science, 'finding an appropriate instrumental variable is rather an art than a science.' This is because having a strictly exogenous variable in observational survey data is hard and sometimes even requires imagination. Traditionally, social scientists rely on some naturally random assignments (Angrist and Krueger 2001), for instance, geographic or spatial variables (Card 1999) and weather variables (Dell et al. 2009), or deliberately engineered changes, for instance, policy changes (Angrist and Lavy 1999) and economic shocks (Autor et al. 2013), and even demographic (Oreopoulos 2006), biological and health (Fletcher and Wolfe 2009) 20 as the instrumental variables in their studies.

Finally, an unavoidable situation for social scientists is that the variables they have to rely on are weak instrumental variables (Stock et al. 2002) 21 . According to Assumption II.7, a valid instrument should satisfy the conditions of both relevance and exogeneity, a weak instrument is defined as its relevance with the treatment is weak. Therefore, the instrument does not explain much of the variation in the treatment, which is exogenous, leading to two consequences: first, we may not reliably estimate the LATE due to endogeneity, and our estimation will be biased towards the model which we didn't use any instrument (compared to a strong IV); second, with less precise estimates on the LATE, the confidence interval of our estimation will be wider. To illustrate this, consider the identification models between the outcome and the treatment, and between the treatment and the instrument:

$$\begin{cases} Y = g ( A ) + \epsilon \\ A = h ( X ) + \eta \end{cases}$$

In IV settings, the operation we estimate the LATE for the compliers is

20 When using demographic, biological, and health variables as the IV, the core research questions usually set on treatment effect of a policy on the socioeconomic outcome, while when using policy and economic shocks, the research questions are usually the relationship between socioeconomic variables.

21 The weak instrumental variable issue discusses more about the statistical error between the estimator from the observational data and the statistical estimand, which we will mainly focus on Section III. We decide to put it here to ensure coherence in the discussion of the IV method.

a two-step model (if we set the model in a parametric way, we use two ordinary least square models, and this is called the two-stage-leastsquare (2SLS) ). We first estimate ˆ = ˆ ( A h X ), and put the predicted ˆ A in the first equation to yield the IV estimator ˆ g IV ( A ) = E Y [ | ˆ ]. A If the IV is correctly identified, we have ˆ A → h X ( ), meanwhile E ϵ [ | X = 0]. Therefore,

(II.5)

$$\dot { \hat { g } } _ { I V } ^ { ^ { * } } ( A ) = E [ Y \ | \ \hat { A } ] \to E [ Y \ | \ h ( X ) ] = E [ g ( A ) + \epsilon \ | \ h ( X ) ] = g ( A )$$

According to the definition of weak IV, which has a weak relevance with the treatment, we may assume that ˆ = ˆ ( A h X ) → 0. Still, E ϵ X [ | = 0]. Therefore,

$$( \text{II.6} ) & = E [ Y \ | \ \hat { A } \to 0 ] = E [ Y \ | \ \hat { A } \to 0 ] = E [ g ( A ) + \epsilon \ | \ \hat { A } \to 0 ] \\ & = E [ g ( A ) \ | \ \hat { A } \to 0 ] + \underbrace { E [ \epsilon \ | \ \hat { A } \to 0 ] } _ { \colon = 0 } = E [ g ( E [ A ] ) \ | \ \hat { A } \to 0 ] = g ( E [ A ] )$$

Therefore, the estimator from the weak IV approaches to the result of purely the second-stage model in the two-stage estimation if the relevance between the instrument and the treatment approaches zero. Combine Equations II.5 and II.6, we have the bias from weak IV, which is:

$$\text{Bias} = g ( E [ A ] ) - g ( A ).$$

## -Fixed Effects

As discussed above, the essence of the LATE method is to isolate the exogenous variation in the endogenous variables (the treatment) and obtain an unbiased and consistent estimator for the treatment effect. The fixed effect method could also be used to make the unconfoundedness assumption hold by isolating the exogenous variation, as it focuses on the variations within the entity. We assume the values of the outcome variable follow a normal distribution:

$$y _ { i j } \sim N ( \alpha _ { j } + f ( a _ { i j } ), \sigma _ { y } ^ { 2 } )$$

Suffix i refers to the individuals, or the granular level of data, while suffix j refers to the group unit, which aggregates the individual-level data. alpha j denotes the group-level specific term, which accounts for the variability in the mean response y ij between different groups j f . ( a ij ) represents the effect of the treatment a ij on the response variable y ij , and σ 2 y denotes the variance for the outcome. We could further rewrite the group-level term as 22 :

$$\alpha _ { j } \sim N ( \mu, \sigma _ { \alpha } ^ { 2 } )$$

̸

In which µ denotes the overall mean effect across all groups, and σ 2 α represents the variance of the group-specific effects 23 . Since we set up the distribution for a normal, we have specified the unrelated relationship between the treatment effect f ( a ij ), the group-level effect α j , and the statistical error (variance) of the outcome ϵ ij for any given group j ( i.e., E σ [ y | f ( a i 1 ) , f ( a i 2 ) , ..., f ( a iJ ) , α j ] = 0). Also, we assume that the error terms between any two groups are uncorrelated: E σ [ y it | σ y is ] = 0 if t = . s

With the assumptions hold, f ( a ij ) in the model denotes the LATE for individuals. For instance, consider our aim is to estimate an individual's education a ij on her earnings y ij from the panel data. Therefore, the group j = 1 2 , , ..., J represents the specific waves, and we could simply calculate the ATE of education on earnings simply by demeaning the outcome:

$$( y _ { i j } - \bar { y } _ { i } ) = ( f ( a _ { i j } ) - f ( \bar { a } _ { i } ) ) + ( \sigma _ { y } - \sigma _ { \alpha } )$$

In the expression above, ¯ y i and ¯ a i separately represent the mean outcome and treatment for the specific individual. Since the exogenous variables (i.e., time-invariant individual characteristics, ability and intelligence, and initial conditions, like early childhood health) are controlled by the fixed term, we could make inferences on our treatment and the outcome, supposing there are no time-varying variables affecting the results (which is restricted by E σ [ y it | σ y is ] = 0 in our model setting).

As can be seen from the above example, if we set the group level as time, the fixed effect is indeed the same as the difference-in-difference model that we discussed before. In other words, DID is a specific application of the FE under the panel data settings. We discussed the advantage of DID in addressing the violation of the positiveness assumption, and indeed, it addresses the violation of the unconfoundedness assumption, simply by ignoring the exogeneity fixing the timeinvariant effects.

22 The assigned probability distribution for α j in multilevel models is another 'soft constraint' in this type of models (see Gelman and Hill 2006, ch. 12, p. 257).

23 the representation is equivalent to the linear form expression: y j i = µ + ( f a j i ) + ϵ α j + ϵ y ij , in which we assume the error term ϵ y ij is uncorrelated with the treatment a j i , group-level fixed mean µ , and the group-level error term ϵ α j .

## III. Analyzing the Statistical Error

In the section above, we mainly discussed the methods addressing the causal error, in other words, how we use the statistical estimand to approximate the causal estimand ψ -ψ ∗ . We mentioned in passing some content of this section, which is how to estimate from the statistical estimator from the observational data to approximate the statistical estimand: the statistical error ˆ ψ -ψ . According to Equation I.3, the statistical error can be decomposed in the statistical variance ˆ ψ -E ψ [ ˆ ] and the statistical bias E ψ [ ˆ -ψ ]. In (especially) machine learning literature, we often discuss a balance to be struck between minimizing the bias or minimizing the variance to achieve optimal model performance, and the balance is called variance-bias tradeoff 24 . Although some biased estimators might be useful as they could optimize the model with a minimal mean square error (MSE) 25 , in causal inference methods discussed in the whole thesis, when

24 We will simply find we could express the mean square error (MSE) of the statistical error by the statistical variance and the statistical bias: the Mean Squared Error (MSE) of an estimator ˆ is defined ψ as: MSE( ˆ ) = ψ E ψ [( ˆ -ψ ) 2 ]. We add and subtract the expected value of the estimator E ψ [ ˆ ] inside the squared term:

$$M S E ( \hat { \psi } ) = E [ ( \hat { \psi } - E [ \hat { \psi } ] + E [ \hat { \psi } ] - \psi ) ^ { 2 } ]$$

Then expand the squared term:

$$M S E ( \hat { \psi } ) = E [ ( \hat { \psi } - E [ \hat { \psi } ] ) ^ { 2 } + 2 ( \hat { \psi } - E [ \hat { \psi } ] ) ( E [ \hat { \psi } ] - \psi ) + ( E [ \hat { \psi } ] - \psi ) ^ { 2 } ]$$

Simplify by recognizing that E ψ [ ˆ -E ψ [ ˆ ]] = 0 (since the expectation of a deviation from the mean is zero):

$$\text{MSE} ( \hat { \psi } ) = E [ ( \hat { \psi } - E [ \hat { \psi } ] ) ^ { 2 } ] + 2 E [ \hat { \psi } - E [ \hat { \psi } ] ] ( E [ \hat { \psi } ] - \psi ) + ( E [ \hat { \psi } ] - \psi ) ^ { 2 } \\ \text{MSE} ( \hat { \psi } ) = E [ ( \hat { \psi } - E [ \hat { \psi } ] ) ^ { 2 } ] + ( E [ \hat { \psi } ] - \psi ) ^ { 2 }$$

Recognize that E ψ [ ˆ -E ψ [ ˆ ]] = 0 means the middle term drops out:

$$M S E ( \hat { \psi } ) = V a r ( \hat { \psi } ) + B i a s ( \hat { \psi } ) ^ { 2 }$$

Thus, the decomposition of the Mean Squared Error (MSE) is:

$$M E S ( \hat { \psi } ) = \text{Var} ( \hat { \psi } ) + \text{Bias} ( \hat { \psi } ) ^ { 2 }$$

.

25 A common case for a biased estimator is to estimate the variance of a normal distribution, which yields smaller MSE than use the unbiased estimator. Given a set of independent and identically distributed (i.i.d.) samples X ,X ,...,X 1 2 n from a normal distribution N µ,σ ( 2 ), we can estimate the population variance σ 2 . The sample variance S 2 is an unbiased estimator of the population variance σ 2 . It is defined as:

$$S ^ { 2 } = \frac { 1 } { n - 1 } \sum _ { i = 1 } ^ { n } ( X _ { i } - \bar { X } ) ^ { 2 },$$

where ¯ X is the sample mean: ¯ = X 1 n ∑ n i =1 X . i . The sample variance S 2 satisfies: E S [ 2 ] = σ . 2 . A biased estimator of the population variance is the biased sample variance, defined as:

$$\sigma _ { \text{biased} } ^ { 2 } = \frac { 1 } { n } \sum _ { i = 1 } ^ { n } ( X _ { i } - \bar { X } ) ^ { 2 }.$$

This estimator is biased because: E σ [ 2 biased ] = n -1 n σ . 2 Despite its bias, σ 2 biased can be preferred (for instance, in industrial manufacturing). This is because, in some cases, the latter is easy to compute, while we do not need a precise estimation of the overall variance. Plus, the reduction in variance provided we calculate the statistical estimator, we require it to be an unbiased estimator on the statistical estimand, so we have E ψ [ ˆ ] -ψ ≡ 0.

The analytical tool on how our estimator can infer the estimand unbiased is the asymptotical analysis . We imagine our estimator ˆ = ψ ψ ( P n ) is a nuisance perturbed estimation on the true estimand ψ ( P . Thus, we denote the perturbed distribution as ψ ( ˜ P ϵ ) ∣ ∣ ϵ =0 , where ϵ denotes the infinitesimal error term which approaches zero. Our core questions are two folds: 1) to which degree ˜ P ϵ drifted from the true measure P ; and 2) how the perturbation in the distribution of ˜ P ϵ affects the estimation of ˆ on ψ ψ . We use two functional forms to identify the two questions: we use the score function to determine the direction of the drift in measure, and we use the influence function to identify the effect of the perturbation. As we will elaborate on this section, from the perspective of asymptotical analysis, the best causal estimator we will derive is the unbiased efficient estimator under the condition of regularity and asymptotical linearity: this estimator is not only unbiased to the statistical estimand, but also yields the smallest variance. We will show below that the efficient estimator is also robust (Neyman orthogonal) from the econometric perspective, and the influence function is also efficient.

Due to the unfamiliarity of the materials discussed in this section for readers, we arrange it into the following parts: we will first introduce the characteristics of the regular and asymptotically linear estimators, and we further discuss the detailed characteristics of the score function and the influence function, revealing their connections and how they assist us in finding the efficient estimator. Finally, we derive the efficient influence function and the efficient estimator for the average treatment effect. In the following chapters of this thesis, we will still use the mathematical equations in this section to derive different efficient estimators adapting to various social science and demographic research scenarios.

## A. Introduction to RAL Estimators

For all the RAL estimators, it would be best to find the one with the lowest variance so that the MSE for the estimator towards the statistical estimand would be lowest. Fortunately, for the unbiased estimator, we know the lowest bound for the variance, given by the Cramer-Rao bound :

LEMMA III.1 (Cramer-Rao Bound): The Cramer-Rao Bound (CRB) states that for an unbiased estimator ˆ = ψ ψ ( P n ) of ψ = ψ ( P ) , the variance of ˆ ψ is at least as large as the inverse of the Fisher information:

$$V a r ( \hat { \psi } ) \geq \frac { 1 } { I ( \psi ) }$$

by the biased estimator can outweigh the increase in bias (especially in smaller sample cases), leading to a lower overall MSE.

where I ( ψ ) is the Fisher information given by:

$$I ( \psi ) = \mathbb { E } \left [ \left ( \frac { \partial } { \partial \psi } \log f ( Z ; \psi ) \right ) ^ { 2 } \right ]$$

The proof will be given later in this section. For all the unbiased estimators, if the variance attains the CRB, we call the estimators the efficient estimators . However, as we will elaborate, efficient estimators are sometimes hard to derive due to the complexity of the underlying functional form (mostly the influence function; see below), and we need deterministic or heuristic approaches that will make our estimator more efficient than any other peers we may find.

Besides unbiased, In our study, due to the central limit theorem (Theorem I.1), we further require the estimator to be regular and asymptotically linear (RAL) Regularity . in general means that our estimator processes towards some desirable characteristics as the sample size increases:

$$\sqrt { n } \left ( \psi ( \mathbb { P } _ { n } ) - \psi ( \mathcal { P } ) \right ) \stackrel { \mathcal { P } } { \sim } \mathcal { D },$$

Where P denotes an empirical measure which we will discuss shortly, and D denotes a fixed distribution 26 . Asympotitical linearity suggests the difference between the estimator and the statistical estimand can be approximated by the linear combination of the appropriate normalized sum of the IID random variables. Suppose the function is ϕ ψ ( ; P ; Z i ), for the estimator ψ , measure P and dataset Z i = ( X , Y i i ) 27 ,therefore,

$$\sqrt { n } \left ( ( \psi ( \mathbb { P } _ { n } ) - \psi ( \mathcal { P } ) ) - \frac { 1 } { n } \sum _ { i = 1 } ^ { n } \phi ( \psi ; \mathbb { P } _ { n } ; Z _ { i } ) \right ) \xrightarrow { p r o b. } 0,$$

We call ϕ ψ ( ; P n ; Z i ) as the influence function for the empirical data Z i with

26 In this thesis, the denotation on the convergence uses the expressions in van der Vaart (1998). In short, convergence in distribution (weak convergence) X n D - → X ⇐⇒ X n ⇝ X suggests that a sequence of random variables { X n } converges in distribution to a random variable X if for all points t at which F X ( ) is continuous:lim t n →∞ F X n ( ) = t F X ( ) t , where F X n ( ) and t F X ( ) are the cumulative distribution t functions of X n and X ; convergence in probability X n P - → X (or X n prob. - - - -→ X , to differentiate with the measure P ) suggests the relationship between X n and X , for every ϵ &gt; 0, is: lim n →∞ P ( | X n -X | ≥ ϵ ) = 0. Almost sure convergence X n a.s. - - -→ X suggests X n almost surely converges to XP (lim n →∞ X n = X ) = 1. Meanwhile, we also use big-O probability and small-o probability to denote convergence: big O p : X n = O p ( a n ) indicates that the sequence of random variable X n is bounded by a n in probability: for every ϵ &gt; 0, there exist constants M &gt; 0 and N &gt; 0 such that P ( | X n | ≤ Ma n ) ≥ 1 -ϵ for all n ≥ N. Small o p suggests that X n is asymptotically smaller than a n as n increases, or the difference between X n and a n is negligible: for every ϵ &gt; 0 and δ &gt; 0, P ( | X n | ≥ δa n ) → 0 as n →∞ .

27 The influence function has three entries, ϕ refers to the estimation functional form, P refers to the measure. Z i refers to the measurable set (dataset). Due to our assumptions in Section I that the estimation function and the measure (DGP) will not change simultaneously, thus, if we specify the influence function for an estimator, we omit the measure and the dataset. We use Z i for discrete elements in the measurable set and z if the elements are continuous.

respect to the empirical measure P n . With Theorem I.1, asymptotical linearity could also be expressed as:

$$\sqrt { n } ( \psi ( \mathbb { P } _ { n } ) - \psi ( \mathcal { P } ) ) \stackrel { d } { \to } \mathcal { N } ( 0, \sigma ^ { 2 } ).$$

where σ 2 = Var( ϕ ψ ( ; P n ; Z i )) and we call it the asymptotic variance . Equation III.2 can also be expressed as:

$$\left ( \psi ( \mathbb { P } _ { n } ) - \psi ( \mathcal { P } ) \right ) - \frac { 1 } { n } \sum _ { i = 1 } ^ { n } \phi ( \psi ; \mathbb { P } _ { n } ; Z _ { i } ) = o _ { p } ( n ^ { - 1 / 2 } )$$

which stresses the approximation error (remainder term) becomes negligible faster than n -1 2 / in probability. The property revealed in Equation III.4 is very useful when we try to yield the doubly robust/debiased machine learning estimator (in some literature, we also denote R 2 ( P n ; P ) = ( ( ψ ( P n ) -ψ ( P )) -1 n ∑ n i =1 ϕ ψ ( ; P n ; Z i ) ) and call it as the second-order remainder ) 28 . We will discuss regularity and asymptotically linearity in detail in the following subsections.

## B. Regularity and Score Function

We start with the definition of the score function . Suppose we have two ( σ -finite) measures P ( Z ) and ˜ ( P Z ) where P ( Z ) dominates ˜ ( P Z ) (or ˜ ( P Z ) is absolutely continuous with respect to P ( Z )). Now imagine we have a differentiable path starting from P ( Z ) and ending at ˜ ( P Z ). The most convenient method to define it is linear interpolation: let ˜ ( P ϵ Z ) = ϵ ˜ ( P Z ) + (1 -ϵ ) P ( Z )( ϵ ∈ [0 , 1]) . Therefore, suppose the probability density function for P ( Z ) and ˜ ( P Z ) are respectively p z ( ) and ˜( p z ), then the density function for ˜ ( p ϵ z ) = ϵp z ˜( ) + (1 -ϵ p z ) ( ) . PROOF:

We have the definition of a probability density function p z ( ) for a measure P such that for any measurable set Z :

$$\mathcal { P } ( Z ) = \int _ { Z } p ( z ) d z$$

Thus,

$$\tilde { \mathcal { P } } _ { \epsilon } ( Z ) = \epsilon \int _ { Z } \tilde { p } ( z ) \, d z + ( 1 - \epsilon ) \int _ { Z } p ( z ) \, d z = \int _ { Z } \left ( \epsilon \tilde { p } ( z ) + ( 1 - \epsilon ) p ( z ) \right ) d z$$

And therefore,

$$\tilde { p } _ { \epsilon } ( z ) = \epsilon \tilde { p } ( z ) + ( 1 - \epsilon ) p ( z ).$$

28 Not all unbiased estimators are RAL. For instance, consider the median estimator for all symmetric distributions. It is an unbiased estimator for the population median, but it is not asymptotically linear as ϕ ψ Z ( ; ) for the median is not a smooth function. Another classic example here is the Hodges' estimator.

With the differentiable path, we define the score function corresponding to the path from P ( Z ) towards ˜ ( P Z ) as the rate ϵ 0 change of the log-likelihood at the starting point P ( Z ) (the gradient of the log-likelihood function):

$$s _ { \epsilon _ { 0 } } ( z ) & = \frac { \partial \log \tilde { p } _ { \epsilon } ( z ) } { \partial \epsilon } \Big | _ { \epsilon = \epsilon _ { 0 } }$$

Usually, in our discussion, we set ϵ 0 = 0. According to Equation III.5 29 , we may rewrite ˜ ( P ϵ Z ) and ˜ ( p ϵ z ) as: ˜ ( P ϵ Z ) = ∫ Z (1+ ϵs z ( )) d P ( Z ) and ˜ ( p ϵ z ) = (1+ ϵs z ( )) p . PROOF:

We first prove that the definition of s z ( ) = ∂logp ˜ ( ϵ z ) ∂ϵ ∣ ∣ ∣ ∣ ϵ =0 can be rewrite as:

$$s ( z ) = \frac { \tilde { p } ( z ) } { p ( z ) } - 1$$

Using the chain rule to differentiate:

$$\frac { \partial \log \tilde { p } _ { \epsilon } ( z ) } { \partial \epsilon } \Big | _ { \epsilon = 0 } & = \frac { \partial \log \tilde { p } _ { \epsilon } ( z ) } { \partial \tilde { p } _ { \epsilon } ( z ) } \cdot \frac { \partial \tilde { p } _ { \epsilon } ( z ) } { \partial \epsilon } \Big | _ { \epsilon = 0 } = \frac { 1 } { \tilde { p } _ { \epsilon } ( z ) } \cdot \frac { \partial } { \partial \epsilon } ( \epsilon \tilde { p } ( z ) + ( 1 - \epsilon ) p ( z ) ) \Big | _ { \epsilon = 0 } \\ & = \frac { 1 } { p ( z ) } ( \tilde { p } ( z ) - p ( z ) ) = \frac { \tilde { p } ( z ) } { p ( z ) } - 1$$

Therefore, ˜( p z ) = ( s z ( ) + 1) p z ( ). We could rewrite ˜ ( p ϵ z ) as ˜ ( p ϵ z ) = (1 -ϵ p z ) ( ) + ϵ s z ( ( ) + 1) p z ( ). Simplifying the expression:

$$\tilde { p } _ { \epsilon } ( z ) = \left ( 1 - \epsilon + \epsilon s ( z ) + \epsilon \right ) p ( z ) = \left ( 1 + \epsilon s ( z ) \right ) p ( z )$$

Further, to prove that ˜ ( P ϵ Z ) = ∫ Z (1+ ϵs z ( )) d P ( Z ), we need the Radon-Nikodym theorem :

THEOREM III.1 (Radon-Nikodym Theorem): Let (Ω , F P , ) be a probability space, and let ˜ P ϵ be another probability measure on (Ω , F ) such that measure ˜ P ϵ is absolutely continuous with respect to measure P . Then there exists a P -integrable function f : Ω → [0 , ∞ ) such that for every A ∈ F ,

$$\tilde { \mathcal { P } } _ { \epsilon } ( A ) = \int _ { A } f \, d \mathcal { P }.$$

The function f is called the Radon-Nikodym derivative and is often denoted by d ˜ P ϵ d P .

The theorem gives us a toolbox to 'rescale' the measure ˜ P ϵ with the measure P and the Radon-Nikodym derivative function f . Since Equation III.7, we

29 Generally, if ϵ 0 = 0, we simplify the score function s ϵ ( z ) as s z ( ), or further, s .

the scale is indeed (1 + ϵs z ( ), and thus,

$$\left ( \text{III.8} \right ) \quad \ \tilde { \mathcal { P } } _ { \epsilon } ( Z ) = \int _ { Z } \left ( 1 + \epsilon s ( z ) \right ) p ( z ) d z = \int _ { Z } \left ( 1 + \epsilon s ( z ) \right ) d \mathcal { P } ( Z ).$$

The definition and Equation III.7 indeed reveal that the essence of the score function is a 'direction pointer' (or compass). We could use the score function to specifically point the direction of the digress between the empirical estimator and the true statistical estimand. Suppose we have the estimand as ψ ( P ), the pathway derivative, or the gradient of the estimand in the direction of the score function 30 , can be defined as:

$$\lim _ { \epsilon \to 0 } \frac { \psi ( \tilde { \mathcal { P } } _ { \epsilon } ) - \psi ( \mathcal { P } ) } { \epsilon } = \nabla _ { s } \psi ( \mathcal { P } )$$

Where we could regard ψ ( ˜ P ϵ ) as a perturbed version of the true estimand ψ ( P ), adjusted by a small amount of ϵ in the direction of some perturbation.

An obvious and important characteristic of the score function s z ( ) is that E s z [ ( )] = 0 (so that ˜ ( P ϵ Z ) = ∫ Z (1 + ϵs z ( )) d P ( Z ) can be integrated to 1) and Var( ( s z )) &lt; ∞ 31 . Given this, we have an important corollary. Suppose we have a bivariate joint distribution P ( Y, X ), clearly, we have the Bayesian rule P ( Y, X ) = P ( Y | X ) P ( X ). We define the score function for P ( Y | X ) and P ( X ) separately as s Y X | ( x, y ) and s X ( x ). The score functions satisfy E s [ Y X | ( x, y ) | X ] = 0 and E s [ X ( x )] = 0. Therefore, s Y X | ( x, y ) ⊥ s X ( x ). PROOF:

The target is to prove the expectation E s [ Y X | ( x, y ) s X ( x )] = 0 . Rewrite the expectation with the law of conditional expectation:

$$E [ s _ { Y | X } ( x, y ) s _ { X } ( x ) ] & = E [ E [ s _ { Y | X } ( x, y ) s _ { X } ( x ) \ | \ X ] ] \\ & = E [ s _ { X } ( x ) E [ s _ { Y | X } ( x, y ) \ | \ X ] ] = E [ s _ { X } ( x ) \cdot 0 ] = 0.$$

With the orthogonal relationship between s Y X | ( x, y ) and s X ( x ), we could have s X,Y ( x, y ) = s Y X | ( x, y ) + s X ( x ) and ∇ s X,Y ψ = ∇ s Y X | ψ + ∇ s X ψ . Therefore, we could calculate the score function for the marginal distribution P ( X ) and the conditional distribution P ( Y X | ) within their own models and sum them together to get the score function for the joint distribution if the score function for the joint distribution is hard to capture.

30

The gradient here is defined on the direction of the score function for the estimator. Indeed, the score function is also a gradient, but the gradient on the measurable set Z : s z ( ) = ∂ log ˜ p ϵ ( z ) ∂ϵ ∣ ∣ ϵ =0 = ∇ z log p z ( ). 31 In other words, score function s z ( ) is defined on the L 0 2 space as its square-integral functions are integrable within finite values and the mean is zero. Besides the score function, the influence function, the efficient estimation, and the Cramer-Rao bound are all defined on the L 0 2 space in our discussion.

## PROOF:

The joint distribution P ( Y, X ) can be expressed using the chain rule of probability:

$$p _ { X, Y } ( x, y ) = p _ { Y | X } ( y | x ) p _ { X } ( x ) \iff \log p _ { X, Y } ( x, y ) = \log p _ { Y | X } ( y | x ) + \log p _ { X } ( x )$$

Applying the definition of the score function, we differentiate both sides with respect to ϵ and evaluate at ϵ = 0:

$$s _ { X, Y } ( x, y ) = \frac { \partial \log \tilde { p } _ { X, Y, \epsilon } ( x, y ) } { \partial \epsilon } \Big | _ { \epsilon = 0 } = \frac { \partial \log \tilde { p } _ { Y | X, \epsilon } ( y | x ) } { \partial \epsilon } \Big | _ { \epsilon = 0 } + \frac { \partial \log \tilde { p } _ { X, \epsilon } ( x ) } { \partial \epsilon } \Big | _ { \epsilon = 0 }$$

Thus:

$$( \text{III.10} ) & & s _ { X, Y } ( x, y ) = s _ { Y | X } ( x, y ) + s _ { X } ( x )$$

By the definition of the gradient in the direction of the score function, we have:

$$\nabla _ { s _ { X, Y } } \psi = \lim _ { \epsilon \rightarrow 0 } \frac { \psi ( \tilde { \mathcal { P } } _ { \epsilon } ( X, Y ) ) - \psi ( \mathcal { P } ( X, Y ) ) } { \epsilon }$$

Since p X,Y ( x, y ) = p Y X | ( x, y ) p X ( x ), we have:

$$\nabla _ { s _ { X, Y } } \psi = \lim _ { \epsilon \rightarrow 0 } \frac { \psi ( \tilde { \mathcal { P } } _ { \epsilon } ( X, Y ) \tilde { \mathcal { P } } _ { \epsilon } ( X ) ) - \psi ( \mathcal { P } ( Y | X ) \mathcal { P } ( X ) ) } { \epsilon }$$

Given that ψ is influenced by the score functions s Y X | and s X independently, we have:

$$\nabla _ { s _ { X, Y } } \psi = \lim _ { \epsilon \to 0 } \frac { \psi ( \tilde { \mathcal { P } } _ { \epsilon } ( Y | X ) ) - \psi ( \mathcal { P } ( Y | X ) } { \epsilon } + \lim _ { \epsilon \to 0 } \frac { \psi ( \tilde { \mathcal { P } } _ { \epsilon } ( X ) ) - \psi ( \mathcal { P } ( X ) ) } { \epsilon }$$

Thus, we get:

$$\left ( \text{III.11} \right ) & & \nabla _ { s _ { X, Y } } \psi = \nabla _ { s _ { Y | X } } \psi + \nabla _ { s _ { X } } \psi$$

In Equations III.10 and III.11, we are actually factorizing the score function. An advantage of factorizing is that suppose the estimator is perturbed but only on one dimension after our appropriate factorization, then we can only calculate the change in the score function on that dimension and keep the others unchanged. We will use this advantage when we yield the efficient estimator later.

Indeed, we define the closure of the linear span of the score functions at P ( Z ) the tangent space of P ( Z ): T ( P ( Z )). As we showed in the factorization steps, we could use the set of the linear combination of the score functions to define the tangent space 32 :

$$\mathcal { T } ( \mathcal { P } _ { Z } ) = \left \{ h \colon h ( z ) = \sum _ { i } \alpha _ { i } s _ { \epsilon } ( z ), \quad \alpha _ { i } \in \mathbb { R }, \quad \epsilon \in [ 0, 1 ] \right \}$$

Since any score function belonging to the tangent space of the specific measure can be used during factorization, therefore, we are actually factorizing the tangent space. In the above bivariate example, we could write the factorization as: T X,Y = T Y X | ⊕T X , where T Y X | is defined as the tangent space associated with the conditional distribution P ( Y X | ): T Y X | = { h Y X | ( x, y ) : E h [ Y X | ( x, y ) | X = x ] = 0 for all x } . Correspondingly, T ( X ) denotes the tangent space associated with the marginal distribution P ( X ): T ( X ) = { h X ( x ) : E h [ X ( X )] = 0 for all x } . The symbol ⊕ denotes the direct sum, indicating that these spaces are orthogonal.

If the tangent space contains all the square-integrable and mean-zero functions (the same as the L 0 2 space, see Footnotes 31 and 32), the model which specifies the measure P ( Z ) : M P : ( Z ) ∈ M is a saturated model . For saturated models, the tangent space exists for all directions, meaning that if we move towards any direction, we are still in the model. Suppose our models are fully nonparametric; then, our model is saturated since no restrictions impede us. Otherwise, the model is not saturated. For instance, a general causal model is saturated, but an RCT causal model is not because the propensity in the randomized trial is fixed.

Finally, we could use the idea in the score function to define regularity. If we define ˜ P n := ˜ P ϵ =1 / √ n . Suppose ˜ P n is the empirical distribution of n samples drawn from the perturbed distribution ˜ P n . According to Equation III.1, regularity is defined as n increases, divergence between the empirical estimator ˆ ( ˜ ψ P n ) and the statistical estimand ψ ( ˜ P n ) converges in a fixed distribution. In this setting, with n increases, ϵ shrinks, and ˜ P ϵ gets closer to P , and the estimates are increasingly accurate and reflective of the true distribution as the sample size grows and the perturbation diminishes.

## C. Asymptotical Linearity and Influence Function

Contradicting to the property of regularity, asymptotical linearity is more intuitive. Multidimensional CLT has suggested that, the divergence between the sample mean of a series of IID random vectors ( ¯ X n = 1 n ∑ n k =1 X k ) and the expec-

32 Since the score functions are defined on the L 0 2 space, we could define the tangent space as the set of the square-integrable functions with respect to P ( Z ) whose means are 0:

$$\mathcal { T } ( \mathcal { P } ( Z ) ) = \left \{ h \in L _ { 2 } ( \mathcal { P } ( Z ) ) \, \colon E _ { \mathcal { P } ( Z ) } h ( x ) = 0 \right \}.$$

tation ( µ = E X [ ]) converges to the normal distribution with mean 0 and variance Σ = E X [( i -µ )( X i -µ ) T ] 33 .

$$\sqrt { n } ( \bar { X } _ { n } - \mu ) \stackrel { d } { \to } N ( 0, \Sigma )$$

Notice the similarities between Equations III.3 and III.12, the relationship between asymptotical linearity and asymptotical normality can't be more clear: in fact, if an estimator is asymptotically linear, then the CLT can be applied to the linear part of the estimator. We rewrite Equation III.2 into a linear form:

$$\psi ( \mathbb { P } _ { n } ) = \psi ( \tilde { \mathcal { P } } ) + \frac { 1 } { n } \sum _ { i = 1 } ^ { n } \phi ( \psi ; \mathbb { P } _ { n } ; Z _ { i } ) + o _ { p } ( n ^ { - 1 / 2 } )$$

33 Proving the multidimensional CLT requires prior knowledge on characteristic functions, Levy's Continuity Theorem, and Cramer-Wold device. The characteristic function of the real-valued random variable defines the probability distribution- in other words, two distinct distributions with the same characteristic function are identically the same distribution. Levy's continuity theorem states that if the characteristic functions of a sequence of random variables pointwise converge towards the characteristic function of a limiting random variable, then the sequence of random variables converges in distribution towards the limiting random variable. Cramer-Wold device states that a sequence of random variables converging to a limiting random variable is equivalent to the scalar of the sequence of random variables converging to the scalar of the limiting random variable. With these backgrounds, we have the proof:

We first define the normalized sample mean vector: given that ¯ X n = 1 n ∑ n i =1 X i , where X i are i.i.d. random variables with mean µ = E X [ i ] and covariance matrix Σ = E X [( i -µ )( X i -µ ) T ], let Z n = √ n X ( ¯ n -µ ). Then we find the characteristic function of Z n :

$$\varphi _ { Z _ { n } } ( t ) = E \left [ e ^ { i t ^ { T } Z _ { n } } \right ] = E \left [ e ^ { i t ^ { T } \sqrt { n } ( \hat { X } _ { n } - \mu ) } \right ] = E \left [ e ^ { i t ^ { T } \sqrt { n } ( \frac { 1 } { n } \sum _ { i \in 1 } ^ { n } X _ { i } - \mu ) } \right ]$$

Using the Taylor expansion of the exponential function e ix ≈ 1 + ix -x 2 2 , we get:

$$E \left [ e ^ { i t ^ { T } \frac { X _ { i } - \mu } { \sqrt { n } } } \right ] \approx 1 + i t ^ { T } \frac { E [ Y _ { i } ] } { \sqrt { n } } - \frac { 1 } { 2 } t ^ { T } \frac { E [ ( X _ { i } - \mu ) ( X _ { i } - \mu ) ^ { T } ] } { n } t$$

Since E X [ i -µ ] = 0 and E X [( i -µ )( X i -µ ) T ] = Σ, this simplifies to:

$$E \left [ e ^ { i t ^ { T } \frac { X _ { i } - \mu } { \sqrt { n } } } \right ] \approx 1 - \frac { 1 } { 2 } t ^ { T } \frac { \Sigma } { n } t$$

Therefore,

$$\varphi _ { Z _ { n } } ( t ) \approx \left ( 1 - \frac { 1 } { 2 } t ^ { T } \frac { \Sigma } { n } t \right ) ^ { n } \to e ^ { - \frac { 1 } { 2 } t ^ { T } \Sigma t } \, \text{as } n \to \infty$$

By Levy's Continuity Theorem, this implies that Z n converges in distribution to a multivariate normal distribution:

$$Z _ { n } \stackrel { d } { \rightarrow } \mathcal { N } ( 0, \Sigma )$$

Therefore, we have proved the multidimensional Central Limit Theorem:

$$\sqrt { n } ( \bar { X } _ { n } - \mu ) \stackrel { d } { \rightarrow } \mathcal { N } ( 0, \Sigma )$$

For the specific details, see van der Vaart 1998: Ch.2, pp.12-16.

The term 1 n ∑ n i =1 ϕ ψ ( ; P n ; Z i ) can be treated as a sum of i.i.d random variables (for instance, if the estimator is the sample mean ψ ( P n ) = ¯ X n , the term ϕ ψ ( ; P n ; X i ) = ( X i -µ ) so that we have Equation III.13 written as ¯ X n = µ + 1 n ∑ n i =1 ( X i -µ )). As we noted in Equation III.3, the variance of the term is the covariance for the converged normal distribution. In asymptotical linearity, the term is called influence function , as it shows the sensitivity of the estimator to small changes or perturbations in the data. It is named the influence function as it measures how each data point in the dataset ( Z i ∈ Z ) influences the overall estimate. If the estimator is unbiased (which for all the RAL estimators are true), the influence function should have a zero mean and finite variance (defined on the L 0 2 space, see Footnote 31).

Secondly, we notice that the convergence rate in both the CLT and the asymptotical linearity is 1 / √ n (as we introduced before, the term o p ( n -1 2 / is called the second-order remainder). The convergence rate in the CLT is understandable: the expectation and variance of the sum of the n random variables S n = ∑ n i =1 X i are separately E S [ n ] = nµ and Var( S n ) = nσ 2 . Therefore, during standardization, we have: S n -E S [ n ] √ Var( S n ) = S n -nµ √ nσ 2 . Since the standardization involves dividing by √ nσ 2 , we have the convergence rate n -1 2 / to balance the increase in total variability as the sample size increases. Similarly, in asymptotically linear settings, since 1 n ∑ n i =1 ϕ ψ ( ; P n ; Z i ) involves averaging the influence function over the sample, which is expected to converge to the normal distribution of mean and variance of ϕ ψ ( ; P n ; Z i ), and to balance the decrease in the asymptotical variance as the sample size n increases, the convergence rate n -1 2 / emerges naturally like we elaborated on the CLT case.

The influence function and the second-order remainder are the core parts of asymptotic analysis, and analogizing them with the CLT expression is only one perspective. The expression ϕ ψ ( ; P n ; Z i ) can not only be regarded as the influence function, it is also the pathway derivative, the gradient, and the Neyman orthogonal score.

Consider the core question of asymptotical analysis. Our target is to find the unbiased estimator ψ ( P n ) for the statistical estimand ψ ( P ) with the least statistical variance. Consider we have a continuous, smooth path defined as {P } ˜ ϵ starting from P ( Z ) = P where ϵ = 0 ending at ˜ ( P Z ) = ˜ where P ϵ = 1. We suppose ψ ( P ) is the estimand ( ψ ( P )), ψ ( ˜ ) is our estimator ( P ψ ( P n )). In the last session, we have shown that we could go from ψ ( ˜ ) to approach P ψ ( P ) by making ϵ → 0 (as n → ∞ ). Now, we consider the static 'snapshot,' which allows us to use ψ ( P ) to estimate ψ ( ˜ ). P An intuitive idea here is to construct the functional (distributional) Taylor expansion :

LEMMA III.2: Distributional Taylor Expansion: For a function f ( x ) which is differentiable in [ X ,X 0 1 ] , the distributional Taylor expansion can be expressed as:

$$( \text{III.1} ) \ f ( X _ { 0 } ) & \approx f ( X _ { 1 } ) + \nabla f ( X _ { 1 } ) ( X _ { 0 } - X _ { 1 } ) + \frac { 1 } { 2 } ( X _ { 0 } - X _ { 1 } ) \nabla ^ { 2 } f ( X _ { 1 } ) ( X _ { 0 } - X _ { 1 } ) ^ { T } \dots \\ o r,$$

$$f ( X _ { 0 } ) \approx f ( X _ { 1 } ) + \frac { \partial f ( x ) } { \partial x } \Big | _ { x = X _ { 1 } } ( X _ { 0 } - X _ { 1 } ) + \frac { 1 } { 2 } \frac { \partial ^ { 2 } f ( x ) } { \partial x ^ { 2 } } \Big | _ { x = X _ { 1 } } ( X _ { 0 } - X _ { 1 } ) ^ { 2 } \dots$$

We can apply Equation III.14 in our analysis. For the function ψ ( ˜ P ϵ ) in the domain ϵ ∈ [0 , 1], we have 34 : (III.15)

$$\psi ( \mathcal { P } ) \approx \underbrace { \psi ( \tilde { \mathcal { P } } ) } _ { \text{ naive plug-in estimator} } + \underbrace { \underbrace { \partial \psi ( \tilde { \mathcal { P } } _ { \epsilon } ) } _ { \text{ first-order bias correction} } } _ { \text{ first-order bias correction} } + \underbrace { \frac { 1 } { 2 } \frac { \partial ^ { 2 } \psi ( \tilde { \mathcal { P } } _ { \epsilon } ) } { \partial \epsilon ^ { 2 } } } _ { \text{ second-order remainder} } + \cdots$$

Therefore,

$$\psi ( \tilde { \mathcal { P } } ) - \psi ( \mathcal { P } ) \approx \underbrace { \frac { \partial \psi ( \tilde { \mathcal { P } } _ { \epsilon } ) } { \partial \epsilon } \Big | _ { \epsilon = 1 } } _ { \text{first-order bias correction} } - \underbrace { \frac { 1 } { 2 } \frac { \partial ^ { 2 } \psi ( \tilde { \mathcal { P } } _ { \epsilon } ) } { \partial \epsilon ^ { 2 } } \Big | _ { \epsilon = 1 } } _ { \text{first-order bias correction} }$$

To make the above decomposition more intuitive, we illustrate the ideas in Figure 2, which shows how to use results of ψ ( ˜ ) to approach P ψ ( P ). We will discuss the details as follows.

First-order bias correction and influence functions

We first analyze the first-order bias correction term. Due to the chain rule, we have:

## (III.16)

$$\stackrel { ( \Omega \colon \iota ) } { \partial \epsilon } \psi ( \tilde { \mathcal { P } } _ { \epsilon } ) \Big | _ { \epsilon = 1 } = \frac { \partial \psi ( \tilde { \mathcal { P } } _ { \epsilon } ) } { \partial \tilde { \mathcal { P } } _ { \epsilon } } \cdot \frac { \partial \tilde { \mathcal { P } } _ { \epsilon } } { \partial \epsilon } | _ { \epsilon = 1 } = \sum _ { i = 1 } ^ { n } \frac { \partial \psi ( \tilde { \mathcal { P } } _ { \epsilon } ) } { \partial \tilde { \mathcal { P } } _ { \epsilon } ( Z _ { i } ) } \cdot \frac { \partial \tilde { \mathcal { P } } _ { \epsilon } ( Z _ { i } ) } { \partial \epsilon } | _ { \epsilon = 1 } = \sum _ { i = 1 } ^ { n } \frac { \partial \psi ( \tilde { \mathcal { P } } _ { \epsilon } ) } { \partial p _ { \epsilon } ( Z _ { i } ) } \cdot ( p _ { 1 } ( Z _ { i } ) - p _ { 0 } ( Z _ { i } ) )$$

as p ϵ ( Z i ) refers to the probability mass function for measure P at the specific point Z i . The equation shows how each probability mass may affect the change in the estimator ψ ( ˜ P ϵ ). We may recall that we introduced the influence function because we would like to reveal how a slight perturbation at any data point in the dataset Z affects the estimator ψ ( ˜ P ϵ ). In order to understand the relationship between ∂ ∂ϵ ψ ( ˜ P ϵ ) and the influence function ϕ ψ ( ; P n , Z i ), we introduce the concept of Gateaux derivative :

34 This is also called 'von Mier Expansion' in some literature.

LEMMA III.3: Definition of the influence function with Gateaux derivative: The Gateaux derivative of ψ at P in the direction of another distribution G is defined as:

$$\mathbf D _ { \psi } ( \mathcal { P } ; \mathcal { G } ) = \lim _ { \epsilon \to 0 } \frac { \psi ( ( 1 - \epsilon ) \mathcal { P } + \epsilon \mathcal { G } ) - \psi ( \mathcal { P } ) } { \epsilon }$$

Now suppose for distribution P we have a small perturbation at a point z , we write the perturbed distribution as:

$$\tilde { \mathcal { P } } _ { \epsilon } = \mathcal { P } + \epsilon ( \delta _ { z } - \mathcal { P } ) \Rightarrow \psi ( \tilde { \mathcal { P } } _ { \epsilon } ) = \psi ( \mathcal { P } + \epsilon ( \delta _ { z } - \mathcal { P } ) )$$

Where δ z stands for the point mass distribution at z 35 . Now that the Gateaux derivative of ψ at P in the direction of δ z is:

$$\text{(III.18)} & & \text{D} _ { \psi } ( \mathcal { P } ; \delta _ { z } ) = \lim _ { \epsilon \rightarrow 0 } \frac { \psi ( \mathcal { P } + \epsilon ( \delta _ { z } - \mathcal { P } ) ) - \psi ( \mathcal { P } ) ) } { \epsilon }$$

Gateaux derivative expressed in Equation III.18 is the mathematical definition of the influence function of ϕ ψ ( ; P ), as it describes the impact of a small change in the data on the given functional ψ . Since the influence function describes an infinitesimal perturbation in the distribution (measure) ˜ P ϵ at point z , we could regard the influence function as a dichotomous representation of the functional derivative δψ ( ˜ P ϵ ) δp ˜ ( ϵ z ) , or the functional derivative as the influence function in the context of continuous distributions. We could have some mathematical transformation with the definition to reveal its relationship with the first-order bias. We first use partial differentiation to express the equation:

$$\phi ( \psi ; \tilde { \mathcal { P } } ; z ) = \lim _ { \epsilon \rightarrow 0 } \frac { \psi ( \tilde { \mathcal { P } } + \epsilon ( \delta _ { z } - \tilde { \mathcal { P } } ) ) - \psi ( \tilde { \mathcal { P } } ) ) } { \epsilon }$$

Integral the influence function on the direction of P - P ˜ , we have:

$$( \text{III.19} ) \quad \ \frac { \partial \psi ( \tilde { \mathcal { P } } + \epsilon ( \mathcal { P } - \tilde { \mathcal { P } } ) ) } { \partial \epsilon } \Big | _ { \epsilon = 0 } = \int \phi ( \psi ; \tilde { \mathcal { P } } ; z ) ( p ( z ) - \tilde { p } ( z ) ) d z$$

Notice on the left side, the numerator is actually ˜ + P ϵ ( P - P ˜ ) = ˜ , the right P ϵ side of Equation III.16, which has a similar structure. Therefore,

$$\frac { \partial } { \partial \epsilon } \psi ( \tilde { \mathcal { P } } _ { \epsilon } ) \Big | _ { \epsilon = 1 } = - \int \phi ( \psi ; \tilde { \mathcal { P } } ; z ) ( p ( z ) - \tilde { p } ( z ) ) d z = - \int \phi ( \psi ; \tilde { \mathcal { P } } ; z ) p ( z ) d z$$

35 Generally, we call δ z the Dirac delta function at z . However, it is not a traditional function; rather, it is a distribution δ a which is concentrated at a single point of z . In some literature, introducing the Dirac delta function is a trick for 'point mass contamination', which introduces a particular kind of path where the destination is a distribution that places all of its mass at the point z .

$$= - \frac { 1 } { n } \sum _ { i = 1 } ^ { n } \phi ( \psi ; \tilde { \mathcal { P } } ; Z _ { i } )$$

Because ∫ ϕ ψ ( ; ˜ ; P z p z dz )˜( ) = 0 36 . We are very familiar with the right side of the equation: this is the expectation of the influence function. This is to say, the influence function can be regarded as a pathwise derivative, and we could approximate ψ ( P ) with the naive plug-in estimator ψ ( P n ) plus the expectation of the influence function for P n ). We define this estimator as the 'one-step' estimator (Fisher and Kennedy 2021):

$$\hat { \psi } ( \mathbb { P } _ { n } ) _ { 1 \text{-step} } \approx \psi ( \mathbb { P } _ { n } ) + \sum _ { i = 1 } ^ { n } \phi ( \psi ; \mathbb { P } _ { n } ; Z _ { i } )$$

The Gateaux derivative also reveals the relationship between the gradient to the score function and the influence function. Recall our definition of the gradient in the direction of the score function on the estimator, we have:

$$\nabla _ { s } \psi ( \mathcal { P } ) = \lim _ { \epsilon \rightarrow 0 } \frac { \psi ( \tilde { \mathcal { P } } _ { \epsilon } ) - \psi ( \mathcal { P } ) } { \epsilon } = \frac { \partial } { \partial \epsilon } \psi ( \tilde { \mathcal { P } } _ { \epsilon } ) \Big | _ { \epsilon = 0 }$$

Similarly, using distributional Taylor expansion, we have:

$$\psi ( \tilde { \mathcal { P } } _ { \epsilon } ) \Big | _ { \epsilon = 0 } \approx \psi ( \mathcal { P } ) + \epsilon \ \frac { \partial \psi } { \partial \epsilon } \Big | _ { \epsilon = 0 } + O ( \epsilon )$$

And use the chain rule to analyze the first-order term :

$$\frac { \partial \psi } { \partial \epsilon } \Big | _ { \epsilon = 0 } = \int \frac { \delta \psi } { \delta p ( z ) } \ \frac { \partial \tilde { p } _ { \epsilon } ( z ) } { \partial \epsilon } \Big | _ { \epsilon = 0 } d z.$$

Based on Equation III.7 describing the relationship between p z ( ) and ˜ ( p ϵ z ):˜ ( p ϵ z ) = (1 + ϵs z ( )) p z ( ), we have:

$$\frac { \partial \tilde { p } _ { \epsilon } ( z ) } { \partial \epsilon } \Big | _ { \epsilon = 0 } = \frac { \partial } { \partial \epsilon } \left ( 1 + \epsilon s ( z ) \right ) p ( z ) \Big | _ { \epsilon = 0 } = s ( z ) p ( z ).$$

36 Similarly, we could also have the derivation: ∂ ∂ϵ ψ ( ˜ P ϵ ) ∣ ∣ ∣ ∣ ϵ =0 = 1 n ∑ n i =1 ϕ ψ ( ; P , Z i ). Notice the left side of the equation, ∂ ∂ϵ ψ ( ˜ P ϵ ) ∣ ∣ ∣ ∣ ϵ =0 , is indeed the gradient on the direction of the score function on the estimator ∇ s ψ ( P ).

Therefore,

$$\dots \frac { d \psi } { d \epsilon } \Big | _ { \epsilon = 0 } = \int \frac { \delta \psi } { \delta p ( z ) } s ( z ) p ( z ) \, d z. \quad \dots \quad.$$

Note the right side of the expression can be expressed as the form of expectation under the original distribution p z ( ):

$$\int \frac { \delta \psi } { \delta p ( z ) } s ( z ) p ( z ) \, d z = \mathbb { E } _ { \mathcal { P } } \left [ \frac { \delta \psi } { \delta p ( z ) } s ( z ) \right ].$$

As we mentioned, the form δψ P ( ) δp z ( ) can be regarded as the continuous form of the influence function ϕ ψ ( ; P ; z ) 37 . Therefore, we finally obtain:

$$\left ( \text{III.21} \right ) & & \nabla _ { s } \psi ( \mathcal { P } ) = \mathbb { E } _ { \mathcal { P } } \left [ \phi ( \psi ; \mathcal { P } ; z ) s ( z ) \right ].$$

Equation III.21 is called the central identity for influence functions . This equation reveals that the gradient of the functional ψ with respect to the score function s z ( ) of the estimator is given by the expectation of the product of the influence function ϕ ψ ( ; P ; z ) and the score function s z ( ) under the distribution P . This relationship indicates that the direction of the gradient is aligned with the interaction between the influence function and the score function (as the expectation of the score and the influence functions under RAL are both zero). This equation is the most useful tool for extracting the influence function or validating if the influence function is correct, especially for efficient influence functions for the saturated models. We will give a specific example of the usage in the next subsection.

Also, since the influence functions Φ are indeed a gradient of the estimator ψ , we may have the 'gradient algebra' rules for the influence functions. For instance, the chain rule: ϕ g ψ ( ( )) = g ′ ( ψ ϕ ψ ) ( ); and the product rule: ϕ ψ ψ ( 1 2 ) = ϕ ψ ( 1 ) ψ 2 + ψ ϕ ψ 1 ( 2 ). Similar rules are also applied to the score functions and the gradient on the score functions.

Finally, the influence function can be regarded as the Neyman's orthogonal score . Originally, Neyman orthogonality refers to a score function that yields robust estimations for small perturbations. Neyman's orthogonal score is defined as a function f ( Z ψ ϵ ; ˆ ; ˆ ) ( Z denotes the measurable set, ˆ denotes the estimator of ψ interest, and ˆ is the perturbation) which is orthogonal to the perturbation when ϵ

37 In some literature, the definition of the influence function with the expectation expression is the score-based definition of the influence function, and the score function can also be regarded as s 0 ( z ) = ∂ ∂ϵ log[ p z ( ) + ϵ p (˜ ( ϵ z ) -p z ( ))] ∣ ∣ ϵ =0 = ˜ ( p ϵ z ) -p z ( ) p z ( ) , which is the same as the definition in Equation III.6.

we are at the true (population) value (Chernozhukov et al. 2018):

$$\frac { \partial } { \partial \epsilon } E [ f ( Z ; \hat { \psi } ; \hat { \epsilon } ) | _ { \hat { \epsilon } = 0 } ] = 0$$

.

The influence function 38 is intrinsically equivalent to Neyman's orthogonal score simply because it is the mean-zero score function which measures the gradient on the estimator towards the perturbation. Recall we have defined the gradient of an estimator on the direction of the score function (Equation III.11): ∇ s ψ ( P ) = lim ϵ → 0 ψ ( ˜ P ϵ ) -ψ ( P ) ϵ . Now suppose s z denotes the score function of the perturbed distribution δ z and ˜ P ϵ = (1 -ϵ ) P + ϵδ z . Thus, we have: (III.22)

$$\overset { \theta } { \nabla } _ { s _ { z } } \psi ( \mathcal { P } ) = \lim _ { \epsilon \rightarrow 0 } \frac { \psi ( \mathcal { P } + \epsilon ( \delta _ { z } - \mathcal { P } ) - \psi ( \mathcal { P } ) } { \epsilon } = \frac { \partial } { \partial \epsilon } [ \psi ( \mathcal { P } + \epsilon ( \delta _ { z } - \mathcal { P } ) ] | _ { \epsilon = 0 } = \phi ( \psi ; \mathcal { P } ; z )$$

Obviously, the right side of the equation is the Gateaux derivative definition of the influence function. Moreover, due to the central identity of the influence function, we know the gradient in the direction of s z should be equivalent to the covariance of the corresponding score function and the influence function. Thus, we get:

$$\nabla _ { s _ { z } } \psi ( \mathcal { P } ) = \phi ( \psi ; \mathcal { P } ; z ) = E [ \phi ( \psi ; \mathcal { P } ; z ) s _ { z } ( Z ) ]$$

We bring up the perspective of Neyman's orthogonality on the influence function since most of the econometric literature on doubly/debiased machine learning starts with Neyman's orthogonality to find the robust estimator of the parameter of interests (for instance, the estimator on the ATE), which is the basis of the doubly/debiased machine learning method(Chernozhukov et al. 2018a; Chernozhukov et al. 2018b). Indeed, we illustrate here that the robust/debiased estimator based on the Neyman orthogonality is mathematically intrinsic to the efficient estimator based on the influence function. As for the double robustness, they apply the sample-splitting method to ensure the asymptotical validity of the convergence process, which we will discuss below.

## Second-order remainder

Once we develop the one-step estimator from the first-order bias correction term, we would like to evaluate the quality of the estimatorwhether it is a 'good' proxy of the true estimand with observational data on hand. Generally, the quality of the estimator is decided by several factors: the sample size (which determines the effect of bias correction since we sum up the influence functions), the smoothness of the functional form (to be discussed below), and the choice of

38 The population true values in general regarded as ϵ 0 , and therefore, the condition is ˆ = ϵ ϵ 0 .

the measure ˜ (when P ϵ = 1).

To introduce the idea of smoothness of the functional form, we turn our focus to the second-order remainder, which is the gap between the one-step estimator and the true estimand. We denote it as R 2 ( P n , P ) = ψ ( P ) -ψ ( P ) one-step . This part is quite similar to the loss function in machine learning terms, as it measures the difference between the predictor (estimator) ψ ( P n )and the true value (estimand) ψ ( P ).

In Figure 2, we simply project the continuous measures of ˜ P ϵ on the x-axis and take two points on the axis separately as ˜ P and P . Therefore, the direction of perturbation parallels to the x-axis. We may also imagine that starting from P , we may have infinite choices of directions (spanning the tangent space), and the hyperplane ψ intersects with each dimension along ϵ . So, we have estimators ψ ( ˜ ) P in all dimensions. If we would like to compare the smoothness of the estimators stretching from different dimensions over different distances, we need a process to 'standardize' the measures onto one dimension. Therefore, we come up with a standardized measure P :

$$\mathbf P \coloneqq \mathcal { P } + \frac { \Delta } { | | \tilde { \mathcal { P } } - \mathcal { P } | | _ { 2 } } ( \tilde { \mathcal { P } } - \mathcal { P } )$$

where ||P-P|| ˜ 2 is the L 2 distance between the two measures ˜ and P P : ||P-P|| ˜ 2 = √ ∫ (˜( p z ) -p z ( )) 2 dz . As we can tell from the definition, ∆ shows the absolute distance from the standardized measure to the original measure: ∆ = || P -P|| 2 . The standardization process makes it possible to compare and analyze the smoothness of the estimators uniformly across different dimensions and distances shown below.

With the definition of the standardized measure, we could define the j -th order smoothness as the j -th order derivative of the estimator ψ ( P (if it is j -th order differentiable) with respect to the absolute distance ∆ converges to a constant number if ∆ approaches zero:

$$\frac { \partial ^ { j } } { \partial \Delta ^ { j } } \psi ( \mathbf P ) | _ { \Delta \rightarrow 0 } = O ( 1 ).$$

Fisher and Kennedy (2021) have a very accurate metaphor to help understand the idea of 'smoothness': it is like the 'magician tablecloth trick.' If the estimators are more smooth, changing the parameter ∆ (like pulling the tablecloth) will have less impact on the estimator. Since many (nonparametric) machine learning models use gradient-based methods for optimization, the smooth differentiable property of the estimator (and the convexity of the function) ensures the model converges to the global minimum, and facilitates the gradient optimization.

Specifically for the second-order smoothness, it allows the functional form of the estimator ψ to change nonlinearly with the perturbation, which is to say, the second derivative of ψ with respect to ∆ exists and is bounded. If the second-order smoothness is allowed, usually researchers presume that the linear approximations are insufficient to capture the true estimand. However, adding the term will help improve the error bound of our estimator. Accounting for the curvature of the function:

$$\left | \psi ( \tilde { P } _ { \epsilon } ) - \psi ( P ) - \psi ^ { \prime } ( P ) ( \tilde { P } _ { \epsilon } - P ) - \frac { 1 } { 2 } \psi ^ { \prime \prime } ( P ) ( \tilde { P } _ { \epsilon } - P ) ^ { 2 } \right | = O ( \epsilon ^ { 2 } )$$

the O ϵ ( 2 ) term indicates that the error decreases quadratically with the size of the perturbation, providing a much tighter bound compared to the linear O ϵ ( ) term from the first-order approximation:

$$\left | \psi ( \tilde { P } _ { \epsilon } ) - \psi ( P ) - \psi ^ { \prime } ( P ) ( \tilde { P } _ { \epsilon } - P ) \right | = O ( \epsilon )$$

However, in the efficient estimation we are discussing in this thesis, we assume that the second-order remainder term becomes asymptotically negligible faster than the rate of 1 / √ n . In other words, as the sample size n increases, R 2 diminishes faster than 1 / √ n and makes its impact on the overall estimation asymptotically insignificant:

ASSUMPTION III.1 (Negligiblity of the Second-Order Remainder): For the secondorder remainder R 2 , we have the assumption that:

$$R _ { 2 } = o _ { p } \left ( \frac { 1 } { \sqrt { n } } \right ).$$

Which indicates:

$$\forall \epsilon > 0, \ \forall c > 0, \ \Pr \left ( | R _ { 2 } | > \frac { c } { \sqrt { n } } \right ) \to 0 \ \ a s \ \ n \to \infty.$$

If the Assumption is violated, which suggests that the R 2 term converges no faster than the rate of 1 / √ n , meaning that the R 2 term cannot be regarded as insignificant as the sample size grows (under large sample size), so the first-step estimator may not be asymptotically valid to be the approximation on the true estimand. Under such a scenario, sensitivity analyses or higher-order correction terms are needed (for instance, adding the second-order smoothness term as we elaborated above). Luckily, many machine learning models, including but not restricted to neural networks, highly adaptive lasso, gradient boosting, nearest neighbors, and regressive methods, have the property of consistency in the prediction, and we do not need to worry about the convergence rate for the second-order remainder

## Empirical Process Convergence

The second-order remainder term is 1 / √ n -negligible suggests that we would like the empirical data distribution, as the sample size increases, to converge to the true distribution at the pace faster than 1 / √ n for all other parts than the one-step estimator. However, in Equation III.20, we may find that our influence function on the second term of the right side of the equation is also based on its empirical (sample) mean, rather than the true (population) expectation. The process of converging the empirical influence function towards the 'true' influence function is called empirical process , as we decompose the difference between the onestep estimator and the true estimand, excluding the R 2 terms which under our assumption is negligible:

$$\text{assumption is negligible} \\ \hat { \psi } ( \mathbb { P } _ { n } ) _ { 1 - \text{step} } - \psi ( \mathbb { P } ) & = \sum _ { i = 1 } ^ { n } \phi ( \psi ; \mathbb { P } _ { n } ; Z _ { i } ) \\ & = \sum _ { i = 1 } ^ { n } \left [ \left ( \phi ( \psi ; \mathbb { P } _ { n } ; Z _ { i } ) - \phi ( \psi ; \mathbb { P } ; Z _ { i } ) \right ) + \left ( \phi ( \psi ; \mathbb { P } ; Z _ { i } ) - E _ { \mathbb { P } } [ ( \phi ( \psi ; \mathbb { P } _ { n } ; Z ) ] \right ) \\ & \quad + \left ( E _ { t i l d e \mathcal { P } } [ ( \phi ( \psi ; \mathbb { P } ; Z ) ] - E _ { \mathbb { P } } [ ( \phi ( \psi ; \mathbb { P } _ { n } ; Z ) ] \right ) \right ] \\ & = \sum _ { i = 1 } ^ { n } \phi ( \psi ; \mathbb { P } _ { n } ; Z _ { i } ) - E _ { \mathbb { P } } \phi ( \psi ; \mathbb { P } _ { n } ; Z ) ] \\ As the second term represents the influence function from the \text{true} \text{simand,} it$$

As the second term represents the influence function from the true estimand, it should be zero due to the property of the influence function. This equation thus indicates that as the empirical measure P n converges to the true distribution ˜ P the difference between the empirical influence function ϕ ψ ( ; P n ; Z i ) and the true influence function ϕ ψ ( ; ˜ ; P Z ) approaches zero 40 . In other words, the empirical influence function ϕ ψ ( ; P n ; Z i ), which serves as an estimator for the true influence function ϕ ψ ( ; ˜ ; P Z ), must also have the property of consistency 41 as the sample size increases. This consistency ensures that the empirical influence function accurately reflects the true influence function as the number of samples grows, thereby contributing to the overall convergence of the one-step estimator to the true parameter.

39 The property is called L 2 consistency , which means that as the sample size increases, prediction error of the model (measured in terms of L 2 norm or simply MSE) converges to the lowest possible prediction error for the underlying data distribution.

40 Recall that in our definition of the influence function for the empirical measure P n = ˜ P ϵ approaches the influence function of the measure ˜ (which P ϵ = 1). This is consistent with our Figure 2 illustration on the Distributional Taylor Expansion, as we only have the derivative at ϵ = 1 end.

41 This is specifically L 2 consistency.

With the influence function of the empirical data converging to the influence function of the estimated, we need statistical methods to control the speed of the convergence at 1 / √ n so that our one-step estimator can be asymptotically valid and become the unbiased estimator of the true estimand in large sample size case. To ensure the convergence speed is no slower than 1 / √ n , we either need sample splitting techniques to process the original data, or we could assume that the empirical influence function falls into the Donsker class 42 . In this thesis, we only use sample-splitting techniques to ensure that our estimator is valid on the estimand.

Sample splitting is the most common technique to control the speed of the empirical process and make the efficient estimator a valid proxy for the true estimand. In machine learning, the method is also called cross-validation . Moreover, since we need to split and fit the model in at least two sub-datasets, this process is also called double machine learning , and the estimator with this strategy is the double machine learning (DML) estimator (intrinsically, it is the efficient estimator).

We can prove the effect of sample splitting on controlling the convergence rate of the empirical influence function towards the influence function from the estimand with Theorem III.12. To make the case simple, consider splitting the dataset ( Z , ..., Z 1 n ) only into two subsets: S 1 = { Z , Z , . . . , Z 1 2 n 1 } be the first subset, and S 2 = { Z n 1 +1 , Z n 1 +2 , . . . , Z n } be the second subset.

We first use S 1 to construct an initial estimator ˆ ψ 0 of the parameter ψ . Which, presumably, is the consistent estimator on the estimand ˆ . ψ We may write the one-step estimator with sample splitting as:

$$\hat { \psi } _ { 1 \text{-step} } = \hat { \psi } _ { 0 } + \frac { 1 } { n _ { 2 } } \sum _ { i \in S _ { 2 } } \phi ( \hat { \psi } _ { 0 }, \mathbb { P } _ { n _ { 2 } }, Z _ { i } )$$

Notice that for the naive estimator, we use the result from S 1 , while for the firstorder bias correction term, we use estimators from S 2 (as n 2 is the size of the split S 2 , and P n 2 = 1 n 2 ∑ i ∈ S 2 δ Z i is the empirical measure based on S 2 ). Obviously, the one-step estimator is consistent, as the law of large numbers indicates ˆ ψ 0 prob. - - -→ ψ and the correction term 1 n 2 ∑ i ∈ S 2 ϕ ψ , ( ˆ 0 P n 2 , Z i ) will converge in probability to zero

42 Aclass of functions Φ is a Donsker class if the empirical process indexed by Φ converges in distribution to a Gaussian process: G n (Φ) = √ n ( P n -P )Φ, as P n denotes the empirical measure and P is the true underlying probability measure. The Donsker class has the property G n ⇝ G , converging to a Gaussian Process. In some ways, we could rewrite the empirical process as √ n ( P n -P )( ˆ Φ -Φ), if we regard the convergence of the empirical measure P n and the convergence of the influence function Φ together at ˆ the rate of 1 / √ n . Obviously, the Donsker class property ensures that the empirical process converges uniformly and at a controlled rate. When the influence function of an estimator belongs to a Donsker class, this guarantees that the empirical process does not exhibit erratic behavior and converges smoothly.

if ϕ ψ, ( P n , Z i ) is well-defined. Meanwhile, since the naive estimator and the firstorder bias correction term from independent datasets, the construction of ˆ ψ 0 is irrelevant to the data used to make the correction.

Now we turn to the asymptotical normality. Due to the CLT, we have:

$$\frac { 1 } { \sqrt { n _ { 2 } } } \sum _ { i \in S _ { 2 } } \phi ( \psi, \mathbb { P } _ { n _ { 2 } }, Z _ { i } ) \stackrel { d } { \rightarrow } \mathcal { N } ( 0, \sigma ^ { 2 } )$$

Since ˆ ψ 0 prob. - - -→ ψ and ϕ ψ , ( ˆ 0 P n 2 , Z i ) ≈ ϕ ψ, ( P n 2 , Z i ) , we have:

$$\frac { 1 } { \sqrt { n _ { 2 } } } \sum _ { i \in S _ { 2 } } \phi ( \hat { \psi } _ { 0 }, \mathbb { P } _ { n _ { 2 } }, Z _ { i } ) \xrightarrow { d } \mathcal { N } ( 0, \sigma ^ { 2 } ) \quad ( \text{Slutsky} ^ { \prime } s \text{ Theorem} )$$

Therefore, the one-step estimator with sample splitting method is asymptotically normal with mean ψ and variance σ /n 2 2 , and the control on the convergence rate at 1 / √ n is accomplished.

In summary, with the second (and higher) order remainders and the empirical process becomes insignificant under the large sample size scenario, we can use the one-step estimator: the empirical naive estimator plus the first-order bias correction term to robustly and efficiently infer the true estimand.

## D. Efficient Influence Functions and Efficient Estimators

## Efficient influence functions

So far, we have discussed the score function, the influence function, and the estimator from the general RAL conditions. With the tools we have, we could have proof of the CRB brought up at the start of this section:

## PROOF ON CRAMER-RAO BOUND:

Since our estimator ˆ = ψ ψ ( P n ) is an unbiased estimator of ψ = ψ ( P ), therefore, E ψ [ ˆ ] = ψ . Our target is to get the upper bound for the variance of ˆ . ψ Recall that we could use Cauchy inequality to grab the upper bound for the variance if the expectation is known: E XY [ ] 2 ≤ E X [ ] 2 E Y [ ] 2 , we use the score function pointing to ψ (from ˆ ) ψ as the auxiliary function, as the expectation of the score function is zero. We denote the score function as s Z ( ; ψ ) = ∇ ψ (˜( p z ) = ∂ log ˜ ( p ϵ z ) ∂ϵ ∣ ∣ ϵ =0 to show the effect of ψ . Therefore,

$$c o v ( \hat { \psi } ; S ( Z ; \psi ) ] ^ { 2 } \leq \text{Var} ( \hat { \psi } ) E [ S ( Z ; \psi ) ] ^ { 2 } \Rightarrow \text{Var} ( \hat { \psi } ) \geq \frac { c o v ( \hat { \psi } ; s ( Z ; \hat { \psi } ) ) } { E [ S ( Z ; \psi ) ^ { 2 } ] }$$

Since we define the score function on the direction to ˆ , ψ the two components in the covariance ˆ and ψ s Z ( ; ψ ) are indeed in the same direction, and therefore the covariance is 1, we denote the expectation of the square of the score function as the Fisher information: I ( ψ ) = E s Z ψ [ ( ; ) 2 ]; therefore, we have the Cramer-Rao bound for the variance of the unbiased estimator ˆ : ψ

$$\text{Var} ( \hat { \psi } ) \geq \frac { 1 } { I ( \psi ) }$$

As we noted earlier in this chapter, if the estimator ˆ attains the Cramer-Rao ψ bound, we call the estimator the efficient estimator. Meanwhile, its associated influence function is called the effient influence function (EIF) . We can either have the EIF to get the efficient estimator, or reversely obtain the efficient estimator and get its EIF. Usually, we adopt the first strategy: we first find a set Φ of possible influence functions and select the most efficient one which satisfies ϕ † = arg min ϕ ∈ Φ Var[ ϕ ] 43 and we finally get the efficient estimator.

A crucial property for the EIF ϕ † is that it lies in the tangent space. This is because only the influence function in the tangent space obtains the lowest variance: ϕ † ⊥ h ⊥ , where h ⊥ stands for the lines orthogonal to any element in the tangent space h ∈ T .

## PROOF:

We suppose we could decompose the influence function into two orthogonal parts: one is the projection on the tangent space h , and one is orthogonal to the tangent space, ϕ -h . For any score function s on the tangent space, according to the definition, we have: E ϕ [( -h s ) ] = 0 44 . Thus, for the variance of the influence function:

$$\text{Var} ( \phi ) & = \text{Var} ( h + ( \phi - h ) ) = \text{Var} ( h ) + \text{Var} ( \phi - h ) \text{ (orthogonality, no covariance)} \\ & = E [ h ^ { 2 } ] - ( E [ h ] ) ^ { 2 } + E [ ( \phi - h ) ^ { 2 } ] + ( E [ \phi - h ) ] ) ^ { 2 } \\ & = E [ h ^ { 2 } ] + E [ ( \phi - h ) ^ { 2 } ]$$

Since h is constant (as it is the projection of the influence functions onto the tangent space), the target function ϕ † = arg min ϕ ( ϕ -h ) 2 ⇒ ϕ † = h . Therefore,

43 It's not hard to understand that the lowest variance for the estimator indicates the lowest variance (square mean) for the influence function so I didn't include it in the main text: we could just get it from Equation III.4 on the definition of asymptotical linearity; or we could still use the Taylor Distributional expansion to understand it, as ˆ is ψ an unbiased estimator for ψ , so the variance only remain in the higher-order derivatives and the increase in variance leads to the increase in the derivatives, which leads the increase in the variance in the influence functions.

44 The influence function and the score function are both Hilbert space (since they are both L 0 2 space as we mentioned before, and we could also use the inner product of Hilbert space to represent it, ⟨ ϕ -h, s ⟩ = 0 , ∀ s ∈ T ( P ). The orthogonal decomposition is justified by the Reisz Representation Theorem, which states that for every continuous linear functional F on the Hilbert space H , there is always an element hin H such that F ( f, h ) = ⟨ f, h ⟩ ∀ f ∈ H .

the influence function should lie on the tangent space to be efficient.

Therefore, we call the EIF the ' canonical gradient .' With this property or the EIF, we have another direct corollary: for the saturated models, there is one influence function, which is the EIF. This is simply because, the score functions of the saturated model point to all directions and therefore could take any values satisfying zero mean and finite variance. Suppose we have two influence functions ϕ 1 and ϕ 2 , therefore, E ϕ s [ 1 ] = E ϕ s [ 2 ] = E ϕ [( 1 -ϕ 2 ) s = 0] ⇒ ϕ 1 = ϕ 2 . However, this is not the case for the non-saturated models, as the restrictions of the nonsaturated models make some directions orthogonal to the tangent space possible.

## Derving EIF for the ATE

The review above in this section gives a rough (but lengthy) introduction to the efficient theory and the efficient estimator, and all the estimators in this thesis are efficient (whose estimation errors reach the CRB or, at least, the most efficient one). The remaining job in the introductory chapter is to derive the efficient estimator for the average treatment effect, which we have discussed in Section II: Y (1) -Y (0).

We start with the saturated model: the ATE in the observational study, where there is no restriction on the statistical estimand to infer the causal estimand, and there's only one influence function- which is the EIF for the estimator. Similar to the operation in Equation II.2, we use ψ a = E E [ X [ Y A | = a, X ]] from the observational study to infer ψ ∗ a = E Y [ ( a )] (therefore, the ATE is ψ 1 -ψ 0 ). From the perspective of the efficient theory, the estimators we applied in Section II, for instance, the IPW estimator, as we have shown, is a RAL estimator but not the efficient one, as it is only the 'naive plug-in estimator' part in the Taylor distributional expansion decomposition (Equation III.15). Our goal is the efficient estimator for ψ a :

$$\psi _ { a } = E [ E _ { X } [ Y | A = a, X ] ] = \sum _ { x } \underbrace { E [ Y | A = a, X ] } _ { \colon = \mu _ { a } ( X ) } p ( x ) = \sum _ { x } \mu _ { a } ( x ) p ( x )$$

As we discussed above, to derive the efficient estimator, we first capture the EIF and then derive the efficient estimator based on the EIF. We need to calculate

$$( \text{III.23} ) \quad \phi ( \psi _ { a } ) = \sum \phi ( \mu _ { a } ( x ) p ( x ) ) = \sum \left ( \phi ( \mu _ { a } ( x ) ) p ( x ) + \mu _ { a } ( x ) \phi ( p ( x ) ) \right ).$$

As the equation shows, we need the EIFs for the expectation ϕ E X ( [ ]) (to yield ϕ p x ( ( )))and the conditional expectation ϕ E Y X ( [ | = x ]) (to yield ϕ µ ( a ( x )).

We first derive the EIF for ψ ( P ) = E P [ X ] = ∫ xd P ( x . ) We start with the

Gateaux derivative definition of the influence function. The perturbed distribution is set as:

$$\tilde { \mathcal { P } } _ { \epsilon } = ( 1 - \epsilon ) \mathcal { P } + \epsilon \delta _ { x }$$

Which is similar to what we did in Equation III.18. Therefore, we construct the nominator ψ ( ˜ P ϵ ):

$$\ddot { \psi } ( \tilde { \mathcal { P } } _ { \epsilon } ) & = \int x d \tilde { \mathcal { P } } _ { \epsilon } = \int x d ( ( 1 - \epsilon ) \mathcal { P } + \epsilon \delta _ { x } ) \\ & = ( 1 - \epsilon ) \int x d \mathcal { P } + \epsilon \int x d \delta _ { x } \\ & = ( 1 - \epsilon ) \psi ( \mathcal { P } ) + \epsilon x \\ & = \psi ( \mathcal { P } ) + \epsilon ( x - \psi ( \mathcal { P } ) ) \\ \cdot \, \dashrightarrow \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \quad \cdot \, \end{cases}$$

Since ∫ dδ x = 1. Thus, we have the influence function ϕ ψ ( = E P [ X ]; P , x ) :

$$\phi ( E _ { \mathcal { P } } [ X ] ) = \frac { \partial } { \partial \epsilon } \psi ( \tilde { \mathcal { P } } _ { \epsilon } ) | _ { \epsilon = 0 } = x - \psi ( \mathcal { P } ) = x - E _ { \mathcal { P } } [ X ]$$

We then derive the EIF for ψ ( P ) = E P [ Y X | = x ] = ∫ y yd P ( y x | ). Still, ˜ P ϵ = (1 -ϵ ) P + ϵδ y x | . The hard part we need to derive here is P ( y x | ). We may recall the Bayesian rule P ( y x | ) = P ( y,x ) P ( x ) . Therefore,

$$\tilde { \mathcal { P } } _ { \epsilon } ( y | x ) = \frac { \tilde { \mathcal { P } } _ { \epsilon } ( y, x ) } { \tilde { \mathcal { P } } _ { \epsilon } ( x ) } = \frac { ( 1 - \epsilon ) \mathcal { P } ( y, x ) + \epsilon \delta _ { y, x } } { ( 1 - \epsilon ) \mathcal { P } ( x ) + \epsilon \delta _ { x } }$$

Given by the equation above. Since the influence function is defined as ∂ψ ( ˜ P ϵ ) ∂ϵ ∣ ∣ ϵ =0 is a gradient, we first deal with the part to be integrated ( d ˜ ( P ϵ y x | )), which we have transformed with the above equation. Recall the gradient algebra rule: for two functions u and v , ( u v ) ′ = u v ′ -uv ′ v 2 . Therefore,

$$\tilde { P } _ { \epsilon } ( y | x ) | _ { \epsilon = 0 } & = \frac { [ ( 1 - \epsilon ) \mathcal { P } ( y, x ) + \epsilon \delta _ { y, x } ] ^ { \prime } [ ( 1 - \epsilon ) \mathcal { P } ( x ) + \epsilon \delta _ { x } ] - [ ( 1 - \epsilon ) \mathcal { P } ( y, x ) + \epsilon \delta _ { y, x } ] [ ( 1 - \epsilon ) \mathcal { P } ( x ) + \epsilon \delta _ { x } ] ^ { \prime } } { [ ( 1 - \epsilon ) \mathcal { P } ( x ) + \epsilon \delta _ { x } ] ^ { 2 } } \\ & = \frac { [ \delta _ { y, x } - \mathcal { P } ( y, x ) ] [ ( 1 - \epsilon ) \mathcal { P } ( x ) + \delta _ { x } ] - [ \delta _ { x } - \mathcal { P } ( x ) ] [ ( 1 - \epsilon ) \mathcal { P } ( y, x ) + \delta _ { y, x } ] } { [ ( 1 - \epsilon ) \mathcal { P } ( x ) + \epsilon \delta _ { x } ] ^ { 2 } } \\ & = \frac { \mathcal { P } ( y, x ) \delta _ { x } - \delta _ { y, x } \mathcal { P } ( x ) } { \mathcal { P } ( x ) ^ { 2 } } \text{ (since $\epsilon = 0$)}$$

$$\mathcal { P } ( x ) ^ { 2 }$$

Therefore, The EIF is:

$$\frac { \partial \tilde { \mathcal { P } } _ { \epsilon } } { \partial \epsilon } | _ { \epsilon = 0 } = \int y d \left [ \frac { \mathcal { P } ( y, x ) \delta _ { x } - \delta _ { y, x } \mathcal { P } ( x ) } { \mathcal { P } ( x ) ^ { 2 } } \right ]$$

We get:

$$& = \frac { \delta _ { x } } { p ( x ) } \left [ \int y d \frac { \delta _ { y, x } } { \delta _ { x } } - \int y d \frac { \mathcal { P } ( y, x ) } { \mathcal { P } ( x ) } \right ] \\ & = \frac { \mathbb { 1 } _ { x } } { p ( x ) } \left [ \int y d \delta _ { y | x } - \int y d \mathcal { P } ( y | x ) \right ] \, ( \text{Since } \frac { \mathcal { P } ( y, x ) } { \mathcal { P } ( x ) } = \mathcal { P } ( y | x ) ; \int d \delta = 1 ) \\ & = \frac { \mathbb { 1 } _ { x } } { p ( x ) } [ y - E _ { \mathcal { P } } [ Y | X = x ] ] \\ \text{We get:} \, & \, \text{UNI} \, \vartheta \, \bigcup$$

(III.25)

$$\left ( \Psi = E _ { \mathcal { P } } [ Y | A = a, X = x ], \mathcal { P }, ( y, a, x ) \right ) = \frac { 1 _ { ( a, x ) } } { p [ A = a, x ] } \left [ y - E _ { \mathcal { P } } [ Y | A = a, X ] \right ] \\ \text{Whish imnific that $A$}, ( r ) \right ) \text{ with $E$mentation}, \text{ if $2$} \text{ } \text{ if $2$} \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{                                                                                                                                                                                                         }$$

Which implies that ϕ µ ( a ( x )). With Equations III.23, III.24 and III.25, we have the EIF for the counterfactual estimand ϕ ψ ( a ) = E E [ X [ Y A | = a, X ]] , P , ( y, a, x )). Since µ a ( x ) = E Y A [ | = a, X = x ], ϕ µ ( a ( x )) = ✶ ( a,x ) p a,x ( ) [ y -µ a ( x ) ] ; p x ( ) = p x ( ), ϕ p X ( ( )) = ✶ x -p x ( ); and π a ( X ) = p A ( = a X | ), ψ a = E Y A [ | = a ]. Therefore, we have:

$$\text{we have} \colon & \\ \phi ( \psi _ { a } ) = \sum _ { x } \left [ \left ( \frac { 1 _ { a, x } } { p ( a, x ) } [ y - \mu _ { a } ( x ) ] p ( x ) \right ) + \left ( \mu _ { a } ( x ) ( 1 _ { x } - p ( x ) ) \right ) \right ] \\ = \frac { 1 ( a ) } { \pi _ { a } ( x ) } [ y - \mu _ { a } ( x ) ] + \mu _ { a } ( x ) - \psi _ { a } \\ = \frac { 1 ( A = a ) } { p ( A = a | X ) } ( Y - E [ Y | A = a, X ] ) + E [ Y | A = a, X ] - E [ Y | A = a ] \\ \text{Thus, the EIF for $\psi_{a}$ with discrete measurable elemtes $\phi(\psi_{a},\mathcal{P}$, (Y_{i},X_{i},A_{i}))$ is} \colon &$$

Thus, the EIF for ψ a with discrete measurable elemtes ϕ ψ , ( a P , ( Y , X , A i i i )) is:

$$( \text{III.26} ) \quad \phi ( \psi = \psi ( \mathcal { P } ), \mathcal { P }, ( Y _ { i }, X _ { i }, A _ { i } ) ) = \frac { A _ { i } } { \pi _ { a } ( X _ { i } ) } [ ( Y _ { i } - \mu _ { a } ( X _ { i } ) ] + \mu _ { a } ( X _ { i } ) - \psi _ { a }$$

In other words, and

$$\phi ( \psi _ { 1 } ) = \frac { \mathbb { 1 } ( A _ { i } = 1 ) } { \pi _ { 1 } ( X _ { i } ) } ( Y _ { i } - \mu _ { 1 } ( X _ { i } ) ) + \mu _ { 1 } ( X _ { i } ) - \psi _ { 1 }$$

$$\phi ( \psi _ { 0 } ) = \frac { \mathbb { 1 } ( A _ { i } = 0 ) } { \pi _ { 0 } ( X _ { i } ) } ( Y _ { i } - \mu _ { 0 } ( X _ { i } ) ) + \mu _ { 0 } ( X _ { i } ) - \psi _ { 0 } \\ \text{d use the central identity for the influence function } ( \text{Equation I}$$

We could use the central identity for the influence function (Equation III.21) to verify if our EIF is correct.

## PROOF:

Suppose we are verifying the EIF for ψ 0 . Using the factorizing of the tangent space technique, we could decompose the gradient of the estimator in any direction into:

$$\nabla _ { s } \psi _ { 0 } = \nabla _ { s _ { Y | A, X } } \psi _ { 0 } + \nabla _ { s _ { A | X } } \psi _ { 0 } + \nabla _ { s _ { X } } \psi _ { 0 }$$

Similarly, for the expectation of the influence and the score function, we could also decompose it through the algebra of expectations:

$$E [ \phi s ] = E [ \phi s _ { Y | A, X } ] + E [ \phi s _ { A | X } ] + E [ \phi s _ { X } ]$$

Therefore, we need to prove that ∇ s Y A,X | ψ 0 = E ϕs [ Y A,X | ], ∇ s A X | ψ 0 = E ϕs [ A X | ], and ∇ s X ψ 0 = E ϕs [ X ], correspondingly. We begin with ∇ s X ψ 0 = E ϕs [ X ].

$$\nabla _ { s _ { X } } \psi _ { 0 } & = \frac { \partial } { \partial \epsilon } \int _ { x } \int _ { y } \tilde { p } _ { \epsilon } ( y | 0, x ) d y \tilde { p } _ { \epsilon } ( x ) d x = \int _ { x } \int _ { y } p ( y | 0, x ) d y \frac { \partial } { \partial \epsilon } ( 1 + \epsilon s ) p ( x ) d x \\ & = \int _ { x } \int _ { y } p ( y | 0, x ) d y s _ { X } p ( x ) d x = \int _ { x } \mu _ { 0 } ( x ) s _ { X } p ( x ) d x$$

And,

$$& \text{and}, \\ & E [ \phi s _ { X } ] = E \left [ \left ( \frac { 1 } { } ( A = 0 ) _ { \pi _ { 0 } ( x ) } [ y - \mu _ { 0 } ( x ) ] + \mu _ { 0 } ( x ) \right ) s _ { X } \right ] \\ & \quad = \int _ { x } \int _ { y } \left [ \left ( \frac { 1 } { } ( A = 0 ) _ { \pi _ { 0 } ( x ) } [ y - \mu _ { 0 } ( x ) ] + \mu _ { 0 } ( x ) \right ) s _ { X } \right ] p ( y | x, a ) \, d y \, p ( x ) \, d x \\ & \quad = \int _ { x } s _ { X } \left [ \int _ { y } \frac { 1 } { } ( A = 0 ) _ { \pi _ { 0 } ( x ) } [ y - \mu _ { 0 } ( x ) ] p ( y | x, a ) \, d y + \mu _ { 0 } ( x ) \int _ { y } p ( y | x, a ) \, d y \right ] p ( x ) \, d x \\ & \text{Notice that } \int _ { y } p ( y | x, a ) d y = 1 \, \text{and} \, \int _ { y } [ y - \mu _ { 0 } ( x ) ] p ( y | x, a ) d y = E [ y - \mu _ { 0 } ( x ) | x, 0 ] = 0,$$

Notice that ∫ y p y x, a dy ( | ) = 1 and ∫ y [ y -µ 0 ( x )] p y x, a dy ( | ) = E y [ -µ 0 ( x ) | x, 0] = 0, thus,

$$E [ \phi s _ { X } ] = \int _ { x } s _ { X } \mu _ { 0 } ( x ) p ( x ) d x = \nabla _ { s _ { X } } \psi _ { 0 }$$

Similarly,

$$\nabla _ { s _ { A | X } } \psi _ { 0 } & = \int _ { x } \int _ { y } p ( y \, | \, 0, x ) d y \frac { \partial } { \partial \epsilon } ( 1 + \epsilon s _ { A | X } ) p ( x ) d ( x ) \\ & = \int _ { x } \mu _ { 0 } s _ { A | X } p ( x ) d x = 0$$

$$E [ \phi s _ { A | X } ] = \int _ { x } \int _ { a } \phi s _ { A | X } ( a | x ) p ( a | x ) p ( x ) d a d x$$

At last,

$$\nabla _ { s _ { Y | A, X } } \psi _ { 0 } & = \int _ { x } \int _ { y } p ( y | x, 0 ) s _ { Y | A, X } d y p ( x ) d x \\ & = \int _ { x } p ( x ) \Big { ( } \int _ { y } p _ { 0 } ( y | x ) s _ { Y | A, X } d y \Big { ) } d x \\ & = \int _ { x } ( y | A = 0 ) s _ { Y | A, X } d x$$

$$E [ \phi s _ { Y | A, X } ] & = \int _ { x } \int _ { y } \left [ \left ( \frac { 1 ( A = 0 ) } { \pi _ { 0 } ( x ) } [ y - \mu _ { 0 } ( x ) ] + \mu _ { 0 } ( x ) \right ) s _ { Y | A, X } \right ] p ( y | x, a ) \, d y \, p ( x ) \, d x \\ & = \int _ { x } \left [ \frac { 1 ( A = 0 ) } { \pi _ { 0 } ( x ) } y s _ { Y | A, X } \right ] p ( x ) d x \\ & = \int _ { x } ( y | A = 0 ) s _ { Y | A, X } p ( x ) d x = \nabla _ { s _ { Y | A, X } } \psi _ { 0 }. \\ T h o w f o w \quad E [ \Delta ] - \nabla _ { a } \, \partial \, \text{and with to on } \text{th to the file} \, \text{www}.$$

Therefore, E ϕs [ ] = ∇ s ψ and our derivation on the EIF is correct.

Further, we consider the non-saturated model scenario. The average treatment effect under the random-controlled trial setting is a non-saturated model since we have placed restrictions on the treatment and control cases. Therefore, unlike the saturated models, we cannot derive the efficient influence function with an arbitrary score function, as the influence function may not be efficient. However, we could start with a known RAL estimator and derive its influence function through the definition of asymptotical linearity (Equation III.2 instead of the Gateaux derivative definition in Equation III.18), and then project it (find its minimized square error) onto the tangent space. As elaborated before, the EIF should be the projection of any influence functions on the tangent space. However, it might be hard to directly find the projection on the tangent space if the underlying estimator is complex. If so, we could still use the factorization technique that first projects the influence function onto tangent subspaces and then sums the factorized influence functions together.

For the ATE under the RCT setting, we start with the IPW estimator. Obviously, as shown in Equation II.3, the IPW estimator is a regular and asymptotical

$$& = \int _ { x } \phi \underbrace { \left ( \int _ { a } s _ { A | X } ( a | x ) p ( a | x ) d a \right ) } _ { = 0 } p ( x ) d x \\ & = 0 = \nabla _ { s _ { A | X } } \psi _ { 0 }$$

linear one. We still use the estimation on the EIF for ψ 0 as an example below:

$$\hat { \psi } _ { 0 } ^ { I P W } = E \left [ \frac { 1 ( A _ { i } = 0 ) } { \pi _ { 0 } ( X _ { i } ) } Y _ { i } \right ] = \sum _ { i = 1 } ^ { n } \left [ \frac { 1 ( A _ { i } = 0 ) } { \pi _ { 0 } ( X _ { i } ) } Y _ { i } \right ]. \\ \intertext { a b o d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d o n d }$$

Therefore, based on Equation III.2, we could derive the influence function as:

$$\phi _ { 0 } ^ { I P W } = \frac { \mathbb { 1 } ( A _ { i } = 0 ) } { \pi _ { 0 } ( X _ { i } ) } Y _ { i } - \psi _ { 0 }. \\ \ n f l u n e c e \ f n { \text{function from the IPW estir} }$$

Then we project the influence function from the IPW estimator into the tangent subspace of T Y A,X | , T A X | , and T X in which s Y A,X | , s A X | , and s X forms (and T Y,A,X = T Y A,X | ⊕ T A X | ⊕T X .

First, we try to project the influence function of the IPW estimator on the tangent space T X . The projection function is defined to find the score function on the tangent space for which its mean square error with the influence function from the IPW is minimal 45 . For any influence function,

$$P r o j _ { \mathcal { T } _ { X } } ( \phi ) = \arg \min _ { h ( X ) \in \mathcal { T } _ { X } } E [ ( \phi - h ( X ) ^ { 2 } ].$$

Let h † ( X ) = arg min h X ( ) ∈T X E ϕ [( -h x ( ) 2 ], therefore, the residual ϕ -h † ( X ) should be orthogonal to functions on T X :

$$E [ ( \phi - h ^ { \dagger } ( X ) ) h ( X ) ] = 0, \, \forall h ( X ) \in \mathcal { T } _ { X }.$$

Since h † ( X ) satisfies the orthogonality condition, it should be the conditional expectation of ϕ given X :

$$h ^ { \dagger } ( X ) = E [ \phi | X ].$$

Therefore, we have the projection of the influence function for the IPW estimator on the tangent space T X as its conditional expectation on the X axis:

$$\ t n e \tan g e n t { \text{space $I_{X}$ as its conjugation expectation on the $A$ axis} } ; \\ \phi _ { 0 } ^ { \dagger } ( \mathcal { T } _ { X } ) & = E [ \phi _ { 0 } ^ { I P W } \, | \, X ] = E \left [ \left ( \frac { 1 ( A _ { i } = 0 ) } { \pi _ { 0 } ( X _ { i } ) } Y _ { i } - \psi _ { 0 } \right ) \, | \, X \right ] \\ & = E \left [ E \left [ \frac { 1 ( A _ { i } = 0 ) } { \pi _ { 0 } ( X _ { i } ) } Y _ { i } \, | \, A _ { i }, X \right ] \, | \, X \right ] - \psi _ { 0 } \\ & = E \left [ E \left [ \frac { Y _ { i } } { \pi _ { 0 } ( X _ { i } ) } \, | \, A _ { i } = 0, X \right ] 1 ( A _ { i } = 0 ) \, | \, X \right ] - \psi _ { 0 } \\ & = E \left [ Y _ { i } | A _ { i } = 0, X \right ] - \psi _ { 0 } = \mu _ { 0 } ( X _ { i } ) - \psi _ { 0 } \\ \quad \ \. \.$$

45 We have a sketch Figure 3 illustrating the projection process for the readers' reference for understanding the algebraic process here.

Now we derive the projection of the influence function to the conditional tangent space. If we project the influence function onto the plane spanned by two vectors X and Y , T X,Y , obviously, the projection should be the conditional expectation of the score function conditioned on the joint distribution of X,Y :

$$h ^ { \dagger } ( X, Y ) = E [ \phi | X, Y ]$$

And if we project an influence function towards the plane spanned by X and Y but conditioned on the line which X = 0, then the projection should be the above result, minus its projection on the X -axis, which is:

$$h ^ { \dagger } ( Y | X = 0 ) = E [ \phi | X, Y ] - E [ \phi | X ]$$

We denote the tangent space in this situation T Y,X 0 . Now we discuss the projections of the influence function of the IPW estimator onto the tangent spaces T A 0 | X and T Y A ,X | 0 . Obviously,

$$\phi _ { 0 \langle \mathcal { T } _ { A _ { 0 } | X } \rangle } ^ { \dagger } = 0,$$

as the tangent space T A 0 | X is orthogonal to ϕ IPW 0 and therefore all projections have no length. For the projection of ϕ IPW 0 onto the tangent space T Y A ,X | 0 , we have:

$$\hat { \phi } _ { 0 ( \mathcal { T } _ { Y | A _ { 0 }, X } ) } ^ { \dagger } & = E [ \hat { \phi } _ { 0 } ^ { I P W } | Y, A _ { 0 }, X ] - E [ \hat { \phi } _ { 0 } ^ { I P W } | A _ { 0 }, X ] \\ & = E \left [ \frac { 1 ( A _ { i } ) = 0 } { \pi _ { 0 } ( X _ { i } ) } Y _ { i } - \psi _ { 0 } \, | \, Y, A _ { 0 }, X \right ] - E \left [ \frac { 1 ( A _ { i } = 0 ) } { \pi _ { 0 } ( X _ { i } ) } Y _ { i } - \psi _ { 0 } \, | \, \ A _ { 0 }, X \right ] \\ & = \left ( \frac { 1 ( A _ { i } = 0 ) } { \pi _ { 0 } ( X _ { i } ) } Y _ { i } - \psi _ { 0 } \right ) - \left ( \frac { 1 ( A _ { i } ) = 0 } { \pi _ { 0 } ( X _ { i } ) } \mu _ { 0 } ( X _ { i } ) - \psi _ { 0 } \right ) \\ \text{Therefore, we sum up the three sub-EIFs on the three tangent subspaces and get} \\ \text{$1$} \sim \text{E} \, \frac { 1 } { 2 } \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \sim \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \text{$1$} \quad \end{$$

Therefore, we sum up the three sub-EIFs on the three tangent subspaces and get the EIF for ψ 0 under the RCT settings:

$$& \text{the EIF for $\psi_{0}$ under the RCT settings} ; \\ & \phi ^ { \dagger } ( \tilde { \psi } _ { 0 } ^ { I P W } ) \\ & \quad = \phi _ { 0 } ^ { \dagger } ( \mathcal { T } _ { X } ) + \phi _ { 0 } ^ { \dagger } ( \mathcal { T } _ { X _ { 0 } | X } ) + \phi _ { 0 } ^ { \dagger } ( \mathcal { T } _ { Y | A _ { 0 }, X } ) \\ & \quad = ( \mu _ { 0 } ( X _ { i } ) - \psi _ { 0 } ) + 0 + \left ( \frac { 1 } { \pi _ { 0 } ( X _ { i } ) } Y _ { i } - \psi _ { 0 } \right ) - \left ( \frac { 1 ( A _ { i } ) = 0 } { \pi _ { 0 } ( X _ { i } ) } \mu _ { 0 } ( X _ { i } ) - \psi _ { 0 } \right ) \\ & \quad = \frac { 1 ( A _ { i } = 0 ) } { \pi _ { 0 } ( X _ { i } ) } ( Y _ { i } - \mu _ { 0 } ( X _ { i } ) ) + \mu _ { 0 } ( X _ { i } ) - \psi _ { 0 } \\ & \text{Which is, unsurprisingly, exactly the EIF for $\psi_{0}$ we obtained from the saturated} \\ & \text{model. Similarly, the EIF for $nsi$}. from the non-saturated model of the RCT will$$

Which is, unsurprisingly, exactly the EIF for ψ 0 we obtained from the saturated model. Similarly, the EIF for psi 1 from the non-saturated model of the RCT will also be the same as the EIF from the saturated model of the observational study.

With the EIFs for ψ 1 and ψ 0 , we can get the EIF for the ATE:

$$\phi ( \psi ) & = \phi ( \psi _ { 1 } ) - \phi ( \psi _ { 0 } ) \\ & = \left [ \frac { 1 ( A _ { i } = 1 ) } { \pi _ { 1 } ( X _ { i } ) } ( Y _ { i } - \mu _ { 1 } ( X _ { i } ) ) + \mu _ { 1 } ( X _ { i } ) - \psi _ { 1 } \right ] - \left [ \frac { 1 ( A _ { i } = 0 ) } { \pi _ { 0 } ( X _ { i } ) } ( Y _ { i } - \mu _ { 0 } ( X _ { i } ) ) + \mu _ { 0 } ( X _ { i } ) - \psi _ { 0 } \right ] \\ ( \text{III.27} ) & \quad \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r } \ \mathfrak { r }$$

(III.27)

$$\stackrel { \searrow \dots \cdot \prime } { = \frac { 1 } { } ( A _ { i } = 1 ) } { \pi ( X _ { i } ) } \left ( Y _ { i } - \mu _ { 1 } ( X _ { i } ) \right ) - \frac { 1 - 1 ( A _ { i } = 1 ) } { 1 - \pi ( X _ { i } ) } \left ( Y _ { i } - \mu _ { 0 } ( X _ { i } ) \right ) + \left ( \mu _ { 1 } ( X _ { i } ) - \mu _ { 0 } ( X _ { i } ) \right ) - \psi \\ \text{As we let $\pi( X:\colon)=\pi_{i}(\X:\colon)=\mathcal{P}(A:\colon=1\mathcal{L}X:\colon)$ and $\psi_{i}-\psi_{n}=\psi_{i}$. With Foundation III 2}$$

As we let π X ( i ) = π 1 ( X i ) = P ( A i = 1 | X i ) and ψ 1 -ψ 0 = ψ . With Equation III.2 at the start of this section, we could derive the efficient estimator for the average treatment effect:

(III.28)

$$& \hat { \psi } = \frac { 1 } { n } \sum _ { i = 1 } ^ { n } \left [ \frac { 1 ( A _ { i } = 1 ) } { \hat { \pi } ( X _ { i } ) } \left ( Y _ { i } - \hat { \mu } _ { 1 } ( X _ { i } ) \right ) - \frac { 1 - 1 ( A _ { i } = 1 ) } { 1 - \hat { \pi } ( X _ { i } ) } \left ( Y _ { i } - \hat { \mu } _ { 0 } ( X _ { i } ) \right ) + ( \hat { \mu } _ { 1 } ( X _ { i } ) - \hat { \mu } _ { 0 } ( X _ { i } ) \right ) \right ] \\ & \mathfrak { P } \dots \, \mathfrak { P } \, \infty \, \colon \, \colon \, \colon \, \colon \, \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \end{array}$$

Equation III.28 is the core equation in the whole thesis. We call the estimator in Equation III.28 in several ways. It is, as we elaborated, the efficient causal estimator, and it is also the doubly robust (DR) causal estimator since the estimator will be consistent as either our specification of π X ( i ) or our specification of µ a ( X i ) is correct. Further, it is also called the debiased machine learning (DML) estimator as it eliminates the bias between the estimator and the true estimand with the double robustness we will give the details below. Meanwhile, as we mentioned earlier, since cross-validation will be used in the algorithm generating the estimator to control the convergence speed of the empirical process, the estimator is also called double machine learning (DML) estimator. Finally, some literature also calls this estimator the Neyman-orthogonal estimator of the ATE, as the EIF satisfies the Neyman orthogonality, yielding this result.

From the perspective of debiased estimation, when the propensity score function π X ( i ) is correctly specified, we have:

$$\mathbb { E } \left [ \frac { A _ { i } } { \pi ( X _ { i } ) } ( Y _ { i } - \mu _ { 1 } ( X _ { i } ) ) \ | \ X _ { i } \right ] = 0$$

and

$$\mathbb { E } \left [ \frac { 1 - A _ { i } } { 1 - \pi ( X _ { i } ) } ( Y _ { i } - \mu _ { 0 } ( X _ { i } ) ) \ | \ X _ { i } \right ] = 0.$$

Therefore, the first two terms of the influence function become mean-zero con- ditional on X , leaving:

$$\hat { \psi } _ { \text{DR} } = \frac { 1 } { n } \sum _ { i = 1 } ^ { n } \left [ \mu _ { 1 } ( X _ { i } ) - \mu _ { 0 } ( X _ { i } ) \right ].$$

So even if µ 1 ( X i ) and µ 0 ( X i ) are misspecified, the terms involving the propensity score correct the bias introduced by the misspecified outcome models, resulting in a consistent estimator for the ATE.

Similarly, if the outcome regression models µ 1 ( X ) and µ 0 ( X ) are correctly specified, we have:

$$\mathbb { E } [ Y _ { i } \ | \ A _ { i } = 1, X _ { i } ] = \mu _ { 1 } ( X _ { i } ) \ \text{ and } \ \mathbb { E } [ Y _ { i } \ | \ A _ { i } = 0, X _ { i } ] = \mu _ { 0 } ( X _ { i } ).$$

In this case, the terms Y i -µ 1 ( X i ) and Y i -µ 0 ( X i ) are mean-zero conditional on A i and X i . Thus, the first two terms of the influence function average out to zero:

$$\hat { \psi } _ { \text{DR} } = \frac { 1 } { n } \sum _ { i = 1 } ^ { n } \left [ \frac { A _ { i } } { \pi ( X _ { i } ) } ( Y _ { i } - \mu _ { 1 } ( X _ { i } ) ) - \frac { 1 - A _ { i } } { 1 - \pi ( X _ { i } ) } ( Y _ { i } - \mu _ { 0 } ( X _ { i } ) ) + \mu _ { 1 } ( X _ { i } ) - \mu _ { 0 } ( X _ { i } ) \right ].$$

Here, even if π X ( i ) is misspecified, the correctly specified outcome regression models ensure that the estimator is consistent for the ATE.

From the perspective of the convergence rate, we will find that the estimator also satisfies double robustness with its higher-order Distributional Taylor Expansion. According to Equation III.15, We can derive the second-order remainder of the estimator ψ a :

$$\text{coupling } \omega \colon \text{$\Gamma$} _ { a } \colon \\ R _ { 2 } & = \psi ( \tilde { \Gamma } _ { \epsilon } ) - \psi ( P ) + \frac { \partial } { \partial \epsilon } \psi ( \tilde { \Gamma } _ { \epsilon } ) | _ { \epsilon = 1 } \\ & = \psi ( \tilde { \Gamma } _ { \epsilon } ) - \psi ( P ) + \frac { 1 } { n } \sum _ { i = 1 } ^ { n } \phi ( \psi ( \tilde { \Gamma } _ { \epsilon } ) ) \\ & = \tilde { \psi } _ { a } - E [ \mu _ { a } ( X _ { i } ) ] + E \left [ \frac { 1 ( A _ { i } = a ) } { \tilde { \pi } _ { a } ( X _ { i } ) } ( Y _ { i } - \tilde { \mu } _ { a } ( X _ { i } ) ) + \tilde { \mu } _ { a } ( X _ { i } ) \right ] - E \left [ \tilde { \psi } ( a ) \right ] \\ & = E \left [ \frac { \pi _ { a } ( X _ { i } ) } { \tilde { \pi } _ { a } ( X _ { i } ) } \mu _ { a } ( X _ { i } ) - \mu _ { a } ( X _ { i } ) \right ] + E \left [ \tilde { \mu } _ { a } ( X _ { i } ) - \frac { \pi _ { a } ( X _ { i } ) } { \tilde { \pi } _ { a } ( X _ { i } ) } \tilde { \mu } _ { a } ( X _ { i } ) \right ] \\ & = E \left [ \frac { 1 } { \tilde { \pi } _ { a } ( X _ { i } ) } \left ( \pi _ { a } ( X _ { i } ) - \tilde { \nu } _ { a } ( X _ { i } ) \right ) ( \mu _ { a } ( X _ { i } ) - \tilde { \mu } _ { a } ( X ) ) \right ] \\ & \leq \| \pi _ { a } ( X _ { i } ) - \tilde { \nu } _ { a } ( X _ { i } ) \| | \| \mu _ { a } ( X _ { i } ) - \tilde { \mu } _ { a } ( X _ { i } ) \| \left ( \text{Cauchy-Schwarz Inequality} \right ) \\ \text{We could see from this formation that either } \pi _ { a } ( X _ { i } ) = \tilde { \nu } _ { a } ( X _ { i } ) \text{ or } \mu _ { a } ( X _ { i } ) - \tilde { \mu } _ { a } ( X _ { i } ),$$

We could see from this formation that either π a ( X i ) = ˆ π a ( X i ) or µ a ( X i ) -ˆ ( µ a X i ), the second-order remainder turns to zero. Therefore, either ˆ π a ( X i ) or ˆ ( µ a X i ) is correctly specified, our estimator is unbiased.

Further, this inequality also suggests the convergence rate of the doubly robust estimator. Since we require the second-order remainder to converge at the rate of o p ( n -1 2 / ), we could require both || π a ( X i ) -ˆ ( π a X i ) || and || µ a ( X i ) -ˆ ( µ a X i ) || to converge at the rate of o p ( n -1 4 / ). Regularized machine learning methods and cross-validation need to be used for both the propensity score model ˆ π a ( X i ) and the outcome model ˆ ( µ a X i ). For instance, we could choose gradient boosting machines or random forests to predict the propensity and generalized additive models (GAM) or random forests to predict the outcome 46 , and then use regularized techniques like lasso, ridge, or elastic nets to control the complexity of the models.

So, practically, we have the algorithm deriving the doubly robust (efficient) causal estimator:

## 1) Set up Input:

- · Dataset { ( X , A , Y i i i ) } n i =1 , where X i represents covariates, A i represents treatment assignment (0 or 1), and Y i represents outcomes.
- · Number of folds for cross-validation k .

## 2) Split Dataset for Cross-Validation

- · Randomly split the dataset into k approximately equal-sized folds. Each fold will be used as a validation set while the remaining k -1 folds will be used for training.
- · Label these folds as { D ,D ,..., D 1 2 k } .

## 3) Cross-Validation Loop

- · Initialize lists to store the fold-specific estimates of ˆ ψ 1 and ˆ . ψ 0
- · For each fold j (from 1 to k ):
- - Training Set: Combine all folds except D j to create the training set { ( X , A , Y i i i ) } i ∈ Training Set .
- - Validation Set: Use fold D j as the validation set { ( X , A , Y i i i ) } i ∈ Validation Set .

## Estimate Propensity Scores in the Training Set

- · Fit a propensity score model ˆ( π X ) using the training set.
- · Calculate the estimated propensity scores ˆ( π X i ) for all i in the validation set.

46 A pretty useful programming package for model choice is called super-learner(Polley et al. 2023), which could be convenient for social scientists choose the appropriate machine learning models for model fitting.

## 4)

## 5) Estimate Outcome Regressions in the Training Set

- · Fit outcome regression models ˆ ( µ 0 X ) and ˆ ( µ 1 X ) using the training set.
- · Calculate the predicted outcomes ˆ ( µ 0 X i ) and ˆ ( µ 1 X i ) for all i in the validation set.

## 6) Calculate the Doubly Robust Estimator in the Validation Set

- · Initialize two variables to accumulate the contributions from the treated and control groups in the validation set: ˆ ψ j 1 and ˆ . ψ j 0
- · For each observation i in the validation set:
- -Compute the contribution for the treated group:

$$\hat { \psi } _ { 1 i } = \frac { A _ { i } } { \hat { \pi } ( X _ { i } ) } \left ( Y _ { i } - \hat { \mu } _ { 1 } ( X _ { i } ) \right ) + \hat { \mu } _ { 1 } ( X _ { i } )$$

- -Compute the contribution for the control group:

$$\hat { \psi } _ { 0 i } = \frac { 1 - A _ { i } } { 1 - \hat { \pi } ( X _ { i } ) } \left ( Y _ { i } - \hat { \mu } _ { 0 } ( X _ { i } ) \right ) + \hat { \mu } _ { 0 } ( X _ { i } )$$

- -Accumulate the contributions:

$$\hat { \psi } _ { 1 } ^ { j } \leftarrow \hat { \psi } _ { 1 } ^ { j } + \hat { \psi } _ { 1 i }$$

$$\hat { \psi } _ { 1 } ^ { j } & \leftarrow \hat { \psi } _ { 1 } ^ { j } + \hat { \psi } _ { 1 i } \\ \hat { \psi } _ { 0 } ^ { j } & \leftarrow \hat { \psi } _ { 0 } ^ { j } + \hat { \psi } _ { 0 i } \\. & \dots.$$

- · Calculate the averages for the validation set:

$$\hat { \psi } _ { 1 } ^ { j } = \frac { 1 } { n _ { j } } \sum _ { i \in \text{Validation Set} } \hat { \psi } _ { 1 i }$$

$$\hat { \psi } _ { 0 } ^ { j } = \frac { 1 } { n _ { j } } \sum _ { i \in \text{Validation Set} } \hat { \psi } _ { 0 i }$$

- · Store the fold-specific estimates.

## 7) Aggregate Results Across Folds

- · Calculate the overall estimates by averaging the fold-specific estimates:

$$\hat { \psi } _ { 1 } = \frac { 1 } { k } \sum _ { j = 1 } ^ { k } \hat { \psi } _ { 1 } ^ { j }$$

$$\hat { \psi } _ { 0 } = \frac { 1 } { k } \sum _ { j = 1 } ^ { k } \hat { \psi } _ { 0 } ^ { j }$$

## 8) Compute the Average Treatment Effect (ATE)

- · Estimate the ATE:

$$\hat { \psi } = \hat { \psi } _ { 1 } - \hat { \psi } _ { 0 }$$

## 9) Output

- · The doubly robust estimator of the average treatment effect ˆ . ψ

## E. Further Discussions

So far, we have introduced methods yielding efficient influence functions and estimators and given the detailed algebraic transformation for the efficient/doubly robust causal estimator. However, in some cases, the underlying estimand is not that precise, and the decomposition of its tangent space is not feasible; also, in some cases, we only have a very small sample size, and the information is limited. If so, we can only use the empirical cases to derive the heuristic approximation for its efficient estimator (which is close but not attained, the Cramer-Rao bound). We will give an example in the next chapter when we need to derive the efficient/doubly robust estimation for the left-truncated-right-censored survival data.

## IV. Conclusion

In this chapter, we introduce the basic ideas and mathematical tools to perform the causal inference in an efficient/doubly robust way from the observational/ survey data in social science. Intrinsically, to correctly identify the causal estimand from the observational data (via the statistical estimand), for all the methods, no matter whether it is a parametric (regressive) model, a nonparametric (machine learning) model, or a semi-parametric (bias reduction) model, the key is the same: to correctly specify/identify the two models: the propensity score model which allocates cases into the treatment and control group and the outcome result model which specifies the factual or counterfactual outcomes.

The doubly robust estimator provides a toolbox that, compared to the previous models, requires the correct specification of both the propensity score model and the outcome result model; we only need to correctly specify one of them. This makes social science research, especially under theoretical-driven studies, much more convenient to yield an unbiased and efficient estimator for the treatment effect. In empirical studies, researchers are always afraid of omitted variable bias when specifying the outcome model, but if the researchers' theories and hypothesis could make the propensities allocating to the treatment and control groups deterministic, then the results are robust (this is pretty like the idea behind the local average treatment effect estimation).

In the chapters below, we will again use the asymptotical analysis methods (with the score and influence function) to yield efficient estimators under different data structures and model settings. In this chapter, we just give readers from social science backgrounds a preliminary introduction (with necessary mathematical transforms) to this area. The method for semiparametric doubly robust target double machine learning (Kennedy 2022) is definitely one of the fastest developing areas in statistics, econometrics, data science, and relevant discipline's methodological discussions. Like other disciplines, social scientists and demographers need the toolbox to have better causal estimations of their research interests.

## REFERENCES

Angrist, J. D. (1990). Lifetime earnings and the vietnam era draft lottery: Evidence from social security administrative records. American Economic Review 80 (3), 313-336.

- Angrist, J. D. and A. B. Krueger (2001). Instrumental variables and the search for identification: From supply and demand to natural experiments. Journal of Economic Perspectives 15 (4), 69-85.
- Angrist, J. D. and V. Lavy (1999). Using maimonides' rule to estimate the effect of class size on scholastic achievement. Quarterly Journal of Economics 114 (2), 533-575.

Autor, D. H., D. Dorn, and G. H. Hanson (2013). The china syndrome: Local labor market effects of import competition in the united states. American Economic Review 103 (6), 2121-2168.

- Baker, S. G., B. S. Kramer, and K. S. Lindeman (2016). Latent class instrumental variables: a clinical and biostatistical perspective. Statistics in medicine 35 (1), 147-160.

Card, D. (1999). The causal effect of education on earnings. In Handbook of Labor Economics , Volume 3, pp. 1801-1863. Elsevier.

- Chernozhukov, V., D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey, and J. Robins (2018). Double/debiased/neyman machine learning of treatment effects. American Economic Review 108 (5), 899-927.
- Chernozhukov, V., M. Demirer, E. Duflo, and I. Fernandez-Val (2018). Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal 21 (1), C1-C68.

Dell, M., B. F. Jones, and B. A. Olken (2009). Temperature and income: Reconciling new cross-sectional and panel estimates. American Economic Review 99 (2), 198-204.

| Diamond, A. and J. S. Sekhon (2013). Genetic matching for estimating causal effects: A general multivariate matching method for achieving balance in ob- servational studies. Review of Economics and Statistics 95 (3), 932-945.   |
|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| Fisher, A. and E. H. Kennedy (2021). Visually communicating and teaching intuition for influence functions. The American Statistician 75 (2), 162-172.                                                                              |
| Fletcher, J. M. and B. L. Wolfe (2009). The effects of teenage childbearing on the short-and long-term health behaviors of mothers. Journal of Population Economics 22 (3), 575-597.                                                |
| Gelman, A. and J. Hill (2006). Data Analysis Using Regression and Multi- level/Hierarchical Models . Cambridge: Cambridge University Press.                                                                                         |
| Haavelmo, T. (1944). The probability approach in econometrics. Economet- rica 12 (Supplement), iii-115.                                                                                                                             |
| Heckman, J., H. Ichimura, J. Smith, and P. Todd (1998). Characterizing selection bias using experimental data. Econometrica 66 (5), 1017-1098.                                                                                      |
| Holland, P. W. (1986). Statistics and causal inference. Journal of the American Statistical Association 81 (396), 945-960.                                                                                                          |
| Imbens, G. W. and J. D. Angrist (1994). Identification and estimation of local average treatment effects. Econometrica 62 (2), 467-475.                                                                                             |
| Kennedy, E. H. (2022). Semiparametric doubly robust targeted double machine learning: A review. arXiv preprint arXiv:2203.06469 . Last revised on January 26, 2023.                                                                 |
| King, G. and R. Nielsen (2019). Why propensity scores should not be used for matching. Political Analysis 27 (4), 435-454.                                                                                                          |
| Neyman, J., K. Iwaszkiewicz, and S. Kolodziejczyk (1935). Statistical problems in agricultural experimentation. Supplement to the Journal of the Royal Statistical Society 2 (2), 107-180.                                          |
| Oreopoulos, P. (2006). Estimating average and local average treatment effects of education when compulsory schooling laws really matter. American Economic Review 96 (1), 152-175.                                                  |
| Pearl, J. (2009). Causality . Cambridge University Press.                                                                                                                                                                           |
| Pearl, J. (2015). Trygve haavelmo and the emergence of causal calculus. Econo- metric Theory 31 (1), 152-179.                                                                                                                       |
| Polley, E., A. E. Hubbard, and M. J. van der Laan (2023). SuperLearner: Super Learner Prediction . R package version 2.0-31.                                                                                                        |

Rosenbaum, P. R. (1989). Optimal matching for observational studies. Journal of the American Statistical Association 84 (408), 1024-1032.

Rosenbaum, P. R. and D. B. Rubin (1983). The central role of the propensity score in observational studies for causal effects. Biometrika 70 (1), 41-55.

Rosenbaum, P. R. and D. B. Rubin (1984). Reducing bias in observational studies using subclassification on the propensity score. Journal of the American Statistical Association 79 (387), 516-524.

Roth, J., P. H. Sant'Anna, A. Bilinski, and J. Poe (2023). What's trending in difference-in-differences? a synthesis of the recent econometrics literature. Journal of Econometrics 235 (2), 2218-2244.

Rubin, D. B. (1980). Bias reduction using mahalanobis-metric matching. Biometrics 36 (2), 293-298.

Rubin, D. B. (1990). Comment: Neyman (1923) and causal inference in experiments and observational studies. Statistical Science 5 (4), 472-480.

Stock, J. H. and F. Trebbi (2003). Who invented instrumental variable regression? Journal of Economic Perspectives 17 (3), 177-194.

Stock, J. H., J. H. Wright, and M. Yogo (2002). A survey of weak instruments and weak identification in generalized method of moments. Journal of Business &amp; Economic Statistics 20 (4), 518-529.

van der Vaart, A. W. (1998). Asymptotic Statistics . Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge: Cambridge University Press.

VanderWeele, T. J. (2015). Explanation in causal inference: methods for mediation and interaction . Oxford University Press.

Figure 1. Illustrations on propensity score matching

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00032v1-figure-d68b6abd-cbf9-40e3-8980-440f6e0c4a88.jpg)

Note: The upper panel shows the ideal scenario for propensity score matching, where for each individual in the treatment group, we could find a corresponding individual in the control group with the same propensity value and achieve one-on-one matching. However, as the lower two panels show, with observational data to train and predict the propensity function, the distributions of the treatment and control individuals are unequal along the propensity score, making the one-on-one ideal matching infeasible. But if for the treatment group and for the control group, the outcome is irrelevant to the propensity scores (the middle panel), using any method to get ˆ π a yields the same unbiased estimation for π a . Otherwise, as the lowest panel indicates, the choice of ˆ π a yields bias.

Figure 2. Illustration on Distributional Taylor Expansion

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00032v1-figure-2a85060c-15b6-4198-97cf-fa7f3f27d15b.jpg)

Note: This is a simplified illustration of the Distributional Taylor expansion. Indeed, ˜ and P P are two measures, but we just simplify them as two values. If so, the direction of the score function should be horizontal, and the direction of the gradient of the score function on the estimator should be the tangent line of ψ at ψ ( P ), which has the same direction as the influence function. Thus, the angle of score and influence function points to the same direction as the gradient.

Figure 3. Illustrations on projections of influence functions

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00032v1-figure-feb93e4d-dcf0-4614-ad81-c7f2bfef0101.jpg)

Note: The yellow plane denotes the tangent space. The figure shows projecting the (inefficient) influence function (not on the tangent space) towards the joint distribution of ( X,Y ) and the marginal distribution of ( Y X | = 0).
Guanghui Pan
Methodological Foundations of Modern Causal Inference in Social Science Research
This paper serves as a literature review of methodology concerning the
(modern) causal inference methods to address the causal estimand with
observational/survey data that have been or will be used in social science
research. Mainly, this paper is divided into two parts: inference from
statistical estimand for the causal estimand, in which we reviewed the
assumptions for causal identification and the methodological strategies
addressing the problems if some of the assumptions are violated. We also
discuss the asymptotical analysis concerning the measure from the observational
data to the theoretical measure and replicate the deduction of the
efficient/doubly robust average treatment effect estimator, which is commonly
used in current social science analysis.