# Ordered probit with control

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())

``````