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

## ----load-package-data--------------------------------------------------------
library(saeproj.multilevel)

data("saeml_modelsvy")
data("saeml_projsvy")

## ----inspect-data-------------------------------------------------------------
dim(saeml_modelsvy)
dim(saeml_projsvy)

head(saeml_modelsvy)
head(saeml_projsvy)

## ----model-survey-domain-count------------------------------------------------
model_domain_count <- table(saeml_modelsvy$kab_kota)

c(
  n_domains = length(model_domain_count),
  min_units_per_domain = min(model_domain_count),
  max_units_per_domain = max(model_domain_count)
)

## ----projection-survey-domain-count-------------------------------------------
proj_domain_count <- table(saeml_projsvy$kab_kota)

c(
  n_domains = length(proj_domain_count),
  min_units_per_domain = min(proj_domain_count),
  max_units_per_domain = max(proj_domain_count)
)

## ----response-variable-check--------------------------------------------------
"Y" %in% names(saeml_modelsvy)

"Y" %in% names(saeml_projsvy)

## ----fit-estimator------------------------------------------------------------
result <- sae_ml_linear(
  formula = Y ~ X1 + X2 + X3 + X4 + Z1 + Z2 + (1 | kab_kota),
  data_model = saeml_modelsvy,
  data_proj = saeml_projsvy,
  domain = "kab_kota",
  cluster_ids = ~1,
  weight = "WEIND",
  strata = "kab_kota",
  summary_function = "mean"
)

## ----estimator-summary--------------------------------------------------------
summary(result)

## ----final-estimates----------------------------------------------------------
head(result$estimates)

## ----estimates-data-frame-----------------------------------------------------
estimates <- as.data.frame(result)

head(estimates)

## ----estimation-components----------------------------------------------------
head(result$estimation_details)

## ----synthetic-component------------------------------------------------------
head(
  result$estimation_details[, c(
    "kab_kota",
    "estimate_synthetic",
    "variance_synthetic"
  )]
)

## ----correction-component-----------------------------------------------------
head(
  result$estimation_details[, c(
    "kab_kota",
    "correction",
    "variance_correction"
  )]
)

## ----model-parameters---------------------------------------------------------
# Fixed-effect estimates
result$model_parameters$fixed_effects

# Random-effect and residual variance components
result$model_parameters$variance_components

# Residual variance
result$model_parameters$residual_variance

## ----random-effects-----------------------------------------------------------
head(result$model_parameters$random_effects$kab_kota)

## ----model-diagnostics--------------------------------------------------------
result$diagnostics

## ----diagnostics-table--------------------------------------------------------
data.frame(
  icc = result$diagnostics$icc,
  singular_fit = result$diagnostics$singular_fit,
  convergence = result$diagnostics$convergence,
  sigma = result$diagnostics$sigma,
  residual_variance = result$diagnostics$residual_variance,
  REML = result$diagnostics$REML,
  AIC = result$diagnostics$AIC,
  BIC = result$diagnostics$BIC
)

## ----fitted-model-------------------------------------------------------------
fit <- result$fitted_model

summary(fit)

## ----residual-plot, fig.width = 7, fig.height = 5-----------------------------
plot(
  fitted(fit),
  resid(fit),
  xlab = "Fitted values",
  ylab = "Residuals",
  main = "Residuals versus Fitted Values"
)

abline(h = 0, lty = 2)

## ----residual-qq-plot, fig.width = 7, fig.height = 5--------------------------
qqnorm(resid(fit))
qqline(resid(fit))

## ----fitted-random-effects----------------------------------------------------
lme4::ranef(fit)

## ----keep-unit, eval = FALSE--------------------------------------------------
# result_unit <- sae_ml_linear(
#   formula = Y ~ X1 + X2 + X3 + X4 + Z1 + Z2 + (1 | kab_kota),
#   data_model = saeml_modelsvy,
#   data_proj = saeml_projsvy,
#   domain = "kab_kota",
#   cluster_ids = ~1,
#   weight = "WEIND",
#   strata = "kab_kota",
#   summary_function = "mean",
#   keep_unit = TRUE
# )
# 
# head(result_unit$unit_projection)
# 
# head(result_unit$unit_model_residual)

## ----direct-estimator, eval = FALSE-------------------------------------------
# result_direct <- sae_ml_linear(
#   formula = Y ~ X1 + X2 + X3 + X4 + Z1 + Z2 + (1 | kab_kota),
#   data_model = saeml_modelsvy,
#   data_proj = saeml_projsvy,
#   domain = "kab_kota",
#   cluster_ids = ~1,
#   weight = "WEIND",
#   strata = "kab_kota",
#   summary_function = "mean",
#   return_direct = TRUE
# )
# 
# head(result_direct$direct_estimator)

## ----multiple-domain-variables, eval = FALSE----------------------------------
# result_multi <- sae_ml_linear(
#   formula = Y ~ X1 + X2 + X3 + X4 + Z1 + Z2 + (1 | kab_kota),
#   data_model = saeml_modelsvy,
#   data_proj = saeml_projsvy,
#   domain = c("prov", "kab_kota"),
#   cluster_ids = ~1,
#   weight = "WEIND",
#   strata = "kab_kota",
#   summary_function = "mean"
# )
# 
# head(result_multi$estimates)

## ----clustered-survey-design, eval = FALSE------------------------------------
# result_clustered <- sae_ml_linear(
#   formula = Y ~ X1 + X2 + X3 + X4 + Z1 + Z2 + (1 | kab_kota),
#   data_model = data_model,
#   data_proj = data_proj,
#   domain = "kab_kota",
#   cluster_ids = "psu_id",
#   weight = "survey_weight",
#   strata = "stratum",
#   summary_function = "mean",
#   nest = TRUE
# )

## ----output-structure---------------------------------------------------------
names(result)