## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  fig.align = "center",
  out.width = "100%",      # <-- Forces the image to scale to the text width
  dpi = 150,                # <-- Optimized for CRAN (read below)
  message = FALSE,
  warning = FALSE
)

## ----eval=FALSE---------------------------------------------------------------
# install.packages("vbm")

## ----results='hide',message=FALSE,warning=FALSE-------------------------------
library(vbm)

## ----eval=FALSE---------------------------------------------------------------
# # Install development version from GitHub
# devtools::install_github("Staniks0/vbm")

## ----load-data----------------------------------------------------------------
library(WeightIt)
library(dplyr)
library(ggplot2)
library(jointVIP)
library(cobalt)
library(vbm)
data(NCDS)
str(ncds)

## ----weighting----------------------------------------------------------------
weightlist_att <- weightit(
  formula = Dany ~ white + maemp + scht +
    qmab + qmab2 +
    qvab + qvab2 +
    paed_u + maed_u +
    agepa + agema +
    sib_u,
  data = ncds,
  method = "ebal",
  estimand = "ATT",
  stabilize = FALSE,
  include.obj = TRUE
)

summary(weightlist_att)
love.plot(weightlist_att, stars="std")


## -----------------------------------------------------------------------------
set.seed(123)
pilot_prop = 0.1
controls <- which(ncds$Dany == 0)
pilot_indices <- sample(controls, length(controls) * pilot_prop)

pilot_df <- ncds[pilot_indices, ]
analysis_data <- ncds[-pilot_indices, ]
analysis_data1 <- ncds[-pilot_indices, ]
post_df <- weightit(
  formula = Dany ~ white + maemp + scht +
    qmab + qmab2 +
    qvab + qvab2 +
    paed_u + maed_u +
    agepa + agema +
    sib_u,
  data = analysis_data,
  method = "ebal",
  estimand = "ATT",
  stabilize = FALSE,
  include.obj = TRUE
)
new_jointVIP <- create_jointVIP(
  treatment = "Dany",           
  outcome = "wage",             
  covariates = c("white", "maemp", "scht", "qmab", "qmab2", "qvab",
                 "qvab2", "paed_u", "maed_u", "agepa", "agema", "sib_u"), 
  pilot_df = pilot_df,
  analysis_df = analysis_data    
)
plot(create_post_jointVIP(new_jointVIP,
                                      post_analysis_df = analysis_data),
     plot_title = "Pre-weighting jointVIP: Educational Attainment Analysis")

## -----------------------------------------------------------------------------
plot(create_post_jointVIP(new_jointVIP,
                                      post_analysis_df = analysis_data,
                                      wts = post_df$weights),
     plot_title = "Post-weighting jointVIP: Educational Attainment Analysis")



## ----estimate-att-------------------------------------------------------------
# Run the sensitivity analysis
vbm_results <- run_sensitivity_analysis(
  weightlist = weightlist_att,
  Y = "wage",
  treatment = "Dany",
  data = ncds,
  approach = "vbm",
  n_bootstrap = 1000,
  seed = 331,
  n_cores = 1,
  benchmarking = TRUE
)

vbm_results$ipw_estimate
vbm_results$ipw_se
vbm_results$difference_in_means
vbm_results$difference_in_means_se

## ----sensitivity-analysis-----------------------------------------------------
head(vbm_results$point_bounds)
head(vbm_results$benchmark_parameters)
plot_sensitivity_analysis(
  vbm_results,
  parameter_level = seq(0, 0.15, by = 0.01),
  benchmarking = TRUE,
  benchmark_variable = c("scht", "qmab2", "agema", "qvab2", "white", "qmab", "qvab"),
  variable_name = c("Selective School", "Math Score at 11", "Mother's Age",
                     "Reading Score at 11", "White", 
                    "Math Score at 7", "Reading Score at 7"),
  header = "Effect of Education Attainment on Hourly Wages by VBM"
)

## ----sensitivity-analysis-msm-------------------------------------------------
msm_results <- run_sensitivity_analysis(
  weightlist = weightlist_att,
  Y = "wage",
  treatment = "Dany",
  data = ncds,
  approach = "msm",
  n_bootstrap = 1000,
  seed = 331,
  n_cores = 1,
  benchmarking = TRUE
)
msm_results$lambda_star
head(msm_results$point_bounds)
msm_results$benchmark_parameters
plot_sensitivity_analysis(
  msm_results,
  parameter_level = seq(1, 5.5, by = 0.25),
  benchmarking = TRUE,
  benchmark_variable = c("maemp", "qvab", "scht", "sib_u", "qvab2", "qmab2", "qmab"),
  variable_name = c("Mother's Employment", "Reading Score at 7", "Selective School", 
                    "No. of Siblings", "Reading Score at 11", "Math Score at 11", 
                    "Math Score at 7"),
  header = "Effect of Education Attainment on Hourly Wages by MSM"
)

## ----sensitivity-analysis-vbm-w-cor-------------------------------------------
# Run the sensitivity analysis
vbm_w_cor_results <- run_sensitivity_analysis(
  weightlist = weightlist_att,
  Y = "wage",
  treatment = "Dany",
  data = ncds,
  approach = "vbm_w_cor",
  R2_seq = seq(0, 0.8, by = 0.01),
  cor_eps_seq = rep(0.8, 81),
  n_bootstrap = 1000,
  n_cores = 1,
  seed = 331
)

# Critical lambda*
vbm_w_cor_results$Rstar
# Point bounds by vbm with less conservative bounds
head(vbm_w_cor_results$point_bounds)
# Point bounds by vbm with worst case bounds
head(vbm_results$point_bounds)
# Point bounds by msm with worst case bounds
head(msm_results$point_bounds)
# Confidence intervals by vbm with less conservative bounds
head(vbm_w_cor_results$confidence_intervals)
# Confidence intervals by vbm with worst case bounds
head(vbm_results$confidence_intervals)
# Confidence intervals by msm with worst case bounds
head(msm_results$confidence_intervals)