DeclareDesign Community

Ordered probit with control


#1

In response to this post, https://declaredesign.org/blog/2019-02-06-ordered-probit.html, @yukitakahashi11 asks, “What happens if we include controls? Do the results still hold?”

This declaration shows that the answer is yes!

library(DeclareDesign)
library(zeligverse)

# This helper function implements an ordered probit model using zelig; calculates 
# quantities of interest and return tidy output

ordered_probit <-
  
  function(data) {
    zelig_out <- zelig(Y ~ Z + X, model = "oprobit", data = data, cite = FALSE)
    
    sim_out <- sim(zelig_out, x = setx(zelig_out, Z = 0), x1 = setx(zelig_out, Z = 1))
    
    zelig_df <- zelig_qi_to_df(sim_out)
    
    predictions <- with(zelig_df,
                        (X1 + 2 * X2 + 3 * X3 + 4 * X4 + 5 * X5)[Z == 1] -
                          (X1 + 2 * X2 + 3 * X3 + 4 * X4 + 5 * X5)[Z == 0])
    
    return_dat <-
      data.frame(term = "Z",
                 estimate = mean(predictions),
                 std.error = sd(predictions),
                 conf.low = quantile(predictions, .025),
                 conf.high = quantile(predictions, .975))
    return(return_dat)
  }


ordered_probit_design <- 
  
  declare_population(N = 200, X = rnorm(N), noise = rnorm(N)) + 
  
  declare_potential_outcomes(Y ~ draw_ordered(1 * Z + 0.5*X + noise, breaks = c(0, 0.75, 1, 1.25))) +
  
  declare_assignment(prob = 0.5) +
  
  declare_estimand(ate = mean(Y_Z_1 - Y_Z_0)) +
  
  declare_reveal(Y) +
  
  declare_estimator(Y ~ Z + X, model = lm_robust, estimand = "ate", label = "OLS") +
  
  declare_estimator(handler = tidy_estimator(ordered_probit), estimand = "ate", label = "Ordered Probit")


diagnosis <- diagnose_design(ordered_probit_design, sims = 500)

summary_df <- get_diagnosands(diagnosis) %>% 
  gather(key, value, mean_estimate, mean_estimand) %>% 
  mutate(Diagnosand = recode(key, "mean_estimand" = "Mean Estimand", "mean_estimate" = "Mean Estimate"))

diagnosis %>%
  get_simulations() %>%
  ggplot(aes(estimate)) +
  geom_histogram(bins = 30, alpha = 0.6) +
  geom_vline(data = summary_df, aes(xintercept = value, linetype = Diagnosand)) +
  facet_wrap(~ estimator_label) +
  theme_bw() +
  theme(strip.background = element_blank())