<a href="https://colab.research.google.com/github/migariane/DeltaMethodInfluenceFunction/blob/main/DeltaMethod.ipynb" target="_parent"></a>

# Delta Method

Code for reproducing all the results in the paper _The Delta-method and influence function in epidemiology: a reproducible tutorial_


```python
%pip install --upgrade --quiet --force-reinstall "sympy>=1.10.0"
```


```python
#We'll use the sympy library for symbolic computation of variances and expected values
import sympy
from sympy import * 
from sympy.stats import Variance, Expectation, Covariance, Binomial
from sympy.stats.rv import RandomSymbol
sympy.__version__   #version of sympy should be >=1.10
```




    '1.10'



## Introduction to Sympy
[Sympy](https://docs.sympy.org/latest/guides/getting_started/install.html) is a computer algebra system for symbolic calculation. In order to use it you need to declare variables using `Symbol`:


```python
x = Symbol("x")
y = Symbol("y")
```

You can use the `Symbol`s to generate algebraic expressions such as:


```python
squarexy = pow((x + y),2) #(x + y)^2
```

which can be expanded, derived and integrated:


```python
expand(squarexy)
```




$\displaystyle x^{2} + 2 x y + y^{2}$




```python
diff(squarexy,x)
```




$\displaystyle 2 x + 2 y$




```python
integrate(squarexy, x)
```




$\displaystyle \frac{x^{3}}{3} + x^{2} y + x y^{2}$



Random variables can be declared with `RandomSymbol`:


```python
#Declaration of random variables can be made via RandomSymbol and their 
#variance and expected values can be estimated and simplified
X                = RandomSymbol('X')

#Properties of expectation can be computed
expectation_of_x = Expectation(3*X + 2)
expectation_of_x.expand()
```




$\displaystyle 2 + 3 \operatorname{E}\left[X\right]$




```python
#As well as properties of variance
variance_of_x    = Variance(3*X + 2)
variance_of_x.expand()
```




$\displaystyle 9 \operatorname{Var}\left(X\right)$



## Delta-method for the mean

First we declare the constants: $\mu$, $\sigma^2$ and $n$ corresponding to the mean, variance and sample size. We also declare the random variable $\bar{X}$ corresponding to the sample mean which has mean $\mu$ and variance $\sigma^2/n$ 


```python
#Declaration of constants (Symbol) and random variables (RandomSymbol)
mu       = Symbol('mu')
sigmasq  = Symbol('sigma^2', positive = True) #Variance of x
n        = Symbol('n', positive = True, integer = True)
Xbar     = RandomSymbol('Xbar') 
```

We write the function $\phi$ as a function of $\mu$ _i.e._ $\phi(\mu) =  \mu$:


```python
#Declaring phi as a function of mu
def phi(mu):
    return mu
```

We calculate the derivative: $\frac{\partial \phi}{\partial \mu}$:




```python
#Obtaining the derivative
classical_derivative = derive_by_array(phi(mu), mu)
```

We then calculate the direction vector for the Hadamard derivative $v = \bar{X} - \mu$:


```python
v = Xbar - mu
```

The Hadamard derivative in the direction of $v$ (that is $ \partial_{\hat{\theta} - \theta}$) is thus:


```python
hadamard = classical_derivative*v
```

And the approximation is obtained by computing the right side of:
$$
\phi(\hat{\theta}) = \phi(\theta)  + \partial_{\hat{\theta} - \theta} \phi(\theta)
$$


```python
phi(mu) + hadamard
```




$\displaystyle \bar{X}$



Hence the Delta-method's approximation to the variance is:


```python
delta_variance = Variance(phi(mu) + hadamard).expand()
```

and we know from classical inference that $\text{Var}(\bar{X}) = \frac{\sigma^2}{n}$


```python
delta_variance.subs(Variance(Xbar), sigmasq/n)
```




$\displaystyle \frac{\sigma^{2}}{n}$



## Delta-method for the variance of the ratio of two sample means

Let $\{X_1, \dots, X_n\}$ denote a random sample of variables $X$ with mean $\mu_X$ and variance $\sigma^2_X$ and  $\{Y_1, \dots, Y_n\}$ denote a random sample of variables $Y$ with mean $\mu_Y$ and variance $\sigma^2_Y$ we are interested in approximating the variance of 

$$
\phi(\mu_X,\mu_Y) = \frac{\mu_X}{\mu_Y}
$$

to so we declare the symbols $\bar{X},\bar{Y},\mu_X,\mu_Y$ and the function $\phi$:


```python
#Declaration of variables
mu_x       = Symbol('mu_x')
mu_y       = Symbol('mu_y')
n          = Symbol('n', positive=True, integer = True) #Sample size of X and Y
sigmasq_x  = Symbol('sigma^2_X', positive = True) 
sigmasq_y  = Symbol('sigma^2_Y', positive = True) 
sigma_xy   = Symbol('sigma_xy') #Covariance
Xbar       = RandomSymbol('Xbar')
Ybar       = RandomSymbol('Ybar')  

#Declaration of function
def phi(mu_x, mu_y):
  return mu_x/mu_y
```

We calculate the derivative (in this case, gradient) in the direction of 
$$
v = \begin{pmatrix}
    \bar{X} - \mu_X\\
    \bar{Y} - \mu_Y
    \end{pmatrix}
$$


```python
#Obtain the direction vector
v = Matrix([Xbar - mu_x, Ybar - mu_y])

#And calculate the gradient
grad = derive_by_array(phi(mu_x,mu_y), [mu_x, mu_y])
grad = Matrix(grad)

#Thus calculating the Hadamard (directional) derivative:
hadamard = grad.dot(v)
hadamard
```




$\displaystyle - \frac{\mu_{x} \left(- \mu_{y} + \bar{Y}\right)}{\mu_{y}^{2}} + \frac{- \mu_{x} + \bar{X}}{\mu_{y}}$



The variance of the hadamard derivative (influence function) is estimated as follows:


```python
variance_expanded = Variance(hadamard).expand().simplify()
variance_expanded
```




$\displaystyle \frac{\mu_{x}^{2} \operatorname{Var}\left(\bar{Y}\right) - 2 \mu_{x} \mu_{y} \operatorname{Cov}\left(\bar{X}, \bar{Y}\right) + \mu_{y}^{2} \operatorname{Var}\left(\bar{X}\right)}{\mu_{y}^{4}}$



We can further simplify by substituting the variances of $\bar{X}$ and $\bar{Y} as well as the covariances$:


```python
variance_expanded = variance_expanded.subs(Covariance(Xbar,Ybar), sigma_xy/n)
variance_expanded = variance_expanded.subs(Variance(Xbar), sigmasq_x/n)
variance_expanded = variance_expanded.subs(Variance(Ybar), sigmasq_y/n)

variance_expanded.simplify()
```




$\displaystyle \frac{\mu_{x}^{2} \sigma^{2}_{Y} - 2 \mu_{x} \mu_{y} \sigma_{xy} + \mu_{y}^{2} \sigma^{2}_{X}}{\mu_{y}^{4} n}$



## Delta-method for the risk ratio
Consider the risk ratio for two probabilities defined by:
$$
RR(p_1,p_2) = \frac{p_1}{p_2}
$$
and the transformation $\phi$ given by:
$$
\phi(p_1,p_2) = \log\big(RR(p_1,p_2)\big)
$$
In this case, a random sample of $N_1$ elements of one category and $N_2$ elements of a second (independent) category with respective probabilities $p_1$ and $p_2$ result in estimators $\hat{p}_1$ and $\hat{p}_2$ which we declare:


```python
#Declaration of variables
p_1       = Symbol('p_1')
p_2       = Symbol('p_2')
p_1hat    = RandomSymbol('phat_1')
p_2hat    = RandomSymbol('phat_2')  
N_1       = Symbol('N_1', positive=True, integer = True) #Sample size asssociated to p1
N_2       = Symbol('N_2', positive=True, integer = True) #Sample size asssociated to p2

#Declaration of relative risk function
def RR(p_1, p_2):
  return p_1/p_2

#Declaration of phi for delta method
def phi(p_1, p_2):
  return log(RR(p_1, p_2));
```

The derivative is given by:


```python
#Obtain the direction vector
v = Matrix([p_1hat - p_1, p_2hat - p_2])

#And calculate the gradient
grad = derive_by_array(phi(p_1,p_2), [p_1, p_2])
grad = Matrix(grad)

#Thus calculating the Hadamard (directional) derivative:
hadamard = grad.dot(v)
hadamard
```




$\displaystyle - \frac{- p_{2} + \hat{p}_{2}}{p_{2}} + \frac{- p_{1} + \hat{p}_{1}}{p_{1}}$



And the variance where we use that $\text{Var}(\hat{p}_i) = p_i(1-p_i)/N_i$ and that the covariance is zero due to independence


```python
variance_expanded = Variance(hadamard).expand().simplify()
variance_expanded
```




$\displaystyle \frac{\operatorname{Var}\left(\hat{p}_{2}\right)}{p_{2}^{2}} - \frac{2 \operatorname{Cov}\left(\hat{p}_{1}, \hat{p}_{2}\right)}{p_{1} p_{2}} + \frac{\operatorname{Var}\left(\hat{p}_{1}\right)}{p_{1}^{2}}$



Which can be simplified into:


```python
variance_expanded = variance_expanded.subs(Covariance(p_1hat,p_2hat), 0.0)
variance_expanded = variance_expanded.subs(Variance(p_1hat), p_1*(1 - p_1)/N_1)
variance_expanded = variance_expanded.subs(Variance(p_2hat), p_2*(1 - p_2)/N_2)

variance_expanded
```




$\displaystyle \frac{1 - p_{2}}{N_{2} p_{2}} + \frac{1 - p_{1}}{N_{1} p_{1}}$



## Counterexample

Consider the attributable fraction among the exposed for an exposure level $x$ given by:
$$
\textrm{AF}_e(x) = \dfrac{RR(\theta,x)-1}{RR(\theta,x)}.
$$
we compute the expansion around $x = 0$:


```python
#Declaration of variables
x         = Symbol('x')
theta     = Symbol('theta', positive = True)

#Declaration of relative risk function
def RR(x, theta):
  return exp(theta/x)

#Declaration of phi for delta method
def AF(x, theta):
  return (RR(x,theta) - 1)/RR(x, theta);

#The derivative at x = 0 exists and is zero for positive theta
derivative = diff(AF(x,theta),x)
limit(derivative, x, 0)
```




$\displaystyle 0$




```python
#This happens for any order (n) of the Taylor expansion as you can verify;
series(AF(x, 1), x0 = 0, n = 10).removeO()
```




$\displaystyle 0$



## Higher order Delta-method
For example, if we wish to estimate the Taylor series for $p\cdot(1-p)$ around $1/2$:


```python
#Get the variables
p    = Symbol('p', positive=True)
phat = RandomSymbol('phat')

#Define the function
def phi(p):
  return p*(1 - p)

#3rd order Taylor series
taylorseries = series(phi(phat), x0 = 1/2, n = 3).removeO()

#get the variance
variance_taylor = Variance(taylorseries).expand()
variance_taylor
```




$\displaystyle \operatorname{Var}\left(\left(\hat{p} - 0.5\right)^{2}\right)$



It is well known that for a sample of size $n$, the variable $B_n = n \cdot \hat{p}$ corresponds to a $\text{Binomial}(n,p)$ this allows us to rewrite the expression above in terms of $B$ to improve the estimation. 

As an example we'll consider $n = 100$ (for now `sympy` can't evaluate with generic $n$)


```python
#CAREFUL! THIS TAKES A WHILE TO EVALUATE
#Declare the variable
n   = 100
B   = Binomial("B_n", n, p)

#Substitute phat = B_n / n
sympy.stats.variance(taylorseries.subs(phat, B/n)).simplify()
```

Finally, we compute the value:

d6ac3229dc660537deb047c74b6106f8ac63277f
ipynb
Jupyter Notebook
DeltaMethod.ipynb
migariane/DeltaMethodInfluenceFunction
fd50ba4fea1d5d6644834e447dc670ffb4a78329
MIT
2022-03-23T09:50:56.000Z
2022-03-23T09:50:56.000Z
DeltaMethod.ipynb
migariane/DeltaMethodInfluenceFunction
fd50ba4fea1d5d6644834e447dc670ffb4a78329
MIT
DeltaMethod.ipynb
migariane/DeltaMethodInfluenceFunction
fd50ba4fea1d5d6644834e447dc670ffb4a78329
MIT
Qwen/Qwen-72B
1. YES
2. YES