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

## ----setup--------------------------------------------------------------------
library(ggstratify)

## ----eval = FALSE-------------------------------------------------------------
# dat <- transform(
#   dat,
#   sex = factor(sex, levels = c("Male", "Female")),
#   severity = factor(severity, levels = c("Mild", "Moderate", "Severe")),
#   age = as.numeric(age)
# )

## ----eval = FALSE-------------------------------------------------------------
# cohort <- read.csv("cohort.csv")
# str(cohort)
# ggstratify(cohort)

## ----eval = FALSE-------------------------------------------------------------
# ggstratify(epi_cohort)

## -----------------------------------------------------------------------------
str(epi_cohort)

## ----eval = FALSE-------------------------------------------------------------
# stat_summary(fun.data = mean_se, geom = "pointrange")

## ----eval = FALSE-------------------------------------------------------------
# mean_ci <- function(x, conf = 0.95) {
#   x <- x[!is.na(x)]
#   n <- length(x)
#   m <- if (n) mean(x) else NA_real_
#   if (n < 2L) return(data.frame(y = m, ymin = NA_real_, ymax = NA_real_))
#   half <- qt(1 - (1 - conf) / 2, n - 1L) * sd(x) / sqrt(n)
#   data.frame(y = m, ymin = m - half, ymax = m + half)
# }

## ----eval = FALSE-------------------------------------------------------------
# ggplot(d, aes(x = severity, y = as.integer(outcome == "Died"))) +
#   stat_summary(fun.data = prop_ci_wilson, fun.args = list(conf = 0.95),
#                geom = "pointrange") +
#   labs(y = "Proportion outcome = Died")

## ----eval = FALSE-------------------------------------------------------------
# prop_ci_exact <- function(x, conf = 0.95) {
#   x <- as.numeric(x); x <- x[!is.na(x)]
#   ci <- binom.test(sum(x), length(x), conf.level = conf)$conf.int
#   data.frame(y = mean(x), ymin = ci[1], ymax = ci[2])
# }

## ----eval = FALSE-------------------------------------------------------------
# prop_ci_wilson <- function(x, conf = 0.95) {
#   x <- as.numeric(x); x <- x[!is.na(x)]
#   n <- length(x); p <- mean(x)
#   z <- qnorm(1 - (1 - conf) / 2)
#   denom <- 1 + z^2 / n
#   centre <- (p + z^2 / (2 * n)) / denom
#   half <- z * sqrt(p * (1 - p) / n + z^2 / (4 * n^2)) / denom
#   data.frame(y = p, ymin = max(0, centre - half), ymax = min(1, centre + half))
# }

## ----eval = FALSE-------------------------------------------------------------
# ggplot(d, aes(x = admit_year, y = los_days, colour = sex)) +
#   stat_summary(aes(group = sex), fun.data = mean_se, geom = "line") +
#   stat_summary(fun.data = mean_se, geom = "pointrange")

## ----eval = FALSE-------------------------------------------------------------
# ggplot(d, aes(x = age, weight = survey_weight)) +
#   geom_histogram()

## ----eval = FALSE-------------------------------------------------------------
# dt[, .svy_row := .I]
# des <- survey::svydesign(ids = ~ward, strata = ~site, weights = ~svy_weight,
#                          nest = TRUE, data = as.data.frame(dt))
# 
# survey::svymean(~crp, des[rows, ])                      # a mean and its SE
# survey::svyciprop(~died, des[rows, ], method = "beta")  # a proportion
# survey::svykm(survival::Surv(fu_days, death) ~ 1, des[rows, ], se = TRUE)

## ----eval = FALSE-------------------------------------------------------------
# theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

## ----eval = FALSE-------------------------------------------------------------
# dt[, age_cat := cut(age, breaks = c(-Inf, 65, Inf))]

## ----eval = FALSE-------------------------------------------------------------
# dt[, crp_missing := factor(is.na(crp), levels = c(FALSE, TRUE),
#                            labels = c("Observed", "Missing"))]

## ----eval = FALSE-------------------------------------------------------------
# dt[, admit_month := as.IDate(cut(as.IDate(admit_date), breaks = "month"))]

## ----eval = FALSE-------------------------------------------------------------
# dt[, admit_month := factor(month.abb[month(admit_date)], levels = month.abb)]

## ----eval = FALSE-------------------------------------------------------------
# dt[, admit_season := factor(
#   c("Spring", "Summer", "Fall", "Winter")[
#     ((month(admit_date) - 3L) %% 12L) %/% 3L + 1L],
#   levels = c("Spring", "Summer", "Fall", "Winter"))]

## ----eval = FALSE-------------------------------------------------------------
# ggplot(d, aes(x = admit_date, y = los_days)) +
#   geom_line(alpha = 0.6) +
#   scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
#   labs(x = "Year") +
#   theme_bw()

## ----eval = FALSE-------------------------------------------------------------
# library(data.table)
# library(ggplot2)
# 
# dt <- as.data.table(epi_cohort)
# 
# d <- dt[sex == "Male"]
# 
# p <- ggplot(d, aes(x = treatment, y = los_days)) +
#   geom_boxplot() +
#   theme_bw()
# 
# p