Advanced Strategies for Inference

Inference w/ non-linearity

\[y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3x_1x_2 + u\]

Calculating measures of uncertainty is tougher with non-linear relationships b/c can’t rely on R canned output

• e.g. what is standard error for marginal effect of \(x_1\) when \(x_2 = 1\)?

Analytic: variance of linear combination

\[ \text{Var} (aX + bY) = a^2\text{Var}(X) + b^2\text{Var}(Y) + 2ab\text{Cov}(X,Y) \]

Standard errors for interactions (see also here)

\[y_i = \beta_0 + \beta_1x_{1i} + \beta_2x_{2i} + \beta_3x_{1i}x_{2i} + u_i\]

Marginal effect of \(x_1\): \(\frac{\partial y_i}{\partial x_{1i}} = \beta_1 + \beta_3x_{2i}\)

Standard errors for quadratic models

Since a quadratic term is a special case of an interaction, this formula holds for those cases as well

• One difference is the “2” multiplier that tags along because of the first derivative of the quadratic term:

\[\sigma_{\frac{\partial y_i}{\partial x_{i}}} = \sqrt{\text{Var}(\hat{\beta}_x) + (2x_{i})^2 \text{Var}(\hat{\beta}_{x^2}) + 4x_{i} \text{Cov}(\hat{\beta}_x, \hat{\beta}_{x^2})}\]

Example

\[y_i = \beta_0 + \beta_1x_{1i} + \beta_2x_{2i} + \beta_3x_{2i}^2 + u_i\]

Interaction example

Call:
lm(formula = econ_mean ~ age_scale + I(age_scale^2) + male + 
    educ_scale + income_scale * know_mean, data = lucid)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.46091 -0.13112 -0.00257  0.11944  0.69693 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)             0.3157633  0.0135866  23.241  < 2e-16 ***
age_scale              -0.1561355  0.0495193  -3.153  0.00163 ** 
I(age_scale^2)          0.1232286  0.0566972   2.173  0.02980 *  
male                    0.0348407  0.0056419   6.175 7.16e-10 ***
educ_scale              0.0257784  0.0111371   2.315  0.02068 *  
income_scale            0.1561915  0.0227687   6.860 7.79e-12 ***
know_mean               0.0007747  0.0165371   0.047  0.96264    
income_scale:know_mean -0.0557548  0.0330965  -1.685  0.09213 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1881 on 4682 degrees of freedom
  (58 observations deleted due to missingness)
Multiple R-squared:  0.05968,   Adjusted R-squared:  0.05827 
F-statistic: 42.45 on 7 and 4682 DF,  p-value: < 2.2e-16

Interaction example

Marginal effect of income_scale:

\[ \frac{\partial y_i}{\partial \text{income_scale_i}} = \beta_5 + \beta_7 \text{know_mean_i} \]

# calculate conditional standard error for marginal effect of income
se_income <- sqrt(   
  
  # (1)^2 * (Var(B_income))
  1^2 * vcov(m5)["income_scale", "income_scale"] +
    
  # (know_mean)^2 * (B_interaction)
  seq(0, 1, 0.05)^2 * vcov(m5)["income_scale:know_mean", "income_scale:know_mean"] + 
    
  # (2) * (1) * (know_mean) * cov(B_income, B_interaction)
  2 * 1 * seq(0, 1, 0.05) * vcov(m5)["income_scale", "income_scale:know_mean"]

)

Plot of conditional marginal effect

Quadratic example

Call:
lm(formula = econ_mean ~ age_scale + I(age_scale^2) + male + 
    educ_scale + income_scale, data = lucid)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.43693 -0.13101 -0.00162  0.11975  0.69088 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     0.322067   0.011080  29.068  < 2e-16 ***
age_scale      -0.158514   0.049534  -3.200  0.00138 ** 
I(age_scale^2)  0.119425   0.056702   2.106  0.03524 *  
male            0.033077   0.005582   5.926 3.33e-09 ***
educ_scale      0.021022   0.010864   1.935  0.05304 .  
income_scale    0.120404   0.010602  11.356  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1882 on 4684 degrees of freedom
  (58 observations deleted due to missingness)
Multiple R-squared:  0.05823,   Adjusted R-squared:  0.05723 
F-statistic: 57.93 on 5 and 4684 DF,  p-value: < 2.2e-16

Plot of conditional marginal effect

marginaleffects package HERE

library(marginaleffects)
options(marginaleffects_print_omit = c("s.value"))

# estimate marginal effects w/ all vars held at means
# if called w/o 'newdata' option, calculates ME for every obs in data
slopes(m5,
       variables = c("male", "educ_scale", "income_scale"),
       newdata = "mean")

         Term Contrast Estimate Std. Error     z Pr(>|z|)   2.5 % 97.5 %
 educ_scale      dY/dX   0.0258    0.01114  2.31   0.0206 0.00395 0.0476
 income_scale    dY/dX   0.1236    0.01066 11.59   <0.001 0.10266 0.1444
 male            1 - 0   0.0348    0.00564  6.18   <0.001 0.02378 0.0459

Type:  response 

At range of values

# estimate marginal effect of income across values of knowledge
m_income <- slopes(m5,
                   variables = c("income_scale"),
                   newdata = datagrid(know_mean = seq(0, 1, 0.1)))
m_income

 know_mean Estimate Std. Error     z Pr(>|z|)  2.5 % 97.5 %
       0.0    0.156     0.0228  6.86   <0.001 0.1116  0.201
       0.1    0.151     0.0199  7.57   <0.001 0.1116  0.190
       0.2    0.145     0.0172  8.43   <0.001 0.1113  0.179
       0.3    0.139     0.0147  9.47   <0.001 0.1106  0.168
       0.4    0.134     0.0127 10.57   <0.001 0.1091  0.159
       0.5    0.128     0.0112 11.43   <0.001 0.1063  0.150
       0.6    0.123     0.0106 11.54   <0.001 0.1019  0.144
       0.7    0.117     0.0111 10.59   <0.001 0.0955  0.139
       0.8    0.112     0.0124  9.01   <0.001 0.0873  0.136
       0.9    0.106     0.0144  7.38   <0.001 0.0779  0.134
       1.0    0.100     0.0168  5.99   <0.001 0.0676  0.133

Term: income_scale
Type:  response 
Comparison: dY/dX

Plot

Average marginal effects

avg_slopes calculates marginal effect at every observation and then takes mean

avg_slopes(m5,
           variables = c("age_scale", "male", "educ_scale", "income_scale")
)

         Term Contrast Estimate Std. Error     z Pr(>|z|)    2.5 %  97.5 %
 age_scale       dY/dX  -0.0479    0.01277 -3.75   <0.001 -0.07294 -0.0229
 educ_scale      dY/dX   0.0258    0.01114  2.31   0.0207  0.00394  0.0476
 income_scale    dY/dX   0.1236    0.01066 11.59   <0.001  0.10266  0.1444
 male            1 - 0   0.0348    0.00564  6.18   <0.001  0.02378  0.0459

Type:  response 

Can also do slope tests

Imagine we want to test two coefficients: null hypothesis is equality

summary(m5)

Call:
lm(formula = econ_mean ~ age_scale + I(age_scale^2) + male + 
    educ_scale + income_scale * know_mean, data = lucid)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.46091 -0.13112 -0.00257  0.11944  0.69693 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)             0.3157633  0.0135866  23.241  < 2e-16 ***
age_scale              -0.1561355  0.0495193  -3.153  0.00163 ** 
I(age_scale^2)          0.1232286  0.0566972   2.173  0.02980 *  
male                    0.0348407  0.0056419   6.175 7.16e-10 ***
educ_scale              0.0257784  0.0111371   2.315  0.02068 *  
income_scale            0.1561915  0.0227687   6.860 7.79e-12 ***
know_mean               0.0007747  0.0165371   0.047  0.96264    
income_scale:know_mean -0.0557548  0.0330965  -1.685  0.09213 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1881 on 4682 degrees of freedom
  (58 observations deleted due to missingness)
Multiple R-squared:  0.05968,   Adjusted R-squared:  0.05827 
F-statistic: 42.45 on 7 and 4682 DF,  p-value: < 2.2e-16

hypotheses(m5, hypothesis = "male - educ_scale = 0")

        Hypothesis Estimate Std. Error     z Pr(>|z|)   2.5 % 97.5 %
 male-educ_scale=0  0.00906     0.0128 0.711    0.477 -0.0159 0.0341

Simulation for inference

In this context, “simulation” describes the process of repeated random sampling from the sampling distribution

Using sim from arm

# draw 1,000 realizations of model parameters from sampling distribution
m4_sims <- arm::sim(m4, n.sims=1000)

# the betas are stored in @coef
m4_sims@coef[1:3,]

     (Intercept)  age_scale I(age_scale^2)       male   educ_scale income_scale
[1,]   0.3201757 -0.1165536     0.06712922 0.03230549  0.023721840    0.1156481
[2,]   0.3084358 -0.1253236     0.07547290 0.03793351  0.019553911    0.1313397
[3,]   0.3538939 -0.2406292     0.17863659 0.04540851 -0.001600411    0.1216964

# the sigmas are stored in @sigma
m4_sims@sigma[1:3]

[1] 0.1872169 0.1880067 0.1916445

What does sim do?

For each sample, sim first draws a random value of \(\hat{\sigma}\) from its sampling distribution, and then uses this value to draw a random sample from the sampling distribution for \(\hat{\boldsymbol{\beta}}\)

sigma.hat <- summary(m4)$sigma
beta.hat <- summary(m4)$coef[, 1, drop=FALSE]
cov.beta <- summary(m4)$cov.unscaled
n <- summary(m4)$df[1] + summary(m4)$df[2]
k <- summary(m4)$df[1]
n.sims <- 1000

sigma <- array(NA, n.sims)
beta <- array(NA, c(n.sims, k))
for (i in 1:n.sims){
  sigma[i] <- sigma.hat*sqrt((n-k)/rchisq(1,n-k))
  beta[i,] <- MASS::mvrnorm (1, beta.hat, cov.beta*sigma[i]^2)
}

sim draws from MV normal

This is the asymptotic distribution of \(\hat{\boldsymbol{\beta}}\) under standard assumptions, and the t distribution converges pretty quickly to normal with increasing sample size

• For small samples, say, less than 100, you might not want to rely on asymptotic normality

• A good alternative is bootstrapping (described below), which does not make distributional assumptions

Plot

Resampling methods

Resampling methods use the actual data you have to approximate the sampling distribution of quantities of interest

Advantages and disadvantages

Jackknife

A leave-one-out estimator computes estimates using all observations except one (there are thus \(N\) possibilities)

Jackknife

Very flexible method, overlap with other approaches

Example with loop

Draw boots random samples from analytical sample, each of size \(N\), run model, store coefs

# number of bootstrap samples
boots <- 1000

# initialize array to store coefs
boots_beta <- array(NA, c(boots, # num of bootstrap samples
                          summary(m4)$df[1]) # number of regression coefficients (K+1)
)

# iterate over samples
for (i in 1:boots){
  
  new <- sample(1:nrow(m4$model), nrow(m4$model), replace = T) # sample rows of data w/ replacement
  
  boots_beta[i, ] <- coef(lm(formula(m4), data = m4$model[new, ])) # run model with bootstrap sample
 
  # repeat boots number of times 
}

Standard errors

The distribution of bootstrap estimates is an estimate of the sampling distribution of the parameter, which we can then use for inference

• The bootstrap standard error is simply the square root of the sample variance of bootstrap estimates
• It is consistent for simple linear models

Confidence intervals

There is a large literature on methods for calculating bootstrap confidence intervals, and various different options

• normal approximation
• percentile intervals
• bias-corrected
• percentile-t

Normal approximation

Use bootstrap standard error and plug into formula for confidence intervals using normal distribution

• Often the default
• Not the best option: worse performance than others

Percentile intervals

Choose confidence level (e.g., \(\alpha = 5%\)), select quantiles of bootstrap estimates associated with \(\alpha/2\) and \((100 - \alpha/2)\) (e.g., 2.5% and 97.5%)

• Easy and popular; does not require calculation of standard errors
• Consistent under weaker assumptions than normal approximation and bootstrap standard errors, and does not require trimming

Bias-corrected (BC) intervals

Adjusts percentile intervals for finite-sample bias using degree of asymmetry in empirical distribution of bootstrap estimates around parameter estimate

• Provides “exact” coverage under assumption that bias depends only on a parameter drawn from a symmetric, invertible distribution
• Is also asymptotically valid for the interval of interest

Percentile-t intervals

Calculate a t-statistic for every bootstrap sample: \(T^* = \frac{\hat{\theta}^* - \hat{\theta}}{se(\hat{\theta}^*)}\)

• \(\hat{\theta}^*\) is the bootstrap parameter estimate for a particular bootstrap sample
• \(se(\hat{\theta}^*)\) is the estimated standard error for a particular bootstrap sample
• \(\hat{\theta}\) is the estimate from your standard procedure (e.g., OLS coefficient)

Percentile-t intervals

Calculate the percentile-t interval from the quantiles of \(T^*\), \(q^*\):

\[ \left( \hat{\theta} - se(\hat{\theta}) q_{(1 - \alpha/2)}^*, \space \hat{\theta} - se(\hat{\theta}) q_{(\alpha/2)}^* \right) \]

• Uses standard error, so need a method for calculating, and accuracy depends on the reliability of the estimate

  • Can use bootstrap to estimate SE on each bootstrap (i.e., nested bootstrap)

• When SE estimate reliable, may be better than other methods

• Faster rate of asymptotic convergence than percentile and BC methods

Using boot package via marginaleffects

# estimate slopes
m5_me <- slopes(m5,
                variables = c("male", "educ_scale", "income_scale"),
                newdata = "mean")

# use inferences() to get bootstrap intervals
inferences(m5_me, method = "boot", R = 1000, conf_type = "perc")
inferences(m5_me, method = "boot", R = 1000, conf_type = "bca")

boot package to get percentile-t

# boot call
m5_boot <- boot(m5$model, 
                function(data, indices) {
                  d <- data[indices, ]
                  fit <- update(m5, . ~ ., data = d)
                  B <- coef(fit)[5]
                  se <- sqrt(diag(vcov(fit)))[5]
                  return(c(B, se))
                },
                R = 1000)

# use boot.ci to get CIs
boot.ci(m5_boot, type = "stud")

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates

CALL : 
boot.ci(boot.out = m5_boot, type = "stud")

Intervals : 
Level    Studentized     
95%   ( 0.0033,  0.0467 )  
Calculations and Intervals on Original Scale

Bootstrap hypothesis testing

We can use the bootstrapped t-statistics to calculate p-values - recall:

\[T^* = \frac{\hat{\theta}^* - \hat{\theta}}{se(\hat{\theta}^*)}\]

• The (two-tailed) p-value is then the percent of \(|T^*| > |T|\), where \(T\) is the t-statistic from your standard procedure (e.g., OLS output for a coefficient)

Wild bootstrap

Our bootstrap methods to this point all sample both \(X\) and \(Y\), but our standard linear regression model treats \(X\) as fixed and imposes \(\text{E}(u | \mathbf{X}) = 0\)

• Integrating these assumptions into the bootstrap will increase precision (and also take into account heteroskedasticity)

Structured data (see Cameron et al. 2008)

If you have data with “structure” that is relevant to the analysis, you need to reproduce that structure in your resampling

How many?

For simulation-based and bootstrapped calculations, how many should you use?

• As always with these things, no firm answer beyond “as many as possible given constraints”

Set seed!

Remember to set a seed value (especially for your final runs that will be archived) (set.seed())

Example

SE for conditional ME of age: \(f(\hat{\boldsymbol{\beta}}) = \hat{\beta}_1 + 2\hat{\beta}_2age_i\)

Non-linear example

m1 <- lm(log(salary) ~ yrs.since.phd + sex + discipline, 
         data = carData::Salaries)
summary(m1)

Call:
lm(formula = log(salary) ~ yrs.since.phd + sex + discipline, 
    data = carData::Salaries)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.84527 -0.15697 -0.00855  0.15550  0.62670 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.125e+01  4.189e-02 268.671  < 2e-16 ***
yrs.since.phd 9.584e-03  9.075e-04  10.560  < 2e-16 ***
sexMale       6.639e-02  3.830e-02   1.734   0.0838 .  
disciplineB   1.447e-01  2.319e-02   6.239 1.14e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2244 on 393 degrees of freedom
Multiple R-squared:  0.2615,    Adjusted R-squared:  0.2559 
F-statistic: 46.39 on 3 and 393 DF,  p-value: < 2.2e-16

Non-linear example

Avg effect of change from 10 to 20 years since PhD on unlogged salaries: \(f = \frac{1}{N} \sum_N e^{\boldsymbol{x}_{i,yrs=20}\hat{\boldsymbol{\beta}}}E[u_i] - e^{\boldsymbol{x}_{i,yrs=10}\hat{\boldsymbol{\beta}}}E[u_i]\)

\(\frac{\partial f}{\partial \boldsymbol{\beta}} =\)

\[ \begin{bmatrix} \frac{1}{N} \sum_N e^{\boldsymbol{x}_{i,yrs=20}\hat{\boldsymbol{\beta}}}E[u_i] - e^{\boldsymbol{x}_{i,yrs=10}\hat{\boldsymbol{\beta}}}E[u_i] \\ \frac{1}{N} \sum_N 20e^{\boldsymbol{x}_{i,yrs=20}\hat{\boldsymbol{\beta}}}E[u_i] - 10e^{\boldsymbol{x}_{i,yrs=10}\hat{\boldsymbol{\beta}}}E[u_i] \\ \frac{1}{N} \sum_N male_ie^{\boldsymbol{x}_{i,yrs=20}\hat{\boldsymbol{\beta}}}E[u_i] - male_ie^{\boldsymbol{x}_{i,yrs=10}\hat{\boldsymbol{\beta}}}E[u_i] \\ \frac{1}{N} \sum_N discipline_ie^{\boldsymbol{x}_{i,yrs=20}\hat{\boldsymbol{\beta}}}E[u_i] - discipline_ie^{\boldsymbol{x}_{i,yrs=10}\hat{\boldsymbol{\beta}}}E[u_i] \end{bmatrix} \]

Calculation

# print

cbind(fd, se)

           fd         
[1,] 9973.553 883.2016

Using car package

Holding all other variables at central tendencies:

library(car)

Eu_1 <- sum(exp(m1$residuals)) / nrow(m1$model)

car <- deltaMethod(m1,

                   "exp(b0 + b1*20 + b2*0 + b3*0)*Eu_1 - exp(b0 + b1*10 + b2*0 + b3*0)*Eu_1",

                  parameterNames = c("b0","b1","b2","b3"))

cbind(car$Estimate, car$SE)

         [,1]     [,2]
[1,] 8765.543 861.2572

marginaleffects package

The delta method is the default approach to calculating standard errors in marginaleffects

• I do not see a simple approach for logged DVs, but I think you can write your own functions (or use car)

• But for most cases, it will work well

• The exception is when you have exact standard errors (assuming normal regression errors), in which case, you would use t for tests and CIs, but delta method uses normal