## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5
)

## ----load---------------------------------------------------------------------
library(MAIHDA)
data("maihda_health_data")

# A two-level outcome: obese (Yes) vs. not (No)
table(maihda_health_data$Obese)

## ----autodetect---------------------------------------------------------------
model_null <- fit_maihda(
  Obese ~ 1 + (1 | Gender:Race:Education),
  data = maihda_health_data
)

## ----explicit, eval=FALSE-----------------------------------------------------
# model_null <- fit_maihda(
#   Obese ~ 1 + (1 | Gender:Race:Education),
#   data   = maihda_health_data,
#   family = "binomial"
# )

## ----summary-null-------------------------------------------------------------
summary_null <- summary(model_null)
print(summary_null)

## ----vpc-response-------------------------------------------------------------
maihda_vpc_response(model_null, seed = 1)

## ----adjusted-----------------------------------------------------------------
health_complete <- maihda_health_data[complete.cases(
  maihda_health_data[, c("Obese", "Age", "Gender", "Race", "Education")]
), ]

model_null2 <- fit_maihda(
  Obese ~ 1 + (1 | Gender:Race:Education),
  data = health_complete, family = "binomial"
)

# Model 2: adjust for an individual-level covariate (age)
model_adj <- fit_maihda(
  Obese ~ Age + (1 | Gender:Race:Education),
  data = health_complete, family = "binomial"
)

pcv <- calculate_pcv(model_null2, model_adj)
print(pcv)

## ----adjusted-reml------------------------------------------------------------
calculate_pcv(model_null2, model_adj, estimation = "REML")

## ----interactions-------------------------------------------------------------
model_int <- fit_maihda(
  Obese ~ Gender + Race + Education + (1 | Gender:Race:Education),
  data = health_complete, family = "binomial"
)

maihda_interactions(model_int)

## ----scale-dependence---------------------------------------------------------
p <- function(a, b) plogis(-2 + 0.7 * a + 0.7 * b)

c(`risk added by B, at A = 0` = p(0, 1) - p(0, 0),
  `risk added by B, at A = 1` = p(1, 1) - p(1, 0))

## ----scale-dependence-zero----------------------------------------------------
rd_int <- function(b0, a, b) {
  (plogis(b0 + a + b) - plogis(b0 + a)) - (plogis(b0 + b) - plogis(b0))
}

vapply(c(`-2.0` = -2, `-0.7` = -0.7, `0.0` = 0), rd_int, numeric(1), a = 0.7, b = 0.7)

## ----interactions-response----------------------------------------------------
maihda_interactions(model_int, scale = "response")

## ----interactions-rope--------------------------------------------------------
table(maihda_interactions(model_int, scale = "response", rope = 0.02)$decision)

## ----da-----------------------------------------------------------------------
da_null <- maihda_discriminatory_accuracy(model_null2)
da_adj  <- maihda_discriminatory_accuracy(model_adj)

da_null
da_adj

## ----auc-direct---------------------------------------------------------------
prob_null <- predict_maihda(model_null2, type = "individual", scale = "response")
y_obs     <- as.numeric(lme4::getME(model_null2$model, "y"))

maihda_auc(prob_null, y_obs)

## ----plot-predicted-----------------------------------------------------------
# Predicted probabilities per stratum with intervals
plot(model_adj, type = "predicted")

## ----plot-vpc-----------------------------------------------------------------
# Latent-scale variance partition
plot(model_adj, type = "vpc")

## ----plot-prediction-deviation, warning = FALSE-------------------------------
# For binomial fits the dashboard highlights the largest absolute
# deviance residuals rather than raw deviations from the mean.
plot(model_adj, type = "prediction_deviation")