Computing and simulating average marginal effect standard error using Delta Method with reproducible codes

Question

I am trying to simulate calculating Average Marginal Effects on a basic linear regression with interaction on a binary variable and compare the empirical standard deviation I get from simulations and the analytic standard error of the Average Marginal Effect I calculate using the Delta Method. But I keep getting a standard error that is way too small.

Below is my attempt with detailed explanation of each step:

1. Data Generation codes:

I am simulating data of the form:

$y = \beta_0 + \beta_1 x_1 + \beta_2 x_1 x_2 + \epsilon$

# import packages
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

def make_data():
    
    # preset beta values
    beta0 = 1
    beta1 = 2
    beta2 = 3
    betas = np.asarray([beta0, beta1, beta2])

    n = 500
    x1 = np.random.choice(2, size=[n, 1])  # x1 is binary!
    x2 = np.random.normal(1, 1, size=[n, 1])
    x1x2 = x1*x2
    X = np.hstack([x1, x2, x1x2])
    epsilon = np.random.normal(0, 0.1, size=[n, 1])
    y = beta0 + x1*beta1 + x1x2*beta2 + epsilon

    data = np.hstack([y, X])
    df = pd.DataFrame(data=data, columns=['y', 'x1', 'x2', 'x1x2'])

    return df, x1, x2, x1x2

df, x1, x2, x1x2 = make_data()

2. Calculate AME (average marginal effect)

I calculate the AME by difference formula for binary variables taken from the official doc of STATA here:

$\frac{1}{n} \sum_{i=1}^n (\hat{\beta}_0 + \hat{\beta}_1 (1) + \hat{\beta}_2 (1) x_2) - (\hat{\beta}_0 + \hat{\beta}_1 (0) + \hat{\beta}_2 (0) x_2) $

def compute_ame(df, x1, x2, x1x2):

    formula = "y ~ x1 + x1x2"
    model = smf.ols(formula, data=df).fit()

    return np.mean( model.predict(df.assign(**{"x1":1, "x1x2":x2})) - model.predict(df.assign(**{"x1":0, "x1x2":0})) )

compute_ame(df, x1, x2, x1x2)  

The output value is around 5, which makes sense since $\beta_1 = 2, \beta_2 = 3$, and $x_2 $ has a mean value of $1$. So from the equation, $y = \beta_0 + \beta_1 x_1 + \beta_2 x_1 x_2 + \epsilon$, turning $x_1$ on increases $y$ by 5.

3. Get the covariance matrix of the model

formula = "y ~ x1 + x1x2"
model = smf.ols(formula, data=df).fit()

vcov = model.cov_params().values

4. Get the gradient of AME with respect to each $\beta$'s

$\frac{\partial}{\partial \beta} \frac{1}{n} \sum_{i=1}^n (\hat{\beta}_0 + \hat{\beta}_1 (1) + \hat{\beta}_2 (1) x_2) - (\hat{\beta}_0 + \hat{\beta}_1 (0) + \hat{\beta}_2 (0) x_2) $

$= \big[ (1-1) \enspace\enspace 1 \enspace\enspace \frac{1}{n} \sum_{i=1}^n x_2 \big]$

We have one row matrix here because we are getting AME for one variable, namely the binary variable $x_1$

grad = np.asarray([0, 1, np.mean(x2)]).reshape([-1, 1])

5. Calculate the standard error:

$(grad)^T (vcov) (grad)$

result = grad.T @ vcov @ grad

6. Run some simulations to get many AME values so we can calculate their standard deviations:

ames = []
for i in range(1000):
    
    df, x1, x2, x1x2 = make_data()
    ame = compute_ame(df, x1, x2, x1x2)
    
    ames.append(ame)
    
np.std(ames)

7. Compare the results:

print(np.std(ames))  # 0.13759581829
print(np.sqrt(result[0][0]))  # 0.0091748250468 from step 5

The analytic standard error is REALLY small.

What am I doing wrong here?

UPDATES (for answer by user4422):

I separate out variable generation from adding Gaussian noise to re-do the simulation as was advised:

def make_data():
    
    # preset beta values
    beta0 = 1
    beta1 = 2
    beta2 = 3
    betas = np.asarray([beta0, beta1, beta2])

    n = 500
    x1 = np.random.choice(2, size=[n, 1])  # x1 is binary!
    x2 = np.random.normal(1, 1, size=[n, 1])
    x1x2 = x1*x2
    X = np.hstack([x1, x2, x1x2])
    # COMMENTED OUT NOISING CODES
    # epsilon = np.random.normal(0, 0.1, size=[n, 1])
    y = beta0 + x1*beta1 + x1x2*beta2 # + epsilon

    data = np.hstack([y, X])
    df = pd.DataFrame(data=data, columns=['y', 'x1', 'x2', 'x1x2'])

    return df, x1, x2, x1x2

def add_noise(df):
    
    n = 500
    epsilon = np.random.normal(0, 0.1, size=n)
    df_noised = df.copy()
    df_noised['y'] += epsilon
    
    return df_noised

The for loop now changes to first generate dataframe and then add new noise in each loop:

ames = []

df, x1, x2, x1x2 = make_data()  # generate variables here first

for i in range(1000):
    
    df_noised = add_noise(df)  # only change the noise component
    ame = compute_ame(df_noised, x1, x2, x1x2)
    ames.append(ame)
    
np.std(ames)

The standard deviations are similar now

print(np.std(ames))  # 0.0091397150909881
print(np.sqrt(result[0][0]))  # 0.0091748250468 from step 5

Answer 1

One thing you are doing wrong is that you do not take the square root of result: you are comparing a variance with a standard deviation.

The second thing, if I understood your code, is that in the computation of result the variance is generated only by the parameter estimation error, while in the computation of np.std(ames) the variance is generated not only by estimation error but also by sampling variability in the variables x1 and x2 (variability in the latter has a double effect in your model, as it contributes also to the average effect). To make a proper comparison you need to re-generate only the error terms epsilon (not the variables) in each simulation.

Let me provide a little more context: the standard errors computed with result are valid in a "fixed regressor" setting (either you have a small sample and the design matrix is fixed, or you have a large sample in which the statistics that depend on the design matrix converge in probability to a constant). If you run a bootstrap (which is what you call a simulation) to reproduce the standard errors of a "fixed regressor" setting, then you only need to randomly draw the errors (the regressors remain fixed). Instead, by randomly drawing also the regressors, you compute the standard errors for a "random regressor" setting. See, for example, this paper (p.254). This is not wrong, but since you want to reproduce the same result with two different methods, you need to keep the "fixed regressor" hypothesis in place.

Also note that, in your setting, not only "random regressors" inflate the variance of the OLS estimates, but they also inflate the variability of the average effect, which depends on the average value of one regressor.