## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(dplyr)
library(rstan)
library(adjustr)
load("eightschools_model.rda")

## ----eval=F-------------------------------------------------------------------
# library(dplyr)
# library(rstan)
# library(adjustr)
# 
# model_code = "
# data {
#     int<lower=0> J;         // number of schools
#     real y[J];              // estimated treatment effects
#     real<lower=0> sigma[J]; // standard error of effect estimates
# }
# parameters {
#     real mu;                // population treatment effect
#     real<lower=0> tau;      // standard deviation in treatment effects
#     vector[J] eta;          // unscaled deviation from mu by school
# }
# transformed parameters {
#     vector[J] theta = mu + tau * eta;        // school treatment effects
# }
# model {
#     eta ~ std_normal();
#     y ~ normal(theta, sigma);
# }"
# 
# model_d = list(J = 8,
#                y = c(28,  8, -3,  7, -1,  1, 18, 12),
#                sigma = c(15, 10, 16, 11,  9, 11, 10, 18))
# eightschools_m = stan(model_code=model_code, chains=2, data=model_d,
#                       warmup=500, iter=1000)

## -----------------------------------------------------------------------------
plot(eightschools_m, pars="theta")

## -----------------------------------------------------------------------------
spec = make_spec(mu ~ normal(0, 20), tau ~ exponential(5))
print(spec)

## ----include=F----------------------------------------------------------------
adjusted = adjust_weights(spec, eightschools_m, keep_bad=TRUE)

## ----eval=F-------------------------------------------------------------------
# adjusted = adjust_weights(spec, eightschools_m)

## -----------------------------------------------------------------------------
print(adjusted)

## -----------------------------------------------------------------------------
summarize(adjusted, mean(mu), var(mu))

## -----------------------------------------------------------------------------
spec = make_spec(eta ~ student_t(df, 0, 1), df=1:10)
print(spec)

## -----------------------------------------------------------------------------
adjusted = adjust_weights(spec, eightschools_m, keep_bad=TRUE)

## -----------------------------------------------------------------------------
print(adjusted)

## -----------------------------------------------------------------------------
spec_plot(adjusted, df, mu)
spec_plot(adjusted, df, theta[3])

## -----------------------------------------------------------------------------
summarize(adjusted, wasserstein(mu))