## ----setup, include = FALSE-----------------------------------------------
Sys.setenv(OMP_NUM_THREADS = "1", OPENBLAS_NUM_THREADS = "1",
           MKL_NUM_THREADS = "1", VECLIB_MAXIMUM_THREADS = "1",
           OMP_THREAD_LIMIT = "1")
options(prompt = "R> ", continue = "+ ", width = 76)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5,
  fig.align = "center"
)

## ----install, eval = FALSE------------------------------------------------
# install.packages("nonprobsvy")

## ----libs, message = FALSE, warning = FALSE-------------------------------
library("nonprobsvy")
library("ggplot2")

## ----jvs-data-------------------------------------------------------------
data("jvs", package = "nonprobsvy")
head(jvs)

## ----jvs-design-----------------------------------------------------------
jvs_svy <- svydesign(ids = ~ 1,
                     weights = ~ weight,
                     strata = ~ size + nace + region,
                     data = jvs)

## ----admin-data-----------------------------------------------------------
data("admin", package = "nonprobsvy")
head(admin)

## ----ipw1-----------------------------------------------------------------
ipw_est1 <- nonprob(
  selection = ~ region + private + nace + size,
  target = ~ single_shift,
  svydesign = jvs_svy,
  data = admin,
  method_selection = "logit"
)

## ----ipw1-print-----------------------------------------------------------
ipw_est1

## ----ipw1-summary---------------------------------------------------------
summary(ipw_est1)

## ----ipw1-extract---------------------------------------------------------
extract(ipw_est1)

## ----ipw2-----------------------------------------------------------------
ipw_est2 <- nonprob(
  selection = ~ region + private + nace + size,
  target = ~ single_shift,
  svydesign = jvs_svy,
  data = admin,
  method_selection = "logit",
  control_selection = control_sel(gee_h_fun = 1, est_method = "gee")
)

## ----ipw2-print-----------------------------------------------------------
ipw_est2

## ----ipw-balance----------------------------------------------------------
data.frame(ipw_mle = check_balance(~ size - 1, ipw_est1, 1)$balance,
           ipw_gee = check_balance(~ size - 1, ipw_est2, 1)$balance)

## ----mi1------------------------------------------------------------------
mi_est1 <- nonprob(
  outcome = single_shift ~ region + private + nace + size,
  svydesign = jvs_svy,
  data = admin,
  method_outcome = "glm",
  family_outcome = "binomial"
)
mi_est1

## ----mi23-----------------------------------------------------------------
mi_est2 <- nonprob(
  outcome = single_shift ~ region + private + nace + size,
  svydesign = jvs_svy,
  data = admin,
  method_outcome = "nn",
  control_outcome = control_out(k = 5)
)
mi_est3 <- nonprob(
  outcome = single_shift ~ region + private + nace + size,
  svydesign = jvs_svy,
  data = admin,
  method_outcome = "pmm",
  family_outcome = "binomial",
  control_outcome = control_out(k = 5)
)

## ----mi23-extract---------------------------------------------------------
rbind("NN" = extract(mi_est2)[, 2:3], "PMM" = extract(mi_est3)[, 2:3])

## ----dr1------------------------------------------------------------------
dr_est1 <- nonprob(
  selection = ~ region + private + nace + size,
  outcome = single_shift ~ region + private + nace + size,
  svydesign = jvs_svy,
  data = admin,
  method_selection = "logit",
  method_outcome = "glm",
  family_outcome = "binomial"
)
dr_est1

## ----dr1-summary----------------------------------------------------------
summary(dr_est1)

## ----dr2------------------------------------------------------------------
set.seed(2026)
dr_est2 <- nonprob(
  selection = ~ region + private + nace + size,
  outcome = single_shift ~ region + private + nace + size,
  svydesign = jvs_svy,
  data = admin,
  method_selection = "logit",
  method_outcome = "glm",
  family_outcome = "binomial",
  control_selection = control_sel(nfolds = 3, nlambda = 10),
  control_outcome = control_out(nfolds = 3, nlambda = 10),
  control_inference = control_inf(bias_correction = TRUE,
                                  vars_combine = TRUE,
                                  vars_selection = TRUE)
)
dr_est2

## ----comparison-of-est, fig.cap = "Figure 1: Comparison of estimates of the share of job vacancies offered on a single-shift."----
df_s <- rbind(extract(ipw_est1), extract(ipw_est2), extract(mi_est1),
              extract(mi_est2), extract(mi_est3), extract(dr_est1),
              extract(dr_est2))

df_s$est <- c("IPW (MLE)", "IPW (GEE)", "MI (GLM)", "MI (NN)",
              "MI (PMM)", "DR", "DR (BM)")

ggplot(data = df_s,
       aes(y = est, x = mean, xmin = lower_bound, xmax = upper_bound)) +
  geom_point() +
  geom_vline(xintercept = mean(admin$single_shift),
             linetype = "dotted", color = "red") +
  geom_errorbar() +
  labs(x = "Point estimator and confidence interval", y = "Estimators") +
  theme_bw()

## ----ipw-boot-------------------------------------------------------------
set.seed(2024)
ipw_est1_boot <- nonprob(
  selection = ~ region + private + nace + size,
  target = ~ single_shift,
  svydesign = jvs_svy,
  data = admin,
  method_selection = "logit",
  control_inference = control_inf(var_method = "bootstrap", num_boot = 50),
  verbose = FALSE
)

## ----ipw-boot-compare-----------------------------------------------------
rbind("IPW analytic variance"  = extract(ipw_est1)[, 2:3],
      "IPW bootstrap variance" = extract(ipw_est1_boot)[, 2:3])

## ----ipw-boot-sample------------------------------------------------------
head(ipw_est1_boot$boot_sample, n = 3)

## ----mi-sel---------------------------------------------------------------
set.seed(2024)
mi_est1_sel <- nonprob(
  outcome = single_shift ~ region + private + nace + size,
  svydesign = jvs_svy,
  data = admin,
  method_outcome = "glm",
  family_outcome = "binomial",
  control_outcome = control_out(nfolds = 3, nlambda = 10, penalty = "lasso"),
  control_inference = control_inf(vars_selection = TRUE),
  verbose = TRUE
)

## ----mi-sel-compare-------------------------------------------------------
rbind("MI without var sel" = extract(mi_est1)[, 2:3],
      "MI with var sel"    = extract(mi_est1_sel)[, 2:3])

## ----mi-sel-coef----------------------------------------------------------
round(coef(mi_est1_sel)$coef_out[, 1], 4)

## ----ipw-coef-------------------------------------------------------------
round(coef(ipw_est1)$coef_sel[, 1], 4)

## ----nobs-----------------------------------------------------------------
nobs(dr_est1)

## ----confint--------------------------------------------------------------
confint(dr_est1, level = 0.99)

## ----weights--------------------------------------------------------------
summary(weights(dr_est1))

## ----method-glm-----------------------------------------------------------
res_glm <- method_glm(
  y_nons = admin$single_shift,
  X_nons = model.matrix(~ region + private + nace + size, admin),
  X_rand = model.matrix(~ region + private + nace + size, jvs),
  svydesign = jvs_svy)
res_glm

## ----method-ps------------------------------------------------------------
method_ps()

## ----method-ps-help, eval = FALSE-----------------------------------------
# ?method_ps()