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

## ----data---------------------------------------------------------------------
data(birthwt, package = "MASS")

# Recode predictors
birthwt$smoke <- factor(birthwt$smoke, labels = c("Non-smoker", "Smoker"))
birthwt$race  <- factor(birthwt$race,
                        labels = c("White", "Black", "Other"))
birthwt$ht    <- factor(birthwt$ht,   labels = c("No", "Yes"))
birthwt$ui    <- factor(birthwt$ui,   labels = c("No", "Yes"))

# Outcome prevalence
mean(birthwt$low)   # 31 % — common outcome, OR is a poor approximation

## ----glm-fit------------------------------------------------------------------
fit_glm <- glm(low ~ smoke + race + age + lwt + ht + ui,
               family = binomial, data = birthwt)

## ----delta-cond---------------------------------------------------------------
res_cond <- prLogisticDelta(fit_glm, standardisation = "conditional")
res_cond

## ----delta-marg---------------------------------------------------------------
res_marg <- prLogisticDelta(fit_glm, standardisation = "marginal")
res_marg

## ----ref-values---------------------------------------------------------------
prLogisticDelta(fit_glm,
                standardisation = "conditional",
                ref_values      = list(age = 25, lwt = 55))

## ----forest-plot, fig.cap = "Forest plot: conditional PR estimates with 95% CI"----
plot(res_cond, main = "Prevalence Ratios — conditional (birthwt)")

## ----bootstrap, cache = TRUE--------------------------------------------------
set.seed(2024)
res_boot_c <- prLogisticBootCond(fit_glm, data = birthwt, R = 999)
res_boot_c

## ----confint-boot-------------------------------------------------------------
# Percentile intervals
confint(res_boot_c, type = "percentile")

## ----glmer, eval = requireNamespace("lme4", quietly = TRUE)-------------------
library(lme4)

# Treat race as a clustering variable (illustrative)
fit_glmer <- glmer(low ~ smoke + age + lwt + ht + (1 | race),
                   family = binomial, data = birthwt)

prLogisticDelta(fit_glmer, standardisation = "marginal")

## ----gee, eval = requireNamespace("geepack", quietly = TRUE)------------------
library(geepack)
data(ohio, package = "geepack")

# Respiratory symptoms in children (4 repeated measures per child)
fit_gee <- geeglm(resp ~ smoke + age,
                  family = binomial,
                  id     = id,
                  corstr = "exchangeable",
                  data   = ohio)

prLogisticGEE(fit_gee)

## ----survey, eval = requireNamespace("survey", quietly = TRUE)----------------
library(survey)
data(api, package = "survey")

apiclus2$target_met <- as.numeric(apiclus2$sch.wide == "Yes")

# Two-stage cluster sample
dclus2 <- svydesign(id   = ~dnum + snum,
                    fpc  = ~fpc1 + fpc2,
                    data = apiclus2)

fit_svy <- svyglm(target_met ~ meals + stype,
                  design = dclus2,
                  family = quasibinomial)

prLogisticSurvey(fit_svy, standardisation = "conditional")

## ----or-vs-pr-----------------------------------------------------------------
OR  <- exp(coef(fit_glm)["smokeSmoker"])
PR  <- coef(res_cond)["smokeSmoker"]

data.frame(
  Measure  = c("Odds Ratio (logistic)", "Prevalence Ratio (conditional)",
               "Prevalence Ratio (marginal)"),
  Estimate = round(c(OR, PR, coef(res_marg)["smokeSmoker"]), 3)
)

## ----session------------------------------------------------------------------
sessionInfo()