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

## ----setup--------------------------------------------------------------------
library(NMAR)

## -----------------------------------------------------------------------------
set.seed(1108)
head(riddles_case1)

## -----------------------------------------------------------------------------
y_original <- riddles_case1$y_true # to compare with true values

## -----------------------------------------------------------------------------
# Exptilt engine configuration
exptilt_config <- exptilt_engine(
  y_dens = 'normal', # Assume f(y|x, delta=1) is Normal
  family = 'logit', # Logit link for response probability
  variance_method = 'bootstrap',
  bootstrap_reps = 5,
  control = list(maxit = 10), # Solver control parameters
  stopping_threshold = 0.01, # Convergence threshold
  standardize = FALSE
)

formula = y ~x
res <- nmar(formula = formula, data = riddles_case1, engine = exptilt_config, trace_level = 0)

## -----------------------------------------------------------------------------
print(res)

## -----------------------------------------------------------------------------
coef(res)

## -----------------------------------------------------------------------------
y_true_mean <- mean(riddles_case1$y_true)
est_mean <- as.numeric(res$y_hat)
se <- res$se

cat('Estimated (NMAR) Mean:', sprintf('%.4f', est_mean),
    ' (SE:', sprintf('%.4f', se), ')\n')
cat('True Y Mean:          ', sprintf('%.4f', y_true_mean), '\n')
cat('Naive (MCAR) Mean:     ', sprintf('%.4f', mean(riddles_case1$y, na.rm = TRUE)), '\n')

## -----------------------------------------------------------------------------
# if survey is installed
if (requireNamespace("survey", quietly = TRUE)) {

  library('survey')
  surv_test_weights <- runif(length(riddles_case1$y), 0.9, 1.1)
  surv_test_data <- svydesign(ids = ~1, data = riddles_case1, weights = ~surv_test_weights)

  exptilt_config <- exptilt_engine(
     standardize = FALSE,
     on_failure = "error",
     bootstrap_reps = 5,
     supress_warnings = FALSE,
     control = list(maxit = 15, method = "Newton"),
     family = "logit",
     y_dens = "normal",
     variance_method = "bootstrap",
     stopping_threshold = 1
)


  formula = y ~ x
  res <- nmar(formula = formula, data = surv_test_data, engine = exptilt_config, trace_level = 3)
  print(res)

}