### R code from vignette source 'prev.rnw'

###################################################
### code chunk number 1: prev.rnw:3-10 (eval = FALSE)
###################################################
## options(width = 90,
##         show.signif.stars = FALSE)
## par(mar = c(3, 3, 1, 1),
##     mgp = c(3, 1, 0) / 1.6,
##     las = 1,
##    lend = "butt",
##     bty = "n")


###################################################
### code chunk number 2: prev.rnw:12-20
###################################################
options(width = 90,
        show.signif.stars = FALSE,
        SweaveHooks=list(fig=function()
                         par(mar = c(3, 3, 1, 1),
                             mgp = c(3, 1, 0) / 1.6,
                             las = 1,
                            lend = "butt",
                             bty = "n")))


###################################################
### code chunk number 3: prev.rnw:42-46
###################################################
library(Epi)
library(popEpi)
library(survival)
library(tidyverse)


###################################################
### code chunk number 4: prev.rnw:48-50
###################################################
clear()
setwd("C:/Bendix/teach/AdvCoh/courses/IDEG2025/pracs")


###################################################
### code chunk number 5: prev.rnw:52-57
###################################################
vrs <-
data.frame(R = do.call(paste, c(R.Version()[c("major","minor")], sep =".")),
         Epi = packageVersion("Epi"),
      popEpi = packageVersion("popEpi"))
print(vrs, row.names = FALSE)


###################################################
### code chunk number 6: prev.rnw:67-68 (eval = FALSE)
###################################################
## nhis <- read.csv("../data/NHIS_IDEG.csv", header = TRUE)


###################################################
### code chunk number 7: prev.rnw:70-73
###################################################
nhis <- read.csv(
"https://bendixcarstensen.com/AdvCoh/courses/IDEG2025/data/NHIS_IDEG.csv",
header = TRUE)


###################################################
### code chunk number 8: prev.rnw:75-80
###################################################
str(nhis)
newn <- tolower(gsub("_A", "", names(nhis)))
cbind(names(nhis), newn)
names(nhis) <- newn
str(nhis)


###################################################
### code chunk number 9: prev.rnw:91-96
###################################################
nhis <- mutate(nhis, dibev = factor(dibev, labels = c("Y","N","R","U")),
                   dibtype = factor(dibtype, labels = c("T1","T2","O","O","O")),
                       agr = cut(agep, seq(0, 100, 10), right = FALSE),
                       sex = factor(sex,  labels = c("M","W","U","U")))
str(nhis)


###################################################
### code chunk number 10: prev.rnw:115-117
###################################################
(tb <- with(nhis, table(dibev, exclude = NULL)))
tb["Y"] / (tb["Y"] + tb["N"]) * 100


###################################################
### code chunk number 11: prev.rnw:136-140
###################################################
with(nhis, table(Age = agr,
            Diabetes = dibev,
             exclude = NULL)) |> addmargins() -> diab
diab


###################################################
### code chunk number 12: prev.rnw:152-154
###################################################
(diab <- addmargins(diab[, 1:2], 2))
cbind(round(diab[,"Y"] / diab[,"Sum"] * 100, 1))


###################################################
### code chunk number 13: prev.rnw:167-169
###################################################
with(nhis, table(agr, dibtype, exclude = NULL)) |> addmargins() -> dtyp
dtyp


###################################################
### code chunk number 14: prev.rnw:180-181
###################################################
round(100 * dtyp[, 1:2] / dtyp[, "Sum"], 1)


###################################################
### code chunk number 15: prev.rnw:202-204
###################################################
nh <- subset(nhis, dibev %in% c("Y", "N") & sex %in% c("M", "W"))
str(nh)


###################################################
### code chunk number 16: both
###################################################
ma <- glm((dibev == "Y") ~ Ns(agep, knots = seq(30, 90,, 4)),
          family = binomial,
          data = nh)
da <- data.frame(agep = 10:95)
head(pa <- ci.pred(ma, da) * 100)
matshade(da[,"agep"], pa, plot = TRUE, lwd = 3, ylim = c(0, 25), yaxs = "i")


###################################################
### code chunk number 17: prev.rnw:249-256
###################################################
Ma <- glm((dibev == "Y") ~ Ns(agep, knots = seq(30, 90,, 4)),
          family = binomial,
          data = subset(nh, sex == "M"))
Wa <- update(Ma, data = subset(nh, sex == "W"))
pM <- ci.pred(Ma, da)
pW <- ci.pred(Wa, da)
mw <- ci.ratio(pM, pW)


###################################################
### code chunk number 18: m-w
###################################################
par(mar = c(3,3,1,3))
matshade(da$agep, cbind(pM, pW) * 100, col = c("blue", "red"),
         plot = TRUE, lwd = 3, ylim = c(0, 25), yaxs = "i",
         xlab = "Age (years)", ylab = "Prevalence of diabetes (%)")
axis(side = 1, at = seq(15, 95, 5), labels = NA, tcl = -0.3)
axis(side = 1, at = seq(10, 90,10), labels = NA, tcl = -0.5)
#
matshade(da$agep, mw * 5, lwd =3)
lines(c(10,100), c(5,5))
axis(side = 4, at = c(1, 1.5, 2) * 5, labels = c(1, 1.5, 2))
axis(side = 4, at = seq(0.8, 2, 0.1) * 5, labels = NA, tcl = -0.3)