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

## ----message=FALSE------------------------------------------------------------
library(survey)
library(svymargins)
library(dplyr)
library(svyVGAM) # for multinomial logistic regression models

## -----------------------------------------------------------------------------
n <- 1000
# Generate data:
set.seed(1234)
d <- data.frame(sex=factor(sample(c("M", "F"), n, replace=TRUE)),
                education=factor(sample(c("low", "middle", "high"), n, replace=TRUE)))
d <- d |> mutate(age=runif(n, 0, 40) + as.numeric(education) * 20,
                 y=rnorm(n, sd=5) + as.numeric(education) + 0.05 * age)
summary(d)

## ----eval=FALSE, echo=FALSE---------------------------------------------------
# cmd <- "svyset [pw=one]
# svy: reg y i.education age i.sex
# margins sex education, grand // vce(unconditional)
# margins education, at(age=(40 50 60 70)) vce(unconditional)"
# # options("RStata.StataPath" =here::here("/usr/bin/stata-se")); options("RStata.StataVersion" = 16); stata(cmd, data.in = d |> mutate(one=1))

## -----------------------------------------------------------------------------
# Create survey design:
my.svy <- svydesign(~ 1, weights=~ 1, data=d)
# Run regression analysis:
res <- svyglm(y ~ education + age + sex, design=my.svy)
summary(res)

## -----------------------------------------------------------------------------
# Define margins as a named list:
target.l <- list(null=list(),
                 sex=list("sex"),
                 educ=list("education"),
                 age=list(age=seq(40,70,10)),
                 educ_age=list("education", age=seq(40,70,10)))

## -----------------------------------------------------------------------------
# Calculate predictive margins:
svymargins(res, groupfactor=target.l)
# Get the output table containing the covariate information from "groups" attribute:
attr(svymargins(res, groupfactor=target.l), "groups")

## -----------------------------------------------------------------------------
marg <- svymargins(res, groupfactor=target.l)
cbind(marg, confint(marg))
marg <- svymargins(res, groupfactor=target.l)
cbind(attr(marg, "groups"), confint(marg))

## -----------------------------------------------------------------------------
marg <- svymargins(res, groupfactor=target.l)

# Difference of adjusted means between high education and population:
svycontrast(marg, quote(`educ,high` - `null`))

# Difference of adjusted means between high and low education:
svycontrast(marg, quote(`educ,high` - `educ,low`))
# ... and confidence interval:
svycontrast(marg, quote(`educ,high` - `educ,low`)) |> confint()

# Ratio of adjusted means between high and low education:
svycontrast(marg, quote(`educ,high` / `educ,low`))

## -----------------------------------------------------------------------------
n <- 10000
# Generate data:
set.seed(1234)
d <- data.frame(sex=factor(sample(c("M", "F"), n, replace=TRUE)),
                education=factor(sample(c("low", "middle", "high"), n, replace=TRUE)))
d <- d |>
  mutate(age=runif(n, 0, 40) + as.numeric(education) * 20,
         pr1=1,
         pr2=exp(-1 + 0.5 * (as.numeric(education)-1) + 0.02 * age),
         pr3=exp(1 + -0.5 * (as.numeric(education)-1) + 0.02 * age)) |> 
  mutate(across(matches("pr[0-9]"), ~ .x / (pr1 + pr2 + pr3))) |> 
  rowwise() |>
  mutate(y=which(rmultinom(1, 1, c(pr1, pr2, pr3))[,1] == 1)) |> 
  ungroup() |>
  mutate(y=factor(y, levels=1:3, labels=LETTERS[1:3]))
summary(d)
# Create survey design:
my.svy <- svydesign(~ 1, weights=~ 1, data=d)
# Run regression analysis:
res <- svy_vglm(y ~ education + age + sex,
                family=multinomial(refLevel=1), design=my.svy)
summary(res)

## -----------------------------------------------------------------------------
# Define margins as a named list:
target.l <- list(null=list(),
                 educ=list("education"),
                 age=list(age=seq(40,70,10)),
                 educ_age=list("education", age=seq(40,70,10)))

# Calculate predictive margins
# for the 1st outcome level "A":
svymargins(res, groupfactor=target.l, y.lev=1)

# Get the output table containing the covariate information from "groups" attribute
# for the 2nd outcome level "B":
attr(svymargins(res, groupfactor=target.l, y.lev="B"), "groups")

## -----------------------------------------------------------------------------
marg <- svymargins(res, groupfactor=target.l, y.lev=1)
cbind(marg, confint(marg))
marg <- svymargins(res, groupfactor=target.l, y.lev="B")
cbind(attr(marg, "groups"), confint(marg))

## -----------------------------------------------------------------------------
marg <- svymargins(res, groupfactor=target.l, y.lev=1)
cbind(marg, confint(marg) |> asymm_ci())
marg <- svymargins(res, groupfactor=target.l, y.lev="B")
cbind(attr(marg, "groups"), confint(marg) |> asymm_ci())

## -----------------------------------------------------------------------------
marg <- svymargins(res, groupfactor=target.l, y.lev=2)

# Adjusted risk (or prevalence) difference between high education and population:
svycontrast(marg, quote(`educ,high` - `null`))

# Adjusted risk (or prevalence) difference between high and low education:
svycontrast(marg, quote(`educ,high` - `educ,low`))
# ... and confidence interval:
svycontrast(marg, quote(`educ,high` - `educ,low`)) |> confint()

# Adjusted risk (or prevalence) ratio between high and low education:
svycontrast(marg, quote(`educ,high` / `educ,low`))

# Population attributable fraction (PAF) of high education:
svycontrast(marg, quote(1 - `educ,high` / `null`))