teffects ra in R

library(conflicted)
library(tidyverse)

── Attaching core tidyverse packages ────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     

library(haven)
library(boot)
library(tictoc)
library(beepr)
library(car)

Loading required package: carData

library(marginaleffects)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

dat <- read_dta("https://www.stata-press.com/data/r14/cattaneo2.dta")

Try to replicate the example at:

https://www.stata.com/manuals/teteffectsra.pdf#teteffectsra

ATE is -239.6392, SE = 23.82402, z = -10.06, 95% CI = [-286.3334, -192.945].

ATT is -223.3017, SE = 22.7422, z = -9.82, 95% CI = [-267.8755, -178.7278].

{marginaleffects}

Fit the model:

gmod <- lm(bweight ~ mbsmoke * (prenatal1 + mmarried + mage + fbaby), data = dat)

ATE:

avg_comparisons(
  gmod,
  variables = "mbsmoke",
  vcov = "HC1"
) |>
  print(style = "data.frame")

The estimate is spot on, the SE is very close to Stata’s 23.82402 (I think HC1 is the flavour of robust SE provided by Stata).

Check ATT:

avg_comparisons(
  gmod,
  variables = "mbsmoke",
  type = "response",
  newdata = subset(dat, mbsmoke == 1),
  vcov = "HC1"
) |>
  print(style = "data.frame")

Again, the estimate is identical and SE close to Stata’s 22.7422.

Replicate teffects “by hand”

boot_ate <- function(d, i) {
  the_dat <- d[i,]
  dat0 <- the_dat |> dplyr::filter(mbsmoke == 0)
  dat1 <- the_dat |> dplyr::filter(mbsmoke == 1)
  lm0 <- lm(bweight ~ prenatal1 + mmarried + mage + fbaby, data = dat0)
  lm1 <- lm(bweight ~ prenatal1 + mmarried + mage + fbaby, data = dat1)
  
  dat0 <- dat0 |>
    mutate(bweight_0 = bweight,
           bweight_1 = predict(lm1, newdata = dat0))
  
  dat1 <- dat1 |>
    mutate(bweight_1 = bweight,
           bweight_0 = predict(lm0, newdata = dat1))
  
  dat_all <- bind_rows(dat0, dat1) |>
    mutate(te = bweight_1 - bweight_0)
  
  res_ATE <- dat_all |>
    summarise(ate = mean(te))
  
  res_ATT <- dat_all |>
    dplyr::filter(mbsmoke == 1) |>
    summarise(att = mean(te))
  
  c(ate = res_ATE$ate, att = res_ATT$att)
}

I’ve packaged it all up a function ready for bootstrapping, which applied to the original data replicates teffects:

boot_ate(dat, 1:nrow(dat))

      ate       att 
-239.6392 -223.3017 

Now bootstrap boot_ate to get SEs:

set.seed(19082024)
tic()
boots <- boot(data = dat, statistic = boot_ate, R = 1499)
toc()

34.71 sec elapsed

beep(2)

res <- summary(boots) |> as.data.frame()
rownames(res) <- names(boots$t0)
resCI <- bind_cols(res, Confint(boots, type = "perc"))
resCI

Stata said:

ATE is -239.6392, SE = 23.82402, 95% CI = [-286.3334, -192.945].

ATT is -223.3017, SE = 22.7422, 95% CI = [-267.8755, -178.7278].