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

## ----load-data----------------------------------------------------------------
library(nonprobsampling)

data("sc")
data("sp1")
data("sp1_bootstrap")
data("sp2")

str(sc)
str(sp1)
str(sp2)

## ----peek-data-sc-------------------------------------------------------------
head(sc)

## ----PSU-survey-designs-------------------------------------------------------
ref1_design <- survey::svydesign(
  ids     = ~psu_sp1,
  strata  = ~strata_sp1,
  weights = ~wts_sp1,
  data    = sp1,
  nest    = TRUE
)

ref2_design <- survey::svydesign(
  ids     = ~psu_sp2,
  strata  = ~strata_sp2,
  weights = ~wts_sp2,
  data    = sp2,
  nest    = TRUE
)

## ----replicate-wts-survey-design----------------------------------------------
ref1_design_rep_wts <- survey::svrepdesign(
   data    = sp1_bootstrap,
   weights = ~wts_sp1,
   repweights = "bw[0-9]+",
   type    = "bootstrap"
)


## ----one-ref-calibration-classic----------------------------------------------

fit_cali <- est_pw(
  data      = list(sc, ref1_design),
  p_formula = ~ agecat + race + education + comorbidity + BMI + diabetes,
  method    = "calibration",
  control   = pw_solver_control(ftol = 1e-6)
)

print(fit_cali)

summary(fit_cali)

## ----one-ref-calibration-autopfml---------------------------------------------
fit_cali_auto <- est_pw(
  data   = list(sc, ref1_design),
  method = "calibration",
  control   = pw_solver_control(ftol = 1e-6)
)
summary(fit_cali_auto)

## ----one-ref-calibration-int--------------------------------------------------

fit_cali_interaction <- est_pw(
  data      = list(sc, ref1_design),
  p_formula = ~ BMI * psa_level + diabetes + pros_enlarged + I(psa_level^2),
  method    = "calibration",
  control   = pw_solver_control(ftol = 1e-6)
)

summary(fit_cali_interaction)

## ----one-ref-calibration-verbose----------------------------------------------
fit_cali <- est_pw(
  data      = list(sc, ref1_design),
  p_formula = ~ agecat + race + education + comorbidity + BMI + diabetes,
  method    = "calibration",
  verbose   = TRUE,
  control   = pw_solver_control(ftol = 1e-6)
)

## ----one-ref-calibration-sc_updated-------------------------------------------
head(fit_cali$sc_updated)

## ----one-ref-calibration-pseudo_weights---------------------------------------
head(fit_cali$pseudo_weights)

## ----one-ref-alp--------------------------------------------------------------
fit_alp <- est_pw(
  data      = list(sc, ref1_design_rep_wts),
  p_formula = ~ agecat + race + education + comorbidity + BMI + diabetes,
  method    = "alp"
)

print(fit_alp)

## ----one-ref-clw--------------------------------------------------------------
fit_clw <- est_pw(
  data      = list(sc, ref1_design),
  p_formula = ~ agecat + race + education + comorbidity + BMI + diabetes,
  method    = "clw"
)

print(fit_clw)

## ----multi-ref-by-size--------------------------------------------------------
fit_multi_order_by_size <- est_pw(
  data      = list(sc, ref1_design, ref2_design),
  p_formula = list(
    ~ agecat + race + education + psa_level  + pros_enlarged + comorbidity + BMI + diabetes,
    ~ agecat + race + diabetes  + BMI + comorbidity
  ),
  sp_order  = "size",
  precali   = TRUE,
  control   = pw_solver_control(ftol=1e-6)
)

summary(fit_multi_order_by_size)

## ----multi-ref-as-given-------------------------------------------------------
fit_multi_as_given <- est_pw(
  data      = list(sc, ref1_design, ref2_design),
  p_formula = list(
    ~ agecat + race + education + psa_level + pros_enlarged + comorbidity + diabetes,
    ~ agecat + race + BMI + diabetes + comorbidity
  ),
  sp_order  = "given",
  precali   = TRUE
)

summary(fit_multi_as_given)

## ----pwmean-numeric-----------------------------------------------------------
est_cali_numeric <- pwmean(fit_cali, y = "psa_level")
print(est_cali_numeric)
summary(est_cali_numeric)

## ----pwmean-factor------------------------------------------------------------
est_cali_factor <- pwmean(fit_cali, y = "BMI")
print(est_cali_factor)
summary(est_cali_factor)

## ----pwmean-domain-numeric----------------------------------------------------
est_by_bmi_numeric <- pwmean(fit_cali, y = "psa_level", zcol = "BMI")
print(est_by_bmi_numeric)
summary(est_by_bmi_numeric)

## ----pwmean-domain-factor-----------------------------------------------------
est_by_bmi_factor <- pwmean(fit_cali, y = "agecat", zcol = "BMI")
print(est_by_bmi_factor)
summary(est_by_bmi_factor)

## ----na-setup-----------------------------------------------------------------
sc_with_na <- sc
sc_with_na$agecat[c(1, 2, 4)] <- NA
sc_with_na$psa_level[c(6, 7, 8, 9)] <- NA

sp1_with_na <- sp1
sp1_with_na$education[c(3, 4)] <- NA

ref1_design_with_na <- survey::svydesign(
  ids     = ~psu_sp1,
  strata  = ~strata_sp1,
  weights = ~wts_sp1,
  data    = sp1_with_na,
  nest    = TRUE
)

## ----fit-with-na--------------------------------------------------------------
fit_na_exclude <- est_pw(
  data      = list(sc_with_na, ref1_design_with_na),
  method    = "calibration",
  p_formula = ~ agecat + race + education + comorbidity + BMI + diabetes,
  na.action = stats::na.exclude,
  control   = pw_solver_control(ftol = 1e-6)
)
summary(fit_na_exclude)


## ----na_summary---------------------------------------------------------------
fit_na_exclude$na_summary

## ----na.action----------------------------------------------------------------
na.action(fit_na_exclude)

## ----na-head------------------------------------------------------------------
head(fit_na_exclude$pseudo_weights)
head(fit_na_exclude$sc_updated)

## ----na-2nd-layer-------------------------------------------------------------
result <- pwmean(fit_na_exclude, "psa_level")
summary(result)

## ----na-info-2nd--------------------------------------------------------------
result$na.action
na.action(result)

## ----fit-na-omit--------------------------------------------------------------
fit_na_omit <- est_pw(
  data      = list(sc_with_na, ref1_design_with_na),
  method    = "calibration",
  na.action = stats::na.omit,
  control   = pw_solver_control(ftol = 1e-6)
)

## ----na-omit-na.action--------------------------------------------------------
na.action(fit_na_omit)

## ----na-omit-head-------------------------------------------------------------
head(fit_na_omit$pseudo_weights)
head(fit_na_omit$sc_updated)

## ----solver-control-----------------------------------------------------------
ctrl <- pw_solver_control(
  method   = "Broyden",
  xtol     = 1e-10,
  ftol     = 1e-4,
  maxit    = 30,
  trace    = TRUE
)

fit_ctrl <- est_pw(
  data      = list(sc, ref1_design),
  p_formula = ~ agecat + race + education + comorbidity + BMI,
  method    = "calibration",
  control   = ctrl
)