knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
set.seed(1)
library(flexsynth)

n <- 800
real <- data.frame(
  id  = seq_len(n),
  age = round(rnorm(n, 62, 11)),
  sex = sample(c("F", "M"), n, replace = TRUE, prob = c(0.45, 0.55))
)
real$sbp <- round(0.6 * real$age + ifelse(real$sex == "M", 5, 0) +
                    rnorm(n, 90, 10))

res <- synth(real, structure = ~ id, m = 10, seed = 1)
res

pooled <- synth_glm(res, sbp ~ age + sex)
pooled

summary(lm(sbp ~ age + sex, data = real))$coefficients[, 1:2]

# A proportion and its variance, supplied directly.
pool_synth(res, function(d) {
  p <- mean(d$sex == "M")
  list(estimate = c(prop_male = p),
       variance = c(prop_male = p * (1 - p) / nrow(d)))
})

compare_estimates(real, res, function(d) lm(sbp ~ age + sex, data = d))

m <- 800
clinic <- data.frame(
  region = sample(c("N", "S", "E", "W"), m, replace = TRUE),
  agecat = sample(c("40s", "50s", "60s", "70s"), m, replace = TRUE)
)
# hypertension risk rises steeply with age, with a small regional bump.
prob_htn <- c("40s" = 0.12, "50s" = 0.40, "60s" = 0.75, "70s" = 0.93)
reg_bump <- c("N" = 0, "S" = 0, "E" = 0.03, "W" = 0.05)
p <- pmin(0.98, prob_htn[clinic$agecat] + reg_bump[clinic$region])
clinic$hypertension <- ifelse(runif(m) < p, "yes", "no")
clinic$id <- seq_len(m)

cres <- synth(clinic, ~ id, m = 1, seed = 1)

disclosure_risk(clinic, cres,
                quasi  = c("region", "agecat", "hypertension"),
                target = "hypertension",
                seed   = 1)