## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4)
has_sf  <- requireNamespace("sf", quietly = TRUE)
has_xgb <- requireNamespace("xgboost", quietly = TRUE)

## -----------------------------------------------------------------------------
library(SurveyNCD)
library(survey)

## -----------------------------------------------------------------------------
d <- ncd_survey
for (v in c("htn", "dm", "arth", "asthma", "angina", "stroke")) {
  d[[v]] <- suppressWarnings(recode_binary(d[[paste0(v, "_raw")]], yes = 1, no = 2))
}
d$dep <- suppressWarnings(recode_binary(d$dep_raw, yes = "Yes", no = "No"))
conds  <- c("htn", "dm", "arth", "asthma", "angina", "stroke", "dep")
labels <- c(htn = "Hypertension", dm = "Diabetes", arth = "Arthritis",
            asthma = "Asthma", angina = "Angina", stroke = "Stroke",
            dep = "Depression")

scored <- multimorbidity_index(d, conditions = conds, na_action = "na")
scored$age_group <- cut(scored$age, c(18, 30, 45, 60, 70), right = FALSE,
                        labels = c("18-29", "30-44", "45-59", "60-69"))
table(scored$mm_category, useNA = "ifany")

## -----------------------------------------------------------------------------
des <- steps_design(scored, weight = "wt")

## -----------------------------------------------------------------------------
ncd_table(des, conds, by = "sex", labels = labels)

## -----------------------------------------------------------------------------
prev <- survey_prevalence(des, c("htn", "dm"), by = c("residence", "sex"))
prev
ggplot2::autoplot(prev)

## -----------------------------------------------------------------------------
mm_prevalence(des)
mm_prevalence(des, by = c("residence", "sex"))

## -----------------------------------------------------------------------------
std <- who_standard_population(breaks = c(18, 30, 45, 60, 70))
std
mm_standardize(des, c("mm_multimorbid", "htn"), age_group = "age_group",
               standard = std, by = "residence")

## ----fig.alt = "Bar chart of the prevalence of the most common pairs of conditions."----
pairs <- mm_patterns(des, conds, top = 8)
pairs
ggplot2::autoplot(pairs)

## -----------------------------------------------------------------------------
mm_bounds(des, conds)

## -----------------------------------------------------------------------------
survey_concentration_index(des, outcome = dm, wealth = wealth_score)

## -----------------------------------------------------------------------------
survey_concentration_index(des, outcome = dm, wealth = wealth_score, type = "erreygers")
survey_concentration_index(des, outcome = dep, wealth = wealth_score, type = "erreygers",
                           by = "residence")

## ----fig.alt = "Concentration curves for diabetes in urban and rural areas lying below the diagonal."----
ggplot2::autoplot(survey_concentration_curve(des, dm, wealth_score, by = "residence"))
survey_slope_index(des, dm, wealth_score)

## -----------------------------------------------------------------------------
ch <- child_anthro
ch$stunting <- who_anthro_score(ch$hw70, indicator = "stunting")
ch$stunted  <- as.numeric(ch$stunting %in% c("Severe stunting", "Moderate stunting"))
ch$stunted[is.na(ch$stunting)] <- NA
cdes <- steps_design(ch, weight = "wt")
survey_prevalence(cdes, "stunted", by = "wealth_quintile")
survey_concentration_index(cdes, outcome = stunted, wealth = wealth_score, type = "wagstaff")

## ----map, eval = has_sf, fig.alt = "Choropleth of multimorbidity prevalence in six fictitious regions."----
reg <- mm_prevalence(des, by = "region")
grid <- sf::st_make_grid(sf::st_as_sfc(sf::st_bbox(c(xmin = 0, ymin = 0, xmax = 3, ymax = 2))),
                         n = c(3, 2))
regions_sf <- sf::st_sf(region = paste("Region", LETTERS[c(4, 5, 6, 1, 2, 3)]),
                        geometry = grid)
survey_map_indicator(as.data.frame(reg), regions_sf, join_by = "region",
                     fill_var = "prevalence", palette = "cividis")

## ----ml, eval = has_xgb, fig.alt = "SHAP summary plot ranking age, BMI and wealth as the main contributors."----
scored$female <- as.numeric(scored$sex == "Female")
scored$urban  <- as.numeric(scored$residence == "Urban")
des_ml <- steps_design(scored, weight = "wt")
fit <- survey_xgboost(des_ml, mm_n_conditions ~ age + bmi + wealth_score + female + urban,
                      params = list(objective = "reg:squarederror", max_depth = 3,
                                    eta = 0.05, nthread = 1),
                      nrounds = 150)
plot_shap_summary(survey_shap(fit), fit$X)