Testing cluster + pweights

Setup

data(ai_data, package = 'aiHelper')
df <- ai_data
set.seed(100)
df$wgt <- runif(nrow(df), 0.5, 1.5)

Fit OLS models

Stata

s <- aiEstimation::ai_toplines(
  outcome = 'ag_fav',
  treat = 'treat',
  model = 'reg',
  df = df,
  covar = c('age', 'partisanship_bin'),
  extra_str = '[pw = wgt], vce(cluster reg_zip5)',
  extra_var = c('wgt', 'reg_zip5')
)

## Going from 10000 to 10000 rows (dropping 0) due to missing values.

R

mod <- glm(ag_fav ~ treat + age + partisanship_bin, 
           weights = wgt,
           data = df)

mr <- list(
  pred = modmarg::marg(
    mod, 
    var_interest = 'treat', 
    weights = mod$prior.weights,
    vcov_mat = aiEstimation::cluster_se(
      mod, cluster = 'reg_zip5', data = mod$data),
    dof = length(unique(mod$data$reg_zip5)) - 1
  )[[1]],
  effects = modmarg::marg(
    mod,
    var_interest = 'treat',
    weights = mod$prior.weights,
    vcov_mat = aiEstimation::cluster_se(
      mod, cluster = 'reg_zip5', data = mod$data),
    dof = length(unique(mod$data$reg_zip5)) - 1,
    type = 'effects'
  )[[1]]
)

## Std. Err. adjusted for 9457 clusters in reg_zip5
## 
## Std. Err. adjusted for 9457 clusters in reg_zip5

Compare

Stata

for(i in seq_along(s$margins)){
  temp <- s$margins[[i]][, c('treat_label', 'margins_b', 'margins_se')]
  temp$treat_label <- sprintf('treat = %s', temp$treat_label)
  names(temp) <- c('Label', 'Margin', 'Standard.Error')
  print(knitr::kable(temp, digits = 7))
}

Label	Margin	Standard.Error
treat = control	0.4434034	0.0070616
treat = treat	0.5813327	0.0070346

Label	Margin	Standard.Error
treat = control	0.0000000	NA
treat = treat	0.1379293	0.0099973

R

for(i in seq_along(mr)){
  # Round to match stata output
  print(knitr::kable(
    mr[[i]][, c('Label', 'Margin', 'Standard.Error')], 
    digits = 7))
}

Label	Margin	Standard.Error
treat = control	0.4434034	0.0070616
treat = treat	0.5813327	0.0070346

Label	Margin	Standard.Error
treat = control	0.0000000	0.0000000
treat = treat	0.1379293	0.0099973

Fit Logit models

Stata

s <- aiEstimation::ai_toplines(
  outcome = 'ag_fav',
  treat = 'treat',
  model = 'logit',
  df = df,
  covar = c('age', 'partisanship_bin'),
  extra_str = '[pw = wgt], vce(cluster reg_zip5)',
  extra_var = c('wgt', 'reg_zip5')
)

## Going from 10000 to 10000 rows (dropping 0) due to missing values.

R

mod <- glm(ag_fav ~ treat + age + partisanship_bin, 
           weights = wgt,
           family = 'binomial',
           data = df)

## Warning: non-integer #successes in a binomial glm!

mr <- list(
  pred = modmarg::marg(
    mod, 
    var_interest = 'treat', 
    weights = mod$prior.weights,
    vcov_mat = aiEstimation::cluster_se(
      mod, cluster = 'reg_zip5', data = mod$data),
    dof = length(unique(mod$data$reg_zip5)) - 1
  )[[1]],
  effects = modmarg::marg(
    mod,
    var_interest = 'treat',
    weights = mod$prior.weights,
    vcov_mat = aiEstimation::cluster_se(
      mod, cluster = 'reg_zip5', data = mod$data),
    dof = length(unique(mod$data$reg_zip5)) - 1,
    type = 'effects'
  )[[1]]
)

## Std. Err. adjusted for 9457 clusters in reg_zip5
## 
## Std. Err. adjusted for 9457 clusters in reg_zip5

Compare

Stata

for(i in seq_along(s$margins)){
  temp <- s$margins[[i]][, c('treat_label', 'margins_b', 'margins_se')]
  temp$treat_label <- sprintf('treat = %s', temp$treat_label)
  names(temp) <- c('Label', 'Margin', 'Standard.Error')
  print(knitr::kable(temp, digits = 7))
}

Label	Margin	Standard.Error
treat = control	0.4434760	0.0070669
treat = treat	0.5812349	0.0070304

Label	Margin	Standard.Error
treat = control	0.0000000	NA
treat = treat	0.1377589	0.0099936

R

for(i in seq_along(mr)){
  # Round to match stata output
  print(knitr::kable(
    mr[[i]][, c('Label', 'Margin', 'Standard.Error')], 
    digits = 7))
}

Label	Margin	Standard.Error
treat = control	0.4434760	0.0070669
treat = treat	0.5812349	0.0070304

Label	Margin	Standard.Error
treat = control	0.0000000	0.0000000
treat = treat	0.1377589	0.0099936