Paper 1, Arm 1: Lifetime abstainers as the reference group
EPI 601, Paper 1
September 14, 2026
What this document is
This is arm 1. It estimates the odds of a higher CKM stage for occasional, light, moderate, and heavy drinking compared with people who have never drunk alcohol. Former drinkers are shown as their own row and are never part of the reference. Parts 1 to 5 and Table 1 below are the shared opening, identical in all three arm documents; Parts 6 to 9 are this arm’s.
Paper 1. Who counts as a non-drinker? Reference-group construction and the alcohol-CKM gradient. NHANES 2011 to March 2020, adults 20 years and older in the fasting subsample, CKM syndrome stage 0 to 4. Target journal Addiction.
What the shared opening builds. The survey design on
the whole file, the four exclusions that define the analytic sample, the
six-category drinking variable drinkcat that all three arms
are cut from, the shared covariate set, and the cell check that decides
whether the outcome is the five-stage ckmstage or the
three-level ckmtier.
What it does not do. The opening fits no model and makes no claim. Parts 6 to 9 do that, and they live in the arm documents.
The covariate set. Age, sex, race and ethnicity,
education, family income to poverty ratio, smoking, self-rated health,
and physical activity, which is the shared set
specs/00-cover.md lists, with Paper 2 adjusting for the
same variables built the same way. One item on the cover’s list is
missing from every model, and it is missing from the file rather than by
choice: diet quality (HEI-2015) needs the dietary recall tables, which
this build does not carry. No component of the outcome is ever a
covariate.
Two definitions, not one. Lifetime abstention is
ALQ101 = 2 and ALQ110 = 2 in 2011 to 2016 and
ALQ111 = 2 in 2017 to March 2020. Those are different
questions, so every result in this paper is reported by instrument era
before it is pooled. era is 1 for 2011 to 2016 and 2 for
2017 to March 2020.
The pre-specified rule fired, so the primary outcome is
ckmtier. On the analytic file this repository
distributes, the thinnest cell of the exposure-by-stage table holds 23
people (heavy drinkers at stage 3), which is below 30, so the rule
written in specs/00-cover.md and in
code/README.md Part 5 collapses the outcome to the three
risk tiers: 1 for stages 0 or 1, 2 for stage 2, 3 for stages 3 or 4.
Every model in every arm runs on ckmtier. Each arm also
reports the same models on the five-stage ckmstage as a
sensitivity analysis, in a table-N-ckmstage file, shown
only because the rule collapsed the primary outcome. Part 5 prints the
table and the verdict for the file you actually read, which is what
governs.
Part 1: Setup and read
Packages
The two paths
Research question for this block: where does the analytic file live, and where does the output go?
# The analytic file. Change "data" to file.path("data", "sample") for a smoke test on the
# 300-row teaching sample, whose numbers are wrong on purpose.
data_file <- here::here("data", "nhanes_ckm_2011_2020.rds")
# Where output is written. The instructor's run writes to output/tentative/; your own run
# writes to here::here("arms", "<yourname>", "output").
out_dir <- here::here("output", "tentative", "arm-01")
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE) # make the directory if it is not there yetRead the file
Research question for this block: did the file this document reads have the rows and columns the codebook describes?
dat <- readRDS("C:\\Users\\safwa\\OneDrive - University at Albany - SUNY\\Fall 2026\\EPI 601\\Paper\\nhanes_ckm_2011_2020.rds")
# the built analytic file, one row per adult 20 and older
dim(dat) # 26,280 rows and 137 columns in the real file## [1] 26280 137##
## 0
## 26280##
## 1 2
## 17048 9232Part 2: Survey design
Research question for this block: what design produced these rows, and is it declared before anything is dropped?
The design is declared on the whole file. Rows are
removed in Part 3 with subset() on the design object, so
variance estimation still sees every stratum and every primary sampling
unit. Never filter the data frame first.
des <- svydesign( # the NHANES design for the fasting subsample
ids = ~SDMVPSU, # masked variance pseudo-PSU, nested in the stratum
strata = ~SDMVSTRA, # masked variance pseudo-stratum
weights = ~wtsafcmb, # combined fasting weight, WTSAF2YR x 2/9.2 plus WTSAFPRP x 3.2/9.2
nest = TRUE, # PSU numbers repeat across strata, so they are nested
data = dat # the whole file, not a subset
)
des # print the design so the weights and nesting are on the record## Stratified 1 - level Cluster Sampling design (with replacement)
## With (140) clusters.
## svydesign(ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~wtsafcmb,
## nest = TRUE, data = dat)## [1] 72Part 3: Analytic sample
Research question for this block: which rows does
Paper 1 analyze, and does the count after each exclusion match
build/build-log.txt?
n_read <- nrow(des) # rows as read: 26,280, every adult 20 and older
n_read # printed, and kept for the run log in Part 9## [1] 26280des_sub <- subset(des, RIDAGEYR >= 20) # aged 20 years and older; the file already is, so this changes nothing
n_age <- nrow(des_sub) # build log: 26,280
n_age # printed## [1] 26280des_sub <- subset(des_sub, is.na(RIDEXPRG) | RIDEXPRG != 1) # not pregnant at examination; RIDEXPRG is asked of women 20 to 44 only
n_preg <- nrow(des_sub) # build log: 26,001
n_preg # printed## [1] 26001des_sub <- subset(des_sub, fastflag == 1) # in the fasting subsample, that is wtsafcmb greater than 0
n_fast <- nrow(des_sub) # build log: 10,679
n_fast # printed## [1] 10679des_sub <- subset(des_sub, insamp == 1) # CKM stage could be determined
n_samp <- nrow(des_sub) # build log: 10,585, the analytic sample of Papers 1 and 2
n_samp # printed## [1] 10585## [1] 72Part 4: Recodes
One variable per update() call on the design object, so
the design and the data can never drift apart. Each block opens with the
codebook line it comes from. Codes that mean refused or don’t know are
set to missing on their own line before the variable is used.
Alcohol items: refused and don’t know to missing
Research question for this block: which alcohol answers are real answers?
# ALQ101, ALQ_G/_H/_I: had at least 12 alcohol drinks in one year. 1 yes, 2 no, 7 refused,
# 9 don't know. Era 1 only.
des_sub <- update(des_sub, alq101 = dplyr::if_else(ALQ101 %in% c(7, 9), NA_real_, as.numeric(ALQ101)))
# ALQ110, ALQ_G/_H/_I: had at least 12 alcohol drinks in lifetime. 1 yes, 2 no, 7 refused,
# 9 don't know. Asked only of those answering no to ALQ101. Era 1 only.
des_sub <- update(des_sub, alq110 = dplyr::if_else(ALQ110 %in% c(7, 9), NA_real_, as.numeric(ALQ110)))
# ALQ111, P_ALQ: ever had a drink of any kind of alcohol. 1 yes, 2 no, 7 refused, 9 don't
# know. Era 2 only, and a stricter abstention question than ALQ101 and ALQ110 together.
des_sub <- update(des_sub, alq111 = dplyr::if_else(ALQ111 %in% c(7, 9), NA_real_, as.numeric(ALQ111)))
# ALQ120Q, ALQ_G/_H/_I: how often drink alcohol over the past 12 months. 0 to 365 with the
# unit in ALQ120U; 777 refused, 999 don't know. Era 1 only.
des_sub <- update(des_sub, alq120q = dplyr::if_else(ALQ120Q %in% c(777, 999), NA_real_, as.numeric(ALQ120Q)))
# ALQ121, P_ALQ: past 12 months how often drank alcoholic beverages, an 11-level frequency
# scale; 77 refused, 99 don't know. Era 2 only.
des_sub <- update(des_sub, alq121 = dplyr::if_else(ALQ121 %in% c(77, 99), NA_real_, as.numeric(ALQ121)))
# ALQ130, ALQ all cycles: average drinks per drinking day, past 12 months. 1 to 15 where 15
# means 15 or more; 777 refused, 999 don't know.
des_sub <- update(des_sub, alq130 = dplyr::if_else(ALQ130 %in% c(777, 999), NA_real_, as.numeric(ALQ130)))
# ALQ151, ALQ all cycles: ever had 4 or 5 or more drinks every day. 1 yes, 2 no, 7 refused,
# 9 don't know. The only past-pattern item present in every cycle; arm 3 leans on it.
des_sub <- update(des_sub, alq151 = dplyr::if_else(ALQ151 %in% c(7, 9), NA_real_, as.numeric(ALQ151)))Drinking days per year
Research question for this block: how often did each person drink in the past twelve months, on one scale, when the two eras asked the question differently?
Era 1 gives a count and a unit (ALQ120Q with
ALQ120U, coded 1 week, 2 month, 3 year). Era 2 gives an
11-level frequency category (ALQ121), which is converted at
the midpoint of each category. The two are put on the same scale here,
and the paper reports by era for exactly this reason.
# ALQ120U, ALQ_G/_H/_I: unit for ALQ120Q. 1 = week, 2 = month, 3 = year. Era 1 only. It
# carries no refused or don't-know code, so it is read as it stands and keeps its NHANES
# name; it is the one alcohol item in this Part with no lowercase twin.
des_sub <- update(des_sub, dpy = dplyr::case_when( # drinking days per year
era == 1 & alq120q == 0 ~ 0, # era 1, no drinking days in the past year
era == 1 & ALQ120U == 1 ~ alq120q * 52, # era 1, the count is per week
era == 1 & ALQ120U == 2 ~ alq120q * 12, # era 1, the count is per month
era == 1 & ALQ120U == 3 ~ alq120q * 1, # era 1, the count is per year
era == 2 & alq121 == 0 ~ 0, # era 2, never in the last year
era == 2 & alq121 == 1 ~ 365, # era 2, every day
era == 2 & alq121 == 2 ~ 350, # era 2, nearly every day
era == 2 & alq121 == 3 ~ 182, # era 2, 3 to 4 times a week, midpoint 3.5 x 52
era == 2 & alq121 == 4 ~ 104, # era 2, 2 times a week
era == 2 & alq121 == 5 ~ 52, # era 2, once a week
era == 2 & alq121 == 6 ~ 30, # era 2, 2 to 3 times a month, midpoint 2.5 x 12
era == 2 & alq121 == 7 ~ 12, # era 2, once a month
era == 2 & alq121 == 8 ~ 9, # era 2, 7 to 11 times in the last year, midpoint 9
era == 2 & alq121 == 9 ~ 4.5, # era 2, 3 to 6 times, midpoint 4.5
era == 2 & alq121 == 10 ~ 1.5 # era 2, 1 to 2 times, midpoint 1.5
))
summary(des_sub$variables$dpy) # 0 to 365 by construction; missing where the item was not answered## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0 1.5 12.0 58.6 52.0 365.0 2042Drinks per week
Research question for this block: how much alcohol does each current drinker report in an average week?
# Drinking days per year times drinks per drinking day, divided by 52.
des_sub <- update(des_sub, dpwk = dpy * alq130 / 52) # average drinks per week
summary(des_sub$variables$dpwk) # right-skewed; the categories below are cut on it## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.019 0.231 1.000 4.387 4.154 164.000 3858quantile(des_sub$variables$dpwk, c(0.5, 0.75, 0.9, 0.95, 0.99), na.rm = TRUE) # where the cut points fall in the sample## 50% 75% 90% 95% 99%
## 1.000 4.154 12.000 20.192 42.000The six-category drinking variable
Research question for this block: who is a lifetime abstainer, who is a former drinker, and how much do the current drinkers drink?
The cut points come from specs/00-cover.md, which
reproduces the shared definitions from the printed cover of the arm
specifications. Its wording, quoted: “Drinking categories among current
drinkers, for comparability with Srivatsa and colleagues (2026):
occasional (one drink a week or fewer), light (2 to 7 a week), moderate
(8 to 14), heavy (more than 14), from drinking days per year times
drinks per drinking day, divided by 52.” The cut points are the same for
women and men. Because dpwk is continuous rather than a
whole number of drinks, the boundaries are read as the intervals that
partition it without a gap: at or below 1, above 1 up to 7, above 7 up
to 14, and above 14.
des_sub <- update(des_sub, drinkcat = factor(dplyr::case_when( # the exposure all three arms are cut from
era == 1 & alq101 == 2 & alq110 == 2 ~ "Lifetime abstainer", # era 1: never 12 in a year and never 12 in a lifetime
era == 2 & alq111 == 2 ~ "Lifetime abstainer", # era 2: never a drink of any kind
era == 1 & (alq101 == 1 | alq110 == 1) & alq120q == 0 ~ "Former drinker", # era 1: ever drank, no drinking days in the past year
era == 2 & alq111 == 1 & alq121 == 0 ~ "Former drinker", # era 2: ever drank, never in the last year
dpy > 0 & dpwk <= 1 ~ "Occasional", # cover: "occasional (one drink a week or fewer)"
dpy > 0 & dpwk > 1 & dpwk <= 7 ~ "Light", # cover: "light (2 to 7 a week)"
dpy > 0 & dpwk > 7 & dpwk <= 14 ~ "Moderate", # cover: "moderate (8 to 14)"
dpy > 0 & dpwk > 14 ~ "Heavy" # cover: "heavy (more than 14)"
), levels = c("Lifetime abstainer", "Former drinker", "Occasional", "Light", "Moderate", "Heavy")))
table(des_sub$variables$drinkcat, des_sub$variables$era, useNA = "ifany") # unweighted counts by era, missing shown##
## 1 2
## Lifetime abstainer 970 317
## Former drinker 1068 741
## Occasional 2134 1265
## Light 1383 819
## Moderate 430 190
## Heavy 300 206
## <NA> 574 188The pooled non-drinker reference (arm 2) and the past-heavy split (arm 3)
Research question for this block: what do the other two arms’ exposures look like when they are cut from the same variable?
des_sub <- update(des_sub, nondrink = factor(dplyr::case_when( # arm 2: lifetime abstainers and former drinkers pooled on purpose
drinkcat %in% c("Lifetime abstainer", "Former drinker") ~ "Non-drinker", # the standard construction the literature uses
drinkcat == "Occasional" ~ "Occasional", # current drinkers keep their categories
drinkcat == "Light" ~ "Light", # current drinkers keep their categories
drinkcat == "Moderate" ~ "Moderate", # current drinkers keep their categories
drinkcat == "Heavy" ~ "Heavy" # current drinkers keep their categories
), levels = c("Non-drinker", "Occasional", "Light", "Moderate", "Heavy")))
# ALQ151 is the only past-pattern item in every cycle, so past-heavy is the only past
# drinking level the data can reconstruct for a former drinker.
des_sub <- update(des_sub, frmheavy = factor(dplyr::case_when( # arm 3: former drinkers split on past heavy drinking
drinkcat == "Former drinker" & alq151 == 1 ~ "Former, past heavy", # ever had 4 or 5 or more drinks every day
drinkcat == "Former drinker" & alq151 == 2 ~ "Former, not past heavy" # never drank that much every day
), levels = c("Former, not past heavy", "Former, past heavy")))
# Arm 3's second estimand: the six categories with the former-drinker row split in two, which
# is as far as Srivatsa's reallocation idea can be taken when the data carry one past-pattern
# item and never say why a person stopped.
des_sub <- update(des_sub, drinkre = factor(dplyr::case_when( # arm 3: seven categories
drinkcat == "Lifetime abstainer" ~ "Lifetime abstainer", # the reference, unchanged
frmheavy == "Former, not past heavy" ~ "Former, not past heavy", # former drinkers who never drank 4 or 5 every day
frmheavy == "Former, past heavy" ~ "Former, past heavy", # former drinkers who did
drinkcat == "Occasional" ~ "Occasional", # current drinkers keep their categories
drinkcat == "Light" ~ "Light", # current drinkers keep their categories
drinkcat == "Moderate" ~ "Moderate", # current drinkers keep their categories
drinkcat == "Heavy" ~ "Heavy" # current drinkers keep their categories
), levels = c("Lifetime abstainer", "Former, not past heavy", "Former, past heavy",
"Occasional", "Light", "Moderate", "Heavy")))
table(des_sub$variables$nondrink, des_sub$variables$era, useNA = "ifany") # arm 2's exposure by era##
## 1 2
## Non-drinker 2038 1058
## Occasional 2134 1265
## Light 1383 819
## Moderate 430 190
## Heavy 300 206
## <NA> 574 188table(des_sub$variables$frmheavy, des_sub$variables$era, useNA = "ifany") # arm 3's split, among former drinkers only##
## 1 2
## Former, not past heavy 806 582
## Former, past heavy 262 158
## <NA> 5791 2986table(des_sub$variables$drinkre, des_sub$variables$era, useNA = "ifany") # arm 3's seven-category exposure by era##
## 1 2
## Lifetime abstainer 970 317
## Former, not past heavy 806 582
## Former, past heavy 262 158
## Occasional 2134 1265
## Light 1383 819
## Moderate 430 190
## Heavy 300 206
## <NA> 574 189Sex
Research question for this block: how is sex carried into the tables and the models?
# RIAGENDR, DEMO all cycles: gender. 1 = male, 2 = female. No refused or don't know code.
des_sub <- update(des_sub, sex = factor(RIAGENDR, levels = c(1, 2), labels = c("Male", "Female")))
table(des_sub$variables$sex) # unweighted counts##
## Male Female
## 5169 5416Race and ethnicity
Research question for this block: which race and ethnicity variable does this paper use, and why that one?
# RIDRETH3, DEMO all cycles: race/Hispanic origin with Non-Hispanic Asian. 1 Mexican
# American, 2 Other Hispanic, 3 Non-Hispanic White, 4 Non-Hispanic Black, 6 Non-Hispanic
# Asian, 7 Other or multiracial. RIDRETH3 rather than RIDRETH1 because code 6 is what
# triggers the Asian-specific body mass index and waist cut-points the outcome is built on.
des_sub <- update(des_sub, racecat = factor(RIDRETH3, levels = c(3, 1, 2, 4, 6, 7),
labels = c("Non-Hispanic White", "Mexican American", "Other Hispanic",
"Non-Hispanic Black", "Non-Hispanic Asian", "Other or multiracial")))
table(des_sub$variables$racecat, useNA = "ifany") # Non-Hispanic White is the reference level##
## Non-Hispanic White Mexican American Other Hispanic Non-Hispanic Black
## 3964 1409 1157 2335
## Non-Hispanic Asian Other or multiracial
## 1334 386Education
Research question for this block: what is each person’s completed education?
# DMDEDUC2, DEMO all cycles: education level, adults 20+. 1 less than 9th grade, 2 9th to
# 11th grade, 3 high school graduate or GED, 4 some college or associate degree, 5 college
# graduate or above, 7 refused, 9 don't know.
des_sub <- update(des_sub, educat = factor(dplyr::if_else(DMDEDUC2 %in% c(7, 9), NA_real_, as.numeric(DMDEDUC2)),
levels = 1:5,
labels = c("Less than 9th grade", "9th to 11th grade", "High school or GED",
"Some college", "College graduate")))
table(des_sub$variables$educat, useNA = "ifany") # refused and don't know are now missing##
## Less than 9th grade 9th to 11th grade High school or GED Some college College graduate
## 967 1328 2369 3219 2697
## <NA>
## 5Income to poverty ratio
Research question for this block: what is each family’s income relative to the poverty threshold?
# INDFMPIR, DEMO all cycles: ratio of family income to poverty, 0 to 5, top-coded at 5.
# No refused or don't know code; the missing values are unreported income.
des_sub <- update(des_sub, pir = INDFMPIR) # carried as reported, on its own scale
summary(des_sub$variables$pir) # the missing count here is what the complete-case models lose most to## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00 1.10 2.12 2.52 4.07 5.00 1059Smoking
Research question for this block: never, former, or current smoker?
# SMQ020, SMQ all cycles: smoked at least 100 cigarettes in life. 1 yes, 2 no, 7 refused,
# 9 don't know. SMQ040, SMQ all cycles: do you now smoke cigarettes. 1 every day, 2 some
# days, 3 not at all; cleaned by the build, and asked only of those answering yes to SMQ020,
# so a missing SMQ040 with SMQ020 = 2 means never-smoker rather than an unanswered question.
des_sub <- update(des_sub, smokecat = factor(dplyr::case_when(
SMQ020 == 2 ~ "Never", # fewer than 100 cigarettes in life
SMQ020 == 1 & SMQ040 == 3 ~ "Former", # smoked 100, does not smoke now
SMQ020 == 1 & SMQ040 %in% c(1, 2) ~ "Current" # smoked 100, smokes every day or some days
), levels = c("Never", "Former", "Current")))
table(des_sub$variables$smokecat, useNA = "ifany") # the missing are the refused and don't know answers to SMQ020##
## Never Former Current <NA>
## 6003 2572 1999 11Self-rated health
Research question for this block: how does each person rate their own health?
The arm specifications name HSD010. HSD010
is missing for every era 2 row, because the 2017 to March 2020 Current
Health Status file does not carry it. HUQ010 is the same
question in the Hospital Utilization file, present in all four cycles,
and the codebook says it is the self-rated health variable the arms
should use. This document uses HUQ010.
# HUQ010, HUQ all cycles: general health condition. 1 excellent, 2 very good, 3 good,
# 4 fair, 5 poor, 7 refused, 9 don't know.
des_sub <- update(des_sub, srh = factor(dplyr::if_else(HUQ010 %in% c(7, 9), NA_real_, as.numeric(HUQ010)),
levels = 1:5, labels = c("Excellent", "Very good", "Good", "Fair", "Poor")))
table(des_sub$variables$srh, useNA = "ifany") # refused and don't know are now missing##
## Excellent Very good Good Fair Poor <NA>
## 1199 2743 4010 2169 455 9Physical activity
Research question for this block: does the person report any vigorous or moderate activity, at work or in their leisure time?
The cover’s definition names four Global Physical Activity
Questionnaire items: PAQ605 vigorous work,
PAQ620 moderate work, PAQ650 vigorous
recreational, and PAQ665 moderate recreational.
PAQ665 is not in this analytic file. The
build’s keep-list carries PAQ605, PAQ610,
PAD615, PAQ620, PAQ625,
PAD630, PAQ650, PAQ655,
PAD660, and PAD680, and the codebook lists the
same ten. physact is therefore built from the three yes/no
items that are present, which is the closest the file allows. It leaves
8 of the 10,585 analytic rows missing rather than the handful the
four-item version would. Anyone who wants the fourth item has to rebuild
the file, and Paper 2 has to use these same three so that the two
papers’ adjustment sets stay identical.
# PAQ605, PAQ all cycles: vigorous work activity. 1 yes, 2 no, 7 refused, 9 don't know.
des_sub <- update(des_sub, paq605 = dplyr::if_else(PAQ605 %in% c(7, 9), NA_real_, as.numeric(PAQ605)))
# PAQ620, PAQ all cycles: moderate work activity. 1 yes, 2 no, 7 refused, 9 don't know.
des_sub <- update(des_sub, paq620 = dplyr::if_else(PAQ620 %in% c(7, 9), NA_real_, as.numeric(PAQ620)))
# PAQ650, PAQ all cycles: vigorous recreational activities. 1 yes, 2 no, 7 refused, 9 don't know.
des_sub <- update(des_sub, paq650 = dplyr::if_else(PAQ650 %in% c(7, 9), NA_real_, as.numeric(PAQ650)))
# Any activity is a yes on any of the three; no activity is a no on all three; anything else
# cannot be resolved and stays missing.
des_sub <- update(des_sub, physact = factor(dplyr::case_when(
paq605 == 1 | paq620 == 1 | paq650 == 1 ~ 1, # any vigorous or moderate activity reported
paq605 == 2 & paq620 == 2 & paq650 == 2 ~ 0 # no to all three items
), levels = c(0, 1), labels = c("No", "Yes")))
table(des_sub$variables$physact, des_sub$variables$era, useNA = "ifany") # unweighted counts by era, missing shown##
## 1 2
## No 3363 1577
## Yes 3493 2144
## <NA> 3 5The outcome as a factor, and the cycle label
Research question for this block: what does a CKM stage look like in a table?
# ckmstage, derived: CKM syndrome stage 0 to 4, hierarchical, the first criterion met
# reading down from 4. See the codebook's derived-outcome row for every threshold.
des_sub <- update(des_sub, stagef = factor(ckmstage, levels = 0:4,
labels = c("Stage 0", "Stage 1", "Stage 2", "Stage 3", "Stage 4")))
# ckmtier, derived: 1 = stages 0 or 1, 2 = stage 2, 3 = stages 3 or 4. The fallback outcome
# that Part 5 switches to when a cell is too thin.
des_sub <- update(des_sub, tierf = factor(ckmtier, levels = 1:3,
labels = c("Low (0 to 1)", "Moderate (2)", "High (3 to 4)")))
# yrstr, derived: the cycle as printed in tables.
des_sub <- update(des_sub, eraf = factor(era, levels = c(1, 2),
labels = c("2011-2016 (ALQ101, ALQ110)", "2017-2020 (ALQ111)")))
table(des_sub$variables$stagef, des_sub$variables$yrstr) # unweighted stage by cycle, the build log prints the same table##
## 2011-2012 2013-2014 2015-2016 2017-2020
## Stage 0 197 210 146 224
## Stage 1 509 552 519 825
## Stage 2 1144 1176 1114 1913
## Stage 3 182 172 189 287
## Stage 4 239 244 266 477Complete cases for the adjusted models
Research question for this block: how many rows can the adjusted models actually use?
The models below are complete-case on the exposure and the shared covariate set. That is the simplest defensible choice and it is what the SAS and Stata twins do by default, so the three languages agree. The count is printed here so that the loss is visible rather than silent, and the income to poverty ratio is most of it.
des_sub <- update(des_sub, modflag = as.integer( # 1 if every model variable is observed
!is.na(drinkcat) & !is.na(RIDAGEYR) & !is.na(sex) & !is.na(racecat) &
!is.na(educat) & !is.na(pir) & !is.na(smokecat) & !is.na(srh) & !is.na(physact)))
table(des_sub$variables$modflag, des_sub$variables$era) # complete cases by era##
## 1 2
## 0 1094 649
## 1 5765 3077## [1] 8842## [1] 1743Part 5: Cell check
Research question for this block: is any exposure-by-stage cell too thin for a five-level ordinal outcome?
The rule is written in code/README.md Part 5 and
promised in all three arm specifications. If any cell of the unweighted
exposure-by-stage table holds fewer than 30 people, or if any stage is
empty in the analytic rows, the outcome becomes ckmtier (0
to 1 low, 2 moderate, 3 to 4 high) instead of ckmstage. The
check runs on the six-category drinkcat, which is the
finest exposure any arm uses, so one verdict covers all three arms.
tab <- table(des_sub$variables$drinkcat, des_sub$variables$ckmstage) # unweighted exposure by stage, analytic rows
tab # printed so the thinnest cell can be read off##
## 0 1 2 3 4
## Lifetime abstainer 89 258 653 138 149
## Former drinker 55 245 867 255 387
## Occasional 277 838 1697 239 348
## Light 228 625 1103 100 146
## Moderate 40 137 348 32 63
## Heavy 31 84 315 23 53## [1] 23##
## 0 1 2 3 4
## 777 2405 5347 830 1226# The rule, written out.
if (min(tab) < 30 | any(table(des_sub$variables$ckmstage) == 0)) outcome <- des_sub$variables$ckmtier else outcome <- des_sub$variables$ckmstage
if (min(tab) < 30 | any(table(des_sub$variables$ckmstage) == 0)) outcname <- "ckmtier (1 low, 2 moderate, 3 high)" else outcname <- "ckmstage (0 to 4)"
des_sub <- update(des_sub, outcome = outcome) # the design now carries the outcome the rule chose
cat("Cell check: the thinnest exposure-by-stage cell holds", min(tab),
"people, so the outcome used in every model below is", outcname, "\n")## Cell check: the thinnest exposure-by-stage cell holds 23 people, so the outcome used in every model below is ckmtier (1 low, 2 moderate, 3 high)##
## 1 2 3
## 3182 5347 2056The cell check runs on all 10,585 analytic rows. The era-specific models below fit on about two thirds and one third of those rows, so a cell that clears 30 in the pooled table can be much thinner inside one era. The era-specific tables are printed next so that a reader can see how thin, and the era estimates are read with that in mind rather than with a second rule.
table(des_sub$variables$drinkcat[des_sub$variables$era == 1], # era 1 exposure by stage
des_sub$variables$ckmstage[des_sub$variables$era == 1])##
## 0 1 2 3 4
## Lifetime abstainer 76 189 476 109 120
## Former drinker 31 142 511 162 222
## Occasional 176 546 1051 147 214
## Light 165 392 680 65 81
## Moderate 32 100 232 24 42
## Heavy 20 46 198 10 26table(des_sub$variables$drinkcat[des_sub$variables$era == 2], # era 2 exposure by stage
des_sub$variables$ckmstage[des_sub$variables$era == 2])##
## 0 1 2 3 4
## Lifetime abstainer 13 69 177 29 29
## Former drinker 24 103 356 93 165
## Occasional 101 292 646 92 134
## Light 63 233 423 35 65
## Moderate 8 37 116 8 21
## Heavy 11 38 117 13 27Table 1: the analytic sample, weighted
Research question for this block: who is in the analytic sample, and how do the six drinking groups differ before any model is fitted?
des_t1 <- subset(des_sub, !is.na(drinkcat)) # Table 1 describes the rows with a drinking category
nrow(des_t1) # analytic rows with the exposure observed## [1] 9823##
## Lifetime abstainer Former drinker Occasional Light Moderate
## 1287 1809 3399 2202 620
## Heavy
## 506tbl1 <- tbl_svysummary( # weighted descriptive table
des_t1, # the design, so every percentage is survey weighted
by = drinkcat, # one column per drinking category
include = c(RIDAGEYR, sex, racecat, educat, pir, smokecat, srh, physact, eraf, stagef), # the shared covariate set plus era and stage
label = list(RIDAGEYR ~ "Age at screening, years", sex ~ "Sex",
racecat ~ "Race and ethnicity", educat ~ "Education",
pir ~ "Family income to poverty ratio", smokecat ~ "Smoking status",
srh ~ "Self-rated health (HUQ010)", physact ~ "Any vigorous or moderate physical activity",
eraf ~ "Alcohol instrument era",
stagef ~ "CKM syndrome stage"),
statistic = list(all_continuous() ~ "{mean} ({sd})"), # mean and standard deviation for the two continuous variables
digits = list(all_continuous() ~ 1), # one decimal place
missing_text = "Missing" # missing shown rather than dropped
) |>
add_overall() |> # an all-categories column
modify_caption("Table 1. The analytic sample by drinking category, NHANES 2011 to March 2020, survey weighted")
tbl1 # print it in the knitted page| Characteristic | Overall N = 215,873,1221 |
Lifetime abstainer N = 21,264,7781 |
Former drinker N = 32,912,2561 |
Occasional N = 74,795,7201 |
Light N = 57,615,1121 |
Moderate N = 17,022,6011 |
Heavy N = 12,262,6561 |
|---|---|---|---|---|---|---|---|
| Age at screening, years | 48.2 (17.1) | 50.4 (18.7) | 56.9 (15.7) | 47.5 (16.8) | 43.9 (16.3) | 46.5 (16.1) | 47.6 (15.2) |
| Sex | |||||||
| Male | 106,500,665 (49%) | 6,672,667 (31%) | 15,128,404 (46%) | 30,891,662 (41%) | 32,720,547 (57%) | 11,098,031 (65%) | 9,989,354 (81%) |
| Female | 109,372,457 (51%) | 14,592,111 (69%) | 17,783,851 (54%) | 43,904,058 (59%) | 24,894,564 (43%) | 5,924,570 (35%) | 2,273,302 (19%) |
| Race and ethnicity | |||||||
| Non-Hispanic White | 141,344,013 (65%) | 10,693,168 (50%) | 21,330,279 (65%) | 47,970,956 (64%) | 39,886,769 (69%) | 12,594,519 (74%) | 8,868,322 (72%) |
| Mexican American | 18,510,973 (8.6%) | 2,413,824 (11%) | 2,805,093 (8.5%) | 6,536,232 (8.7%) | 4,639,686 (8.1%) | 1,188,868 (7.0%) | 927,271 (7.6%) |
| Other Hispanic | 13,803,826 (6.4%) | 1,749,650 (8.2%) | 2,283,331 (6.9%) | 4,964,642 (6.6%) | 3,732,733 (6.5%) | 642,236 (3.8%) | 431,234 (3.5%) |
| Non-Hispanic Black | 24,053,349 (11%) | 2,989,683 (14%) | 3,763,013 (11%) | 8,611,913 (12%) | 5,561,256 (9.7%) | 1,705,035 (10%) | 1,422,449 (12%) |
| Non-Hispanic Asian | 10,667,181 (4.9%) | 2,911,543 (14%) | 1,259,279 (3.8%) | 4,012,077 (5.4%) | 1,889,499 (3.3%) | 310,609 (1.8%) | 284,175 (2.3%) |
| Other or multiracial | 7,493,779 (3.5%) | 506,909 (2.4%) | 1,471,261 (4.5%) | 2,699,901 (3.6%) | 1,905,170 (3.3%) | 581,333 (3.4%) | 329,205 (2.7%) |
| Education | |||||||
| Less than 9th grade | 9,780,494 (4.5%) | 2,426,842 (11%) | 2,700,087 (8.2%) | 2,447,473 (3.3%) | 1,613,198 (2.8%) | 207,616 (1.2%) | 385,277 (3.1%) |
| 9th to 11th grade | 20,149,158 (9.3%) | 2,592,208 (12%) | 4,780,599 (15%) | 6,145,911 (8.2%) | 3,940,649 (6.8%) | 1,335,766 (7.8%) | 1,354,025 (11%) |
| High school or GED | 50,226,040 (23%) | 5,422,331 (26%) | 9,438,354 (29%) | 17,841,630 (24%) | 10,990,912 (19%) | 3,393,720 (20%) | 3,139,093 (26%) |
| Some college | 67,762,881 (31%) | 5,889,535 (28%) | 9,279,385 (28%) | 25,482,045 (34%) | 18,044,581 (31%) | 4,779,621 (28%) | 4,287,714 (35%) |
| College graduate | 67,928,571 (31%) | 4,928,856 (23%) | 6,697,226 (20%) | 22,878,662 (31%) | 23,025,773 (40%) | 7,305,877 (43%) | 3,092,177 (25%) |
| Missing | 25,979 | 5,006 | 16,605 | 0 | 0 | 0 | 4,368 |
| Family income to poverty ratio | 3.0 (1.6) | 2.4 (1.5) | 2.5 (1.6) | 3.0 (1.6) | 3.4 (1.6) | 3.4 (1.6) | 3.0 (1.7) |
| Missing | 16,872,485 | 2,381,469 | 3,241,709 | 5,372,405 | 4,230,194 | 1,021,476 | 625,232 |
| Smoking status | |||||||
| Never | 119,808,299 (56%) | 19,080,571 (90%) | 14,391,978 (44%) | 45,709,014 (61%) | 31,506,689 (55%) | 6,148,574 (36%) | 2,971,473 (24%) |
| Former | 56,741,253 (26%) | 1,289,544 (6.1%) | 12,005,845 (37%) | 17,833,546 (24%) | 15,657,921 (27%) | 6,041,649 (36%) | 3,912,748 (32%) |
| Current | 39,225,229 (18%) | 894,663 (4.2%) | 6,489,108 (20%) | 11,207,656 (15%) | 10,443,054 (18%) | 4,825,830 (28%) | 5,364,918 (44%) |
| Missing | 98,341 | 0 | 25,324 | 45,504 | 7,449 | 6,548 | 13,516 |
| Self-rated health (HUQ010) | |||||||
| Excellent | 28,987,334 (13%) | 2,762,335 (13%) | 2,697,298 (8.2%) | 9,567,472 (13%) | 10,393,106 (18%) | 2,202,529 (13%) | 1,364,594 (11%) |
| Very good | 67,832,117 (31%) | 6,029,091 (28%) | 7,062,538 (21%) | 24,203,290 (32%) | 20,141,441 (35%) | 6,757,202 (40%) | 3,638,556 (30%) |
| Good | 78,912,443 (37%) | 7,437,449 (35%) | 12,904,695 (39%) | 28,074,899 (38%) | 19,829,498 (34%) | 5,956,810 (35%) | 4,709,092 (38%) |
| Fair | 33,617,602 (16%) | 4,144,359 (19%) | 7,892,626 (24%) | 11,343,928 (15%) | 6,352,326 (11%) | 1,831,950 (11%) | 2,052,413 (17%) |
| Poor | 6,334,866 (2.9%) | 881,947 (4.1%) | 2,316,230 (7.0%) | 1,541,961 (2.1%) | 826,520 (1.4%) | 274,111 (1.6%) | 494,098 (4.0%) |
| Missing | 188,760 | 9,596 | 38,869 | 64,171 | 72,221 | 0 | 3,903 |
| Any vigorous or moderate physical activity | 129,096,918 (60%) | 10,068,093 (47%) | 16,575,341 (50%) | 44,289,585 (59%) | 38,236,558 (66%) | 11,366,311 (67%) | 8,561,030 (70%) |
| Missing | 30,175 | 0 | 22,612 | 5,558 | 0 | 0 | 2,004 |
| Alcohol instrument era | |||||||
| 2011-2016 (ALQ101, ALQ110) | 136,683,843 (63%) | 16,294,289 (77%) | 19,586,807 (60%) | 46,426,585 (62%) | 35,671,528 (62%) | 11,259,204 (66%) | 7,445,430 (61%) |
| 2017-2020 (ALQ111) | 79,189,279 (37%) | 4,970,489 (23%) | 13,325,449 (40%) | 28,369,135 (38%) | 21,943,584 (38%) | 5,763,397 (34%) | 4,817,225 (39%) |
| CKM syndrome stage | |||||||
| Stage 0 | 20,343,706 (9.4%) | 1,973,683 (9.3%) | 1,368,863 (4.2%) | 6,993,413 (9.4%) | 7,612,149 (13%) | 1,293,487 (7.6%) | 1,102,112 (9.0%) |
| Stage 1 | 54,718,736 (25%) | 4,696,456 (22%) | 5,220,843 (16%) | 19,976,014 (27%) | 18,569,605 (32%) | 4,077,996 (24%) | 2,177,823 (18%) |
| Stage 2 | 108,421,752 (50%) | 10,418,667 (49%) | 16,359,652 (50%) | 37,610,339 (50%) | 26,765,330 (46%) | 9,664,347 (57%) | 7,603,417 (62%) |
| Stage 3 | 11,408,281 (5.3%) | 1,870,055 (8.8%) | 3,332,233 (10%) | 3,821,162 (5.1%) | 1,521,839 (2.6%) | 454,746 (2.7%) | 408,246 (3.3%) |
| Stage 4 | 20,980,647 (9.7%) | 2,305,917 (11%) | 6,630,666 (20%) | 6,394,792 (8.5%) | 3,146,189 (5.5%) | 1,532,026 (9.0%) | 971,058 (7.9%) |
| 1 Mean (SD); n (%) | |||||||
tbl1_df <- as.data.frame(tbl1) # the same table as a plain data frame
write.csv(tbl1_df, file.path(out_dir, "table-1.csv"), row.names = FALSE) # csv first, the file that is compared
flextable::save_as_docx(gtsummary::as_flex_table(tbl1), # then Word, which carries a write timestamp
path = file.path(out_dir, "table-1.docx"))
list.files(out_dir) # what landed in the output directory## [1] "arm-01-lifetime-abstainers.html" "figure-1.png"
## [3] "run-log.txt" "table-1.csv"
## [5] "table-1.docx" "table-2.csv"
## [7] "table-2.docx" "table-3-ckmstage.csv"
## [9] "table-3-ckmstage.docx" "table-3.csv"
## [11] "table-3.docx"Part 6: Model
What this arm estimates. For each drinking level, the proportional odds of a higher CKM stage compared with people who have never drunk alcohol. Former drinkers are kept in the model as their own row and are never part of the reference. The reference group is what carries the identification here, because the target trial, which would assign a drinking pattern for the decade before examination, cannot be run.
What the arm removes and what it adds. Lifetime abstainers never quit, so the sick-quitter arrow from prior illness into the exposure is gone. In its place is a new backdoor path through whatever makes a person abstain for life, such as religion, culture, income, and illness early enough to have prevented drinking at all. The covariate set closes part of that path and cannot close the unmeasured part.
The outcome the pre-specified rule chose. Part 5
collapsed the five stages to the three risk tiers, because the thinnest
exposure-by-stage cell holds fewer than 30 people. The models below
therefore run on ckmtier. The five-stage version of the
same models is reported as a sensitivity analysis and written to
table-3-ckmstage.csv, shown only because the rule collapsed
the primary outcome. The arm specification asks for a three-tier version
and a binary high-risk version as checks; the rule made the three-tier
version primary, so the five-stage version takes the place of the first
check and the binary high-risk check is unchanged.
What this arm does not yet produce. The specification’s estimand names adjusted stage prevalences alongside the odds ratios. Table 2 reports crude survey-weighted stage distributions, not model-standardised prevalences, and the arm’s Methods should say so until the team decides how to standardise them.
The rows this arm models
Research question for this block: how many people have a drinking category, and how many of those have every covariate?
des_exp <- subset(des_sub, !is.na(drinkcat)) # analytic rows with a drinking category
n_exp <- nrow(des_exp) # rows the descriptive tables use
n_exp # printed## [1] 9823## [1] 762des_mod <- subset(des_exp, modflag == 1) # complete cases on the shared covariate set
n_mod <- nrow(des_mod) # rows the adjusted models use
n_mod # printed## [1] 8842des_e1 <- subset(des_mod, era == 1) # 2011 to 2016, ALQ101 and ALQ110
nrow(des_e1) # era 1 model rows## [1] 5765## [1] 3077##
## 1 2
## Lifetime abstainer 848 260
## Former drinker 969 622
## Occasional 1972 1110
## Light 1292 727
## Moderate 406 169
## Heavy 278 189table(des_mod$variables$drinkcat, des_mod$variables$outcome) # model rows by category and the outcome the rule chose##
## 1 2 3
## Lifetime abstainer 302 558 248
## Former drinker 268 762 561
## Occasional 1019 1539 524
## Light 780 1021 218
## Moderate 171 316 88
## Heavy 107 296 64The weighted stage distribution, pooled
Research question for this block: what share of each drinking group is at each CKM stage, in the population the weights represent?
sb_p <- svyby(~stagef, ~drinkcat, des_exp, svymean, na.rm = TRUE) # weighted stage distribution by category
sb_p # proportions and their standard errors## drinkcat stagefStage 0 stagefStage 1 stagefStage 2 stagefStage 3
## Lifetime abstainer Lifetime abstainer 0.09281 0.2209 0.4899 0.08794
## Former drinker Former drinker 0.04159 0.1586 0.4971 0.10125
## Occasional Occasional 0.09350 0.2671 0.5028 0.05109
## Light Light 0.13212 0.3223 0.4646 0.02641
## Moderate Moderate 0.07599 0.2396 0.5677 0.02671
## Heavy Heavy 0.08988 0.1776 0.6200 0.03329
## stagefStage 4 se.stagefStage 0 se.stagefStage 1 se.stagefStage 2
## Lifetime abstainer 0.10844 0.012255 0.01317 0.01880
## Former drinker 0.20146 0.007539 0.01184 0.01637
## Occasional 0.08550 0.009882 0.01210 0.01264
## Light 0.05461 0.010140 0.01538 0.01769
## Moderate 0.09000 0.015695 0.02414 0.02703
## Heavy 0.07919 0.019391 0.02415 0.02903
## se.stagefStage 3 se.stagefStage 4
## Lifetime abstainer 0.010453 0.012735
## Former drinker 0.007936 0.013119
## Occasional 0.004390 0.006181
## Light 0.004115 0.007205
## Moderate 0.006066 0.021664
## Heavy 0.009167 0.014851ci_p <- confint(sb_p) # 95% confidence limits, ordered stage by stage within category
round(100 * as.matrix(sb_p[, 2:6]), 2) # the same proportions as percentages, for reading## stagefStage 0 stagefStage 1 stagefStage 2 stagefStage 3 stagefStage 4
## Lifetime abstainer 9.28 22.09 48.99 8.79 10.84
## Former drinker 4.16 15.86 49.71 10.12 20.15
## Occasional 9.35 26.71 50.28 5.11 8.55
## Light 13.21 32.23 46.46 2.64 5.46
## Moderate 7.60 23.96 56.77 2.67 9.00
## Heavy 8.99 17.76 62.00 3.33 7.92The weighted stage distribution, era 1
Research question for this block: does the picture look the same when lifetime abstention means never 12 drinks in a year and never 12 in a lifetime?
des_x1 <- subset(des_exp, era == 1) # era 1 rows with a drinking category
nrow(des_x1) # era 1 descriptive rows## [1] 6285sb_1 <- svyby(~stagef, ~drinkcat, des_x1, svymean, na.rm = TRUE) # weighted stage distribution, era 1
sb_1 # proportions and standard errors## drinkcat stagefStage 0 stagefStage 1 stagefStage 2 stagefStage 3
## Lifetime abstainer Lifetime abstainer 0.09613 0.2303 0.4724 0.09220
## Former drinker Former drinker 0.04185 0.1570 0.4973 0.10076
## Occasional Occasional 0.08601 0.2715 0.5064 0.05016
## Light Light 0.14827 0.3064 0.4714 0.02769
## Moderate Moderate 0.08545 0.2635 0.5309 0.02838
## Heavy Heavy 0.08043 0.1464 0.6681 0.02886
## stagefStage 4 se.stagefStage 0 se.stagefStage 1 se.stagefStage 2
## Lifetime abstainer 0.10891 0.01426 0.01522 0.02096
## Former drinker 0.20309 0.01025 0.01535 0.02110
## Occasional 0.08602 0.01004 0.01332 0.01334
## Light 0.04624 0.01178 0.01722 0.01895
## Moderate 0.09178 0.01880 0.03137 0.02754
## Heavy 0.07624 0.02001 0.02624 0.03117
## se.stagefStage 3 se.stagefStage 4
## Lifetime abstainer 0.012642 0.013198
## Former drinker 0.008987 0.016283
## Occasional 0.005567 0.007183
## Light 0.005784 0.006936
## Moderate 0.006308 0.016751
## Heavy 0.012940 0.019284The weighted stage distribution, era 2
Research question for this block: and when lifetime abstention means never a drink of any kind?
des_x2 <- subset(des_exp, era == 2) # era 2 rows with a drinking category
nrow(des_x2) # era 2 descriptive rows## [1] 3538sb_2 <- svyby(~stagef, ~drinkcat, des_x2, svymean, na.rm = TRUE) # weighted stage distribution, era 2
sb_2 # proportions and standard errors## drinkcat stagefStage 0 stagefStage 1 stagefStage 2 stagefStage 3
## Lifetime abstainer Lifetime abstainer 0.08194 0.1898 0.5474 0.07397
## Former drinker Former drinker 0.04122 0.1610 0.4968 0.10195
## Occasional Occasional 0.10576 0.2599 0.4971 0.05260
## Light Light 0.10586 0.3481 0.4535 0.02434
## Moderate Moderate 0.05750 0.1928 0.6397 0.02346
## Heavy Heavy 0.10447 0.2259 0.5458 0.04014
## stagefStage 4 se.stagefStage 0 se.stagefStage 1 se.stagefStage 2
## Lifetime abstainer 0.10688 0.02324 0.02509 0.03987
## Former drinker 0.19908 0.01096 0.01855 0.02594
## Occasional 0.08463 0.02008 0.02338 0.02508
## Light 0.06821 0.01895 0.02818 0.03441
## Moderate 0.08653 0.02807 0.03119 0.05662
## Heavy 0.08374 0.03851 0.04591 0.05670
## se.stagefStage 3 se.stagefStage 4
## Lifetime abstainer 0.016813 0.03314
## Former drinker 0.014469 0.02182
## Occasional 0.007127 0.01128
## Light 0.005304 0.01538
## Moderate 0.012961 0.05514
## Heavy 0.011998 0.02339The proportional-odds model, era 1
Research question for this block: in 2011 to 2016, what are the odds of a higher CKM tier for each drinking level compared with lifetime abstainers?
fit_e1 <- svyolr(factor(outcome, ordered = TRUE) ~ drinkcat + RIDAGEYR + sex + racecat +
educat + pir + smokecat + srh + physact, design = des_e1) # era 1, the outcome Part 5 chose
or_e1 <- as.data.frame(broom::tidy(fit_e1, conf.int = TRUE, exponentiate = TRUE))[, c("term", "estimate", "conf.low", "conf.high")] # exponentiated coefficients with 95% confidence limits
or_e1$model <- rep("Era 1: 2011-2016", nrow(or_e1)) # label the block for Table 3
or_e1$nrows <- rep(nrow(des_e1), nrow(or_e1)) # the rows this model used
or_e1 # every term, so the covariate estimates are on the record too## term estimate conf.low conf.high model nrows
## 1 drinkcatFormer drinker 1.2573 1.0029 1.5763 Era 1: 2011-2016 5765
## 2 drinkcatOccasional 0.8947 0.7421 1.0787 Era 1: 2011-2016 5765
## 3 drinkcatLight 0.7123 0.5701 0.8899 Era 1: 2011-2016 5765
## 4 drinkcatModerate 0.9492 0.6696 1.3456 Era 1: 2011-2016 5765
## 5 drinkcatHeavy 1.0153 0.7508 1.3731 Era 1: 2011-2016 5765
## 6 RIDAGEYR 1.0851 1.0777 1.0925 Era 1: 2011-2016 5765
## 7 sexFemale 0.5427 0.4685 0.6286 Era 1: 2011-2016 5765
## 8 racecatMexican American 0.8929 0.7253 1.0992 Era 1: 2011-2016 5765
## 9 racecatOther Hispanic 0.9316 0.7241 1.1985 Era 1: 2011-2016 5765
## 10 racecatNon-Hispanic Black 1.2766 1.0804 1.5085 Era 1: 2011-2016 5765
## 11 racecatNon-Hispanic Asian 0.9016 0.7156 1.1359 Era 1: 2011-2016 5765
## 12 racecatOther or multiracial 1.4187 0.9106 2.2103 Era 1: 2011-2016 5765
## 13 educat9th to 11th grade 1.4274 0.9943 2.0493 Era 1: 2011-2016 5765
## 14 educatHigh school or GED 1.2348 0.9065 1.6820 Era 1: 2011-2016 5765
## 15 educatSome college 1.5067 1.1404 1.9906 Era 1: 2011-2016 5765
## 16 educatCollege graduate 1.0714 0.7640 1.5027 Era 1: 2011-2016 5765
## 17 pir 0.8776 0.8248 0.9338 Era 1: 2011-2016 5765
## 18 smokecatFormer 1.0222 0.8418 1.2412 Era 1: 2011-2016 5765
## 19 smokecatCurrent 1.0881 0.8908 1.3291 Era 1: 2011-2016 5765
## 20 srhVery good 1.9596 1.4335 2.6788 Era 1: 2011-2016 5765
## 21 srhGood 3.0381 2.2374 4.1254 Era 1: 2011-2016 5765
## 22 srhFair 4.9598 3.4795 7.0701 Era 1: 2011-2016 5765
## 23 srhPoor 7.5058 4.4367 12.6980 Era 1: 2011-2016 5765
## 24 physactYes 0.8379 0.6979 1.0060 Era 1: 2011-2016 5765
## 25 1|2 24.9695 15.2015 41.0142 Era 1: 2011-2016 5765
## 26 2|3 833.5012 478.8815 1450.7226 Era 1: 2011-2016 5765The proportional-odds model, era 2
Research question for this block: and in 2017 to March 2020, where the abstention question is stricter?
The era 2 fits are the only ones in this document wrapped in
try(). One era on its own is the smallest sample any model
here uses, and the 300-row teaching sample under
data/sample/ cannot fit it. The guard prints the condition
message, so a failure on the real file says why rather than leaving an
empty row in Table 3. Every other fit runs bare and will stop the knit
if it fails.
fit_e2 <- try(svyolr(factor(outcome, ordered = TRUE) ~ drinkcat + RIDAGEYR + sex + racecat +
educat + pir + smokecat + srh + physact, design = des_e2), silent = TRUE) # guarded, see the box above
if (inherits(fit_e2, "try-error")) cat("Era 2 model did not fit on this file:", conditionMessage(attr(fit_e2, "condition")), "\n") # the reason travels with the failure
or_e2 <- data.frame(term = character(0), estimate = numeric(0), conf.low = numeric(0), conf.high = numeric(0)) # stays empty if the fit failed
if (!inherits(fit_e2, "try-error")) or_e2 <- as.data.frame(broom::tidy(fit_e2, conf.int = TRUE, exponentiate = TRUE))[, c("term", "estimate", "conf.low", "conf.high")] # exponentiated coefficients with 95% limits, kept only if the fit succeeded
or_e2$model <- rep("Era 2: 2017-2020", nrow(or_e2)) # label the block for Table 3
or_e2$nrows <- rep(nrow(des_e2), nrow(or_e2)) # the rows this model used
or_e2 # every term## term estimate conf.low conf.high model nrows
## 1 drinkcatFormer drinker 0.8512 0.5033 1.4394 Era 2: 2017-2020 3077
## 2 drinkcatOccasional 0.7861 0.4923 1.2551 Era 2: 2017-2020 3077
## 3 drinkcatLight 0.6650 0.4203 1.0523 Era 2: 2017-2020 3077
## 4 drinkcatModerate 1.1350 0.6482 1.9874 Era 2: 2017-2020 3077
## 5 drinkcatHeavy 0.5670 0.2942 1.0929 Era 2: 2017-2020 3077
## 6 RIDAGEYR 1.0711 1.0656 1.0766 Era 2: 2017-2020 3077
## 7 sexFemale 0.6525 0.4889 0.8708 Era 2: 2017-2020 3077
## 8 racecatMexican American 0.8135 0.5701 1.1608 Era 2: 2017-2020 3077
## 9 racecatOther Hispanic 0.7190 0.4732 1.0925 Era 2: 2017-2020 3077
## 10 racecatNon-Hispanic Black 1.3136 0.9797 1.7613 Era 2: 2017-2020 3077
## 11 racecatNon-Hispanic Asian 0.8397 0.5947 1.1855 Era 2: 2017-2020 3077
## 12 racecatOther or multiracial 1.0502 0.5894 1.8715 Era 2: 2017-2020 3077
## 13 educat9th to 11th grade 1.5206 0.9569 2.4163 Era 2: 2017-2020 3077
## 14 educatHigh school or GED 1.2201 0.8432 1.7653 Era 2: 2017-2020 3077
## 15 educatSome college 1.2049 0.7911 1.8350 Era 2: 2017-2020 3077
## 16 educatCollege graduate 0.9199 0.6142 1.3777 Era 2: 2017-2020 3077
## 17 pir 0.9552 0.8760 1.0417 Era 2: 2017-2020 3077
## 18 smokecatFormer 1.1682 0.8944 1.5259 Era 2: 2017-2020 3077
## 19 smokecatCurrent 1.0116 0.8041 1.2726 Era 2: 2017-2020 3077
## 20 srhVery good 2.1886 1.5286 3.1337 Era 2: 2017-2020 3077
## 21 srhGood 3.4091 2.4951 4.6581 Era 2: 2017-2020 3077
## 22 srhFair 5.4757 3.8494 7.7891 Era 2: 2017-2020 3077
## 23 srhPoor 9.4153 4.3838 20.2216 Era 2: 2017-2020 3077
## 24 physactYes 0.7345 0.6010 0.8976 Era 2: 2017-2020 3077
## 25 1|2 15.0540 8.0955 27.9938 Era 2: 2017-2020 3077
## 26 2|3 413.1957 230.9575 739.2299 Era 2: 2017-2020 3077The proportional-odds model, pooled
Research question for this block: what does the association look like when the two instruments are pooled, which is what the literature reports?
fit_po <- svyolr(factor(outcome, ordered = TRUE) ~ drinkcat + RIDAGEYR + sex + racecat +
educat + pir + smokecat + srh + physact, design = des_mod) # the same model, all model rows
or_po <- as.data.frame(broom::tidy(fit_po, conf.int = TRUE, exponentiate = TRUE))[, c("term", "estimate", "conf.low", "conf.high")] # exponentiated coefficients with 95% confidence limits
or_po$model <- rep("Pooled 2011-2020", nrow(or_po)) # label the block for Table 3
or_po$nrows <- rep(nrow(des_mod), nrow(or_po)) # the rows this model used
or_po # every term## term estimate conf.low conf.high model nrows
## 1 drinkcatFormer drinker 1.0911 0.8780 1.3560 Pooled 2011-2020 8842
## 2 drinkcatOccasional 0.8690 0.7241 1.0429 Pooled 2011-2020 8842
## 3 drinkcatLight 0.7067 0.5833 0.8560 Pooled 2011-2020 8842
## 4 drinkcatModerate 1.0239 0.7625 1.3750 Pooled 2011-2020 8842
## 5 drinkcatHeavy 0.8176 0.6106 1.0947 Pooled 2011-2020 8842
## 6 RIDAGEYR 1.0792 1.0744 1.0840 Pooled 2011-2020 8842
## 7 sexFemale 0.5827 0.5046 0.6728 Pooled 2011-2020 8842
## 8 racecatMexican American 0.8660 0.7219 1.0390 Pooled 2011-2020 8842
## 9 racecatOther Hispanic 0.8321 0.6643 1.0424 Pooled 2011-2020 8842
## 10 racecatNon-Hispanic Black 1.2884 1.1127 1.4920 Pooled 2011-2020 8842
## 11 racecatNon-Hispanic Asian 0.8745 0.7212 1.0605 Pooled 2011-2020 8842
## 12 racecatOther or multiracial 1.2413 0.8725 1.7660 Pooled 2011-2020 8842
## 13 educat9th to 11th grade 1.4263 1.0688 1.9033 Pooled 2011-2020 8842
## 14 educatHigh school or GED 1.2142 0.9522 1.5485 Pooled 2011-2020 8842
## 15 educatSome college 1.3579 1.0702 1.7229 Pooled 2011-2020 8842
## 16 educatCollege graduate 0.9934 0.7637 1.2922 Pooled 2011-2020 8842
## 17 pir 0.9091 0.8649 0.9555 Pooled 2011-2020 8842
## 18 smokecatFormer 1.0869 0.9282 1.2727 Pooled 2011-2020 8842
## 19 smokecatCurrent 1.0613 0.9074 1.2414 Pooled 2011-2020 8842
## 20 srhVery good 2.0356 1.6069 2.5786 Pooled 2011-2020 8842
## 21 srhGood 3.1207 2.4884 3.9137 Pooled 2011-2020 8842
## 22 srhFair 5.0913 3.9236 6.6065 Pooled 2011-2020 8842
## 23 srhPoor 8.1221 5.2249 12.6257 Pooled 2011-2020 8842
## 24 physactYes 0.8022 0.7009 0.9182 Pooled 2011-2020 8842
## 25 1|2 20.3477 13.7757 30.0550 Pooled 2011-2020 8842
## 26 2|3 622.4415 414.5413 934.6075 Pooled 2011-2020 8842## Call:
## svyolr(factor(outcome, ordered = TRUE) ~ drinkcat + RIDAGEYR +
## sex + racecat + educat + pir + smokecat + srh + physact,
## design = des_mod)
##
## Coefficients:
## Value Std. Error t value
## drinkcatFormer drinker 0.087203 0.110894 0.78636
## drinkcatOccasional -0.140357 0.093066 -1.50815
## drinkcatLight -0.347218 0.097836 -3.54897
## drinkcatModerate 0.023655 0.150426 0.15725
## drinkcatHeavy -0.201403 0.148905 -1.35257
## RIDAGEYR 0.076233 0.002272 33.55411
## sexFemale -0.540082 0.073375 -7.36055
## racecatMexican American -0.143825 0.092898 -1.54820
## racecatOther Hispanic -0.183762 0.114963 -1.59844
## racecatNon-Hispanic Black 0.253429 0.074838 3.38635
## racecatNon-Hispanic Asian -0.134075 0.098375 -1.36291
## racecatOther or multiracial 0.216177 0.179879 1.20179
## educat9th to 11th grade 0.355063 0.147209 2.41197
## educatHigh school or GED 0.194112 0.124055 1.56473
## educatSome college 0.305946 0.121477 2.51855
## educatCollege graduate -0.006644 0.134166 -0.04952
## pir -0.095304 0.025405 -3.75137
## smokecatFormer 0.083331 0.080517 1.03494
## smokecatCurrent 0.059526 0.079944 0.74460
## srhVery good 0.710790 0.120647 5.89149
## srhGood 1.138052 0.115524 9.85126
## srhFair 1.627530 0.132919 12.24448
## srhPoor 2.094590 0.225080 9.30598
## physactYes -0.220363 0.068888 -3.19886
##
## Intercepts:
## Value Std. Error t value
## 1|2 3.013 0.199 15.139
## 2|3 6.434 0.207 31.022Sensitivity: the same models on the five-stage outcome, era 1
Research question for this block: how much of this arm’s answer depends on the collapse from five stages to three tiers that the cell-check rule forced?
fs_e1 <- svyolr(factor(ckmstage, ordered = TRUE) ~ drinkcat + RIDAGEYR + sex + racecat +
educat + pir + smokecat + srh + physact, design = des_e1) # era 1, five-stage outcome
os_e1 <- as.data.frame(broom::tidy(fs_e1, conf.int = TRUE, exponentiate = TRUE))[, c("term", "estimate", "conf.low", "conf.high")] # exponentiated coefficients with 95% confidence limits
os_e1$model <- rep("Era 1: 2011-2016", nrow(os_e1)) # label the block for the companion table
os_e1$nrows <- rep(nrow(des_e1), nrow(os_e1)) # the rows this model used
os_e1 # every term## term estimate conf.low conf.high model nrows
## 1 drinkcatFormer drinker 1.2359 0.9746 1.5673 Era 1: 2011-2016 5765
## 2 drinkcatOccasional 0.9361 0.7866 1.1142 Era 1: 2011-2016 5765
## 3 drinkcatLight 0.7025 0.5671 0.8702 Era 1: 2011-2016 5765
## 4 drinkcatModerate 0.9600 0.7063 1.3049 Era 1: 2011-2016 5765
## 5 drinkcatHeavy 0.9434 0.6863 1.2969 Era 1: 2011-2016 5765
## 6 RIDAGEYR 1.0821 1.0755 1.0887 Era 1: 2011-2016 5765
## 7 sexFemale 0.5366 0.4653 0.6187 Era 1: 2011-2016 5765
## 8 racecatMexican American 1.0078 0.8347 1.2168 Era 1: 2011-2016 5765
## 9 racecatOther Hispanic 0.9475 0.7546 1.1897 Era 1: 2011-2016 5765
## 10 racecatNon-Hispanic Black 1.2865 1.0993 1.5055 Era 1: 2011-2016 5765
## 11 racecatNon-Hispanic Asian 0.8591 0.7103 1.0390 Era 1: 2011-2016 5765
## 12 racecatOther or multiracial 1.4399 0.9093 2.2801 Era 1: 2011-2016 5765
## 13 educat9th to 11th grade 1.4468 1.0537 1.9866 Era 1: 2011-2016 5765
## 14 educatHigh school or GED 1.3051 1.0073 1.6910 Era 1: 2011-2016 5765
## 15 educatSome college 1.5504 1.2317 1.9516 Era 1: 2011-2016 5765
## 16 educatCollege graduate 1.1066 0.8269 1.4808 Era 1: 2011-2016 5765
## 17 pir 0.8888 0.8365 0.9444 Era 1: 2011-2016 5765
## 18 smokecatFormer 1.0668 0.9035 1.2596 Era 1: 2011-2016 5765
## 19 smokecatCurrent 1.1337 0.9405 1.3665 Era 1: 2011-2016 5765
## 20 srhVery good 2.0286 1.5463 2.6614 Era 1: 2011-2016 5765
## 21 srhGood 3.2397 2.4746 4.2413 Era 1: 2011-2016 5765
## 22 srhFair 5.2621 3.7893 7.3075 Era 1: 2011-2016 5765
## 23 srhPoor 8.8647 5.4116 14.5211 Era 1: 2011-2016 5765
## 24 physactYes 0.8521 0.7207 1.0074 Era 1: 2011-2016 5765
## 25 0|1 3.2652 2.0591 5.1777 Era 1: 2011-2016 5765
## 26 1|2 26.3200 16.2027 42.7546 Era 1: 2011-2016 5765
## 27 2|3 881.4178 514.4747 1510.0788 Era 1: 2011-2016 5765
## 28 3|4 1743.7931 1006.5342 3021.0741 Era 1: 2011-2016 5765Sensitivity on the five-stage outcome, era 2
Research question for this block: and the same question inside the second instrument era?
fs_e2 <- try(svyolr(factor(ckmstage, ordered = TRUE) ~ drinkcat + RIDAGEYR + sex + racecat +
educat + pir + smokecat + srh + physact, design = des_e2), silent = TRUE) # guarded for the same reason as the era 2 primary model
if (inherits(fs_e2, "try-error")) cat("Era 2 five-stage model did not fit on this file:", conditionMessage(attr(fs_e2, "condition")), "\n") # the reason travels with the failure
os_e2 <- data.frame(term = character(0), estimate = numeric(0), conf.low = numeric(0), conf.high = numeric(0)) # stays empty if the fit failed
if (!inherits(fs_e2, "try-error")) os_e2 <- as.data.frame(broom::tidy(fs_e2, conf.int = TRUE, exponentiate = TRUE))[, c("term", "estimate", "conf.low", "conf.high")] # exponentiated coefficients with 95% limits, kept only if the fit succeeded
os_e2$model <- rep("Era 2: 2017-2020", nrow(os_e2)) # label the block for the companion table
os_e2$nrows <- rep(nrow(des_e2), nrow(os_e2)) # the rows this model used
os_e2 # every term## term estimate conf.low conf.high model nrows
## 1 drinkcatFormer drinker 0.9682 0.5898 1.5895 Era 2: 2017-2020 3077
## 2 drinkcatOccasional 0.8638 0.5400 1.3818 Era 2: 2017-2020 3077
## 3 drinkcatLight 0.8042 0.5484 1.1792 Era 2: 2017-2020 3077
## 4 drinkcatModerate 1.3428 0.7729 2.3329 Era 2: 2017-2020 3077
## 5 drinkcatHeavy 0.6245 0.3337 1.1687 Era 2: 2017-2020 3077
## 6 RIDAGEYR 1.0688 1.0637 1.0739 Era 2: 2017-2020 3077
## 7 sexFemale 0.6258 0.4966 0.7886 Era 2: 2017-2020 3077
## 8 racecatMexican American 0.9352 0.6768 1.2922 Era 2: 2017-2020 3077
## 9 racecatOther Hispanic 0.7669 0.5420 1.0852 Era 2: 2017-2020 3077
## 10 racecatNon-Hispanic Black 1.3948 1.0668 1.8237 Era 2: 2017-2020 3077
## 11 racecatNon-Hispanic Asian 0.9634 0.6748 1.3753 Era 2: 2017-2020 3077
## 12 racecatOther or multiracial 1.0446 0.6227 1.7523 Era 2: 2017-2020 3077
## 13 educat9th to 11th grade 1.4187 0.8825 2.2806 Era 2: 2017-2020 3077
## 14 educatHigh school or GED 1.2381 0.8446 1.8149 Era 2: 2017-2020 3077
## 15 educatSome college 1.2149 0.7992 1.8471 Era 2: 2017-2020 3077
## 16 educatCollege graduate 1.0004 0.6658 1.5034 Era 2: 2017-2020 3077
## 17 pir 0.9553 0.8752 1.0427 Era 2: 2017-2020 3077
## 18 smokecatFormer 1.2014 0.9511 1.5175 Era 2: 2017-2020 3077
## 19 smokecatCurrent 1.0269 0.8041 1.3114 Era 2: 2017-2020 3077
## 20 srhVery good 2.4958 1.7970 3.4662 Era 2: 2017-2020 3077
## 21 srhGood 4.0789 2.9557 5.6289 Era 2: 2017-2020 3077
## 22 srhFair 6.6808 4.7387 9.4187 Era 2: 2017-2020 3077
## 23 srhPoor 11.7355 5.3177 25.8987 Era 2: 2017-2020 3077
## 24 physactYes 0.7524 0.6171 0.9173 Era 2: 2017-2020 3077
## 25 0|1 2.4609 1.2809 4.7276 Era 2: 2017-2020 3077
## 26 1|2 19.5600 9.8626 38.7925 Era 2: 2017-2020 3077
## 27 2|3 537.7007 279.1524 1035.7142 Era 2: 2017-2020 3077
## 28 3|4 961.5256 497.7648 1857.3663 Era 2: 2017-2020 3077Sensitivity on the five-stage outcome, pooled
Research question for this block: and pooled, which is the comparison a reader will make against the primary table?
fs_po <- svyolr(factor(ckmstage, ordered = TRUE) ~ drinkcat + RIDAGEYR + sex + racecat +
educat + pir + smokecat + srh + physact, design = des_mod) # pooled, five-stage outcome
os_po <- as.data.frame(broom::tidy(fs_po, conf.int = TRUE, exponentiate = TRUE))[, c("term", "estimate", "conf.low", "conf.high")] # exponentiated coefficients with 95% confidence limits
os_po$model <- rep("Pooled 2011-2020", nrow(os_po)) # label the block for the companion table
os_po$nrows <- rep(nrow(des_mod), nrow(os_po)) # the rows this model used
os_po # every term## term estimate conf.low conf.high model nrows
## 1 drinkcatFormer drinker 1.1283 0.9080 1.4021 Pooled 2011-2020 8842
## 2 drinkcatOccasional 0.9178 0.7667 1.0986 Pooled 2011-2020 8842
## 3 drinkcatLight 0.7456 0.6267 0.8870 Pooled 2011-2020 8842
## 4 drinkcatModerate 1.0826 0.8238 1.4228 Pooled 2011-2020 8842
## 5 drinkcatHeavy 0.8047 0.5991 1.0808 Pooled 2011-2020 8842
## 6 RIDAGEYR 1.0765 1.0723 1.0808 Pooled 2011-2020 8842
## 7 sexFemale 0.5699 0.5024 0.6464 Pooled 2011-2020 8842
## 8 racecatMexican American 0.9860 0.8369 1.1616 Pooled 2011-2020 8842
## 9 racecatOther Hispanic 0.8703 0.7165 1.0573 Pooled 2011-2020 8842
## 10 racecatNon-Hispanic Black 1.3251 1.1542 1.5214 Pooled 2011-2020 8842
## 11 racecatNon-Hispanic Asian 0.8963 0.7506 1.0702 Pooled 2011-2020 8842
## 12 racecatOther or multiracial 1.2381 0.8812 1.7393 Pooled 2011-2020 8842
## 13 educat9th to 11th grade 1.4191 1.0905 1.8469 Pooled 2011-2020 8842
## 14 educatHigh school or GED 1.2685 1.0222 1.5742 Pooled 2011-2020 8842
## 15 educatSome college 1.4000 1.1331 1.7298 Pooled 2011-2020 8842
## 16 educatCollege graduate 1.0535 0.8315 1.3347 Pooled 2011-2020 8842
## 17 pir 0.9154 0.8710 0.9622 Pooled 2011-2020 8842
## 18 smokecatFormer 1.1281 0.9849 1.2921 Pooled 2011-2020 8842
## 19 smokecatCurrent 1.0944 0.9422 1.2713 Pooled 2011-2020 8842
## 20 srhVery good 2.1770 1.7652 2.6850 Pooled 2011-2020 8842
## 21 srhGood 3.4718 2.8126 4.2856 Pooled 2011-2020 8842
## 22 srhFair 5.6742 4.4334 7.2622 Pooled 2011-2020 8842
## 23 srhPoor 9.7489 6.3317 15.0103 Pooled 2011-2020 8842
## 24 physactYes 0.8216 0.7224 0.9343 Pooled 2011-2020 8842
## 25 0|1 2.9133 2.0005 4.2427 Pooled 2011-2020 8842
## 26 1|2 23.1815 15.6399 34.3597 Pooled 2011-2020 8842
## 27 2|3 711.3437 472.2931 1071.3894 Pooled 2011-2020 8842
## 28 3|4 1351.1460 891.4299 2047.9409 Pooled 2011-2020 8842Check: the binary high-risk outcome
Research question for this block: and if the outcome is only whether a person is at stage 3 or 4, which is the clinically loaded end?
des_mod <- update(des_mod, highrisk = as.integer(ckmstage >= 3)) # 1 for stages 3 and 4, 0 for stages 0 to 2
table(des_mod$variables$highrisk) # unweighted counts behind the binary check##
## 0 1
## 7139 1703fit_hr <- svyglm(highrisk ~ drinkcat + RIDAGEYR + sex + racecat + educat + pir +
smokecat + srh + physact, design = des_mod, family = quasibinomial()) # survey logistic regression
or_hr <- as.data.frame(broom::tidy(fit_hr, conf.int = TRUE, exponentiate = TRUE))[, c("term", "estimate", "conf.low", "conf.high")] # exponentiated coefficients with 95% confidence limits
or_hr$model <- rep("Check: high risk, stages 3 to 4, pooled", nrow(or_hr)) # label the block for Table 3
or_hr$nrows <- rep(nrow(des_mod), nrow(or_hr)) # the rows this model used
or_hr # every term## term estimate conf.low conf.high
## 1 (Intercept) 0.0001143 4.725e-05 2.764e-04
## 2 drinkcatFormer drinker 0.9292347 6.939e-01 1.244e+00
## 3 drinkcatOccasional 0.8213573 6.015e-01 1.122e+00
## 4 drinkcatLight 0.5374299 3.527e-01 8.188e-01
## 5 drinkcatModerate 0.7489528 4.058e-01 1.382e+00
## 6 drinkcatHeavy 0.4320975 2.539e-01 7.355e-01
## 7 RIDAGEYR 1.1327091 1.123e+00 1.143e+00
## 8 sexFemale 0.4967071 3.902e-01 6.323e-01
## 9 racecatMexican American 0.5515584 4.099e-01 7.422e-01
## 10 racecatOther Hispanic 0.9015945 6.410e-01 1.268e+00
## 11 racecatNon-Hispanic Black 1.3010192 1.031e+00 1.641e+00
## 12 racecatNon-Hispanic Asian 0.5043932 3.588e-01 7.091e-01
## 13 racecatOther or multiracial 1.4439351 8.247e-01 2.528e+00
## 14 educat9th to 11th grade 1.4939783 1.005e+00 2.220e+00
## 15 educatHigh school or GED 1.2953715 8.873e-01 1.891e+00
## 16 educatSome college 1.2323411 8.649e-01 1.756e+00
## 17 educatCollege graduate 1.0183261 7.416e-01 1.398e+00
## 18 pir 0.8881717 8.311e-01 9.492e-01
## 19 smokecatFormer 1.4429064 1.186e+00 1.756e+00
## 20 smokecatCurrent 1.7441818 1.330e+00 2.288e+00
## 21 srhVery good 1.3307338 8.699e-01 2.036e+00
## 22 srhGood 2.1002551 1.477e+00 2.986e+00
## 23 srhFair 4.3993052 3.069e+00 6.307e+00
## 24 srhPoor 8.5575035 4.806e+00 1.524e+01
## 25 physactYes 0.8622694 7.045e-01 1.055e+00
## model nrows
## 1 Check: high risk, stages 3 to 4, pooled 8842
## 2 Check: high risk, stages 3 to 4, pooled 8842
## 3 Check: high risk, stages 3 to 4, pooled 8842
## 4 Check: high risk, stages 3 to 4, pooled 8842
## 5 Check: high risk, stages 3 to 4, pooled 8842
## 6 Check: high risk, stages 3 to 4, pooled 8842
## 7 Check: high risk, stages 3 to 4, pooled 8842
## 8 Check: high risk, stages 3 to 4, pooled 8842
## 9 Check: high risk, stages 3 to 4, pooled 8842
## 10 Check: high risk, stages 3 to 4, pooled 8842
## 11 Check: high risk, stages 3 to 4, pooled 8842
## 12 Check: high risk, stages 3 to 4, pooled 8842
## 13 Check: high risk, stages 3 to 4, pooled 8842
## 14 Check: high risk, stages 3 to 4, pooled 8842
## 15 Check: high risk, stages 3 to 4, pooled 8842
## 16 Check: high risk, stages 3 to 4, pooled 8842
## 17 Check: high risk, stages 3 to 4, pooled 8842
## 18 Check: high risk, stages 3 to 4, pooled 8842
## 19 Check: high risk, stages 3 to 4, pooled 8842
## 20 Check: high risk, stages 3 to 4, pooled 8842
## 21 Check: high risk, stages 3 to 4, pooled 8842
## 22 Check: high risk, stages 3 to 4, pooled 8842
## 23 Check: high risk, stages 3 to 4, pooled 8842
## 24 Check: high risk, stages 3 to 4, pooled 8842
## 25 Check: high risk, stages 3 to 4, pooled 8842The pre-specified sensitivity analysis: occasional drinkers as the reference
Research question for this block: if the reference is people who drink one drink a week or fewer rather than people who never drank, does the gradient survive?
The arm specification pre-specifies this analysis. It is the corrected reference in Sarich (2024) and Srivatsa (2026), and it answers a different objection from the one arm 1 answers. A lifetime abstainer may differ from a drinker in ways no covariate reaches; an occasional drinker is a drinker. Only the reference level changes, so the rows, the covariates, and the outcome are the same.
des_oc <- update(des_mod, drinkoc = relevel(drinkcat, ref = "Occasional")) # the same variable, occasional drinkers first
table(des_oc$variables$drinkoc) # the reference level is now Occasional##
## Occasional Lifetime abstainer Former drinker Light Moderate
## 3082 1108 1591 2019 575
## Heavy
## 467fit_oc <- svyolr(factor(outcome, ordered = TRUE) ~ drinkoc + RIDAGEYR + sex + racecat +
educat + pir + smokecat + srh + physact, design = des_oc) # the same model, new reference
or_oc <- as.data.frame(broom::tidy(fit_oc, conf.int = TRUE, exponentiate = TRUE))[, c("term", "estimate", "conf.low", "conf.high")] # exponentiated coefficients with 95% confidence limits
or_oc$term <- sub("^drinkoc", "drinkcat", or_oc$term) # one prefix, so Table 3 can hold both references
or_oc$model <- rep("Sensitivity: occasional drinkers as the reference, pooled", nrow(or_oc)) # label the block for Table 3
or_oc$nrows <- rep(nrow(des_mod), nrow(or_oc)) # the rows this model used
or_oc # every term## term estimate conf.low conf.high
## 1 drinkcatLifetime abstainer 1.1501 0.9584 1.3802
## 2 drinkcatFormer drinker 1.2560 1.0527 1.4987
## 3 drinkcatLight 0.8136 0.7005 0.9450
## 4 drinkcatModerate 1.1788 0.9084 1.5297
## 5 drinkcatHeavy 0.9414 0.7032 1.2602
## 6 RIDAGEYR 1.0792 1.0744 1.0840
## 7 sexFemale 0.5826 0.5046 0.6727
## 8 racecatMexican American 0.8662 0.7220 1.0391
## 9 racecatOther Hispanic 0.8323 0.6644 1.0426
## 10 racecatNon-Hispanic Black 1.2886 1.1128 1.4921
## 11 racecatNon-Hispanic Asian 0.8749 0.7215 1.0609
## 12 racecatOther or multiracial 1.2412 0.8725 1.7658
## 13 educat9th to 11th grade 1.4270 1.0694 1.9043
## 14 educatHigh school or GED 1.2150 0.9527 1.5494
## 15 educatSome college 1.3589 1.0710 1.7241
## 16 educatCollege graduate 0.9936 0.7639 1.2925
## 17 pir 0.9092 0.8651 0.9556
## 18 smokecatFormer 1.0869 0.9283 1.2727
## 19 smokecatCurrent 1.0612 0.9073 1.2412
## 20 srhVery good 2.0354 1.6068 2.5783
## 21 srhGood 3.1227 2.4900 3.9161
## 22 srhFair 5.0928 3.9249 6.6084
## 23 srhPoor 8.1325 5.2316 12.6420
## 24 physactYes 0.8024 0.7011 0.9184
## 25 1|2 23.4404 15.5050 35.4371
## 26 2|3 716.6429 471.5678 1089.0841
## model nrows
## 1 Sensitivity: occasional drinkers as the reference, pooled 8842
## 2 Sensitivity: occasional drinkers as the reference, pooled 8842
## 3 Sensitivity: occasional drinkers as the reference, pooled 8842
## 4 Sensitivity: occasional drinkers as the reference, pooled 8842
## 5 Sensitivity: occasional drinkers as the reference, pooled 8842
## 6 Sensitivity: occasional drinkers as the reference, pooled 8842
## 7 Sensitivity: occasional drinkers as the reference, pooled 8842
## 8 Sensitivity: occasional drinkers as the reference, pooled 8842
## 9 Sensitivity: occasional drinkers as the reference, pooled 8842
## 10 Sensitivity: occasional drinkers as the reference, pooled 8842
## 11 Sensitivity: occasional drinkers as the reference, pooled 8842
## 12 Sensitivity: occasional drinkers as the reference, pooled 8842
## 13 Sensitivity: occasional drinkers as the reference, pooled 8842
## 14 Sensitivity: occasional drinkers as the reference, pooled 8842
## 15 Sensitivity: occasional drinkers as the reference, pooled 8842
## 16 Sensitivity: occasional drinkers as the reference, pooled 8842
## 17 Sensitivity: occasional drinkers as the reference, pooled 8842
## 18 Sensitivity: occasional drinkers as the reference, pooled 8842
## 19 Sensitivity: occasional drinkers as the reference, pooled 8842
## 20 Sensitivity: occasional drinkers as the reference, pooled 8842
## 21 Sensitivity: occasional drinkers as the reference, pooled 8842
## 22 Sensitivity: occasional drinkers as the reference, pooled 8842
## 23 Sensitivity: occasional drinkers as the reference, pooled 8842
## 24 Sensitivity: occasional drinkers as the reference, pooled 8842
## 25 Sensitivity: occasional drinkers as the reference, pooled 8842
## 26 Sensitivity: occasional drinkers as the reference, pooled 8842The other two arms’ references, for Figure 1
Research question for this block: how far do the same four current-drinking categories move when only the reference group changes?
Arm 2 pools lifetime abstainers with former drinkers; arm 3 splits the former-drinker row on whether the person ever drank four or five or more drinks every day. The two models are fitted here, on the same rows and with the same covariates and the same outcome, so that Figure 1 can put the three reference constructions beside each other. This is the only place where this arm touches another arm’s exposure, and it does so by writing the model out rather than by reading another file.
fit_a2 <- svyolr(factor(outcome, ordered = TRUE) ~ nondrink + RIDAGEYR + sex + racecat +
educat + pir + smokecat + srh + physact, design = des_mod) # arm 2: the pooled non-drinker reference
or_a2 <- as.data.frame(broom::tidy(fit_a2, conf.int = TRUE, exponentiate = TRUE))[, c("term", "estimate", "conf.low", "conf.high")] # exponentiated coefficients with 95% confidence limits
or_a2$model <- rep("Arm 2: pooled non-drinker reference", nrow(or_a2)) # label the block for Figure 1
or_a2$nrows <- rep(nrow(des_mod), nrow(or_a2)) # the rows this model used
des_re <- subset(des_mod, !is.na(drinkre)) # arm 3's exposure is missing for former drinkers with no ALQ151 answer
nrow(des_re) # rows the arm 3 model uses## [1] 8841fit_a3 <- svyolr(factor(outcome, ordered = TRUE) ~ drinkre + RIDAGEYR + sex + racecat +
educat + pir + smokecat + srh + physact, design = des_re) # arm 3: former drinkers split on past heavy drinking
or_a3 <- as.data.frame(broom::tidy(fit_a3, conf.int = TRUE, exponentiate = TRUE))[, c("term", "estimate", "conf.low", "conf.high")] # exponentiated coefficients with 95% confidence limits
or_a3$model <- rep("Arm 3: former drinkers split on past heavy drinking", nrow(or_a3)) # label the block for Figure 1
or_a3$nrows <- rep(nrow(des_re), nrow(or_a3)) # the rows this model used
or_a2 # arm 2's terms## term estimate conf.low conf.high model
## 1 nondrinkOccasional 0.8239 0.7131 0.9521 Arm 2: pooled non-drinker reference
## 2 nondrinkLight 0.6696 0.5669 0.7909 Arm 2: pooled non-drinker reference
## 3 nondrinkModerate 0.9682 0.7202 1.3015 Arm 2: pooled non-drinker reference
## 4 nondrinkHeavy 0.7721 0.5661 1.0530 Arm 2: pooled non-drinker reference
## 5 RIDAGEYR 1.0793 1.0744 1.0841 Arm 2: pooled non-drinker reference
## 6 sexFemale 0.5814 0.5027 0.6724 Arm 2: pooled non-drinker reference
## 7 racecatMexican American 0.8661 0.7220 1.0391 Arm 2: pooled non-drinker reference
## 8 racecatOther Hispanic 0.8318 0.6640 1.0421 Arm 2: pooled non-drinker reference
## 9 racecatNon-Hispanic Black 1.2872 1.1108 1.4917 Arm 2: pooled non-drinker reference
## 10 racecatNon-Hispanic Asian 0.8671 0.7164 1.0495 Arm 2: pooled non-drinker reference
## 11 racecatOther or multiracial 1.2431 0.8748 1.7665 Arm 2: pooled non-drinker reference
## 12 educat9th to 11th grade 1.4329 1.0742 1.9114 Arm 2: pooled non-drinker reference
## 13 educatHigh school or GED 1.2200 0.9551 1.5583 Arm 2: pooled non-drinker reference
## 14 educatSome college 1.3651 1.0747 1.7340 Arm 2: pooled non-drinker reference
## 15 educatCollege graduate 0.9991 0.7694 1.2974 Arm 2: pooled non-drinker reference
## 16 pir 0.9095 0.8654 0.9560 Arm 2: pooled non-drinker reference
## 17 smokecatFormer 1.0971 0.9398 1.2807 Arm 2: pooled non-drinker reference
## 18 smokecatCurrent 1.0707 0.9124 1.2565 Arm 2: pooled non-drinker reference
## 19 srhVery good 2.0365 1.6074 2.5801 Arm 2: pooled non-drinker reference
## 20 srhGood 3.1291 2.4953 3.9238 Arm 2: pooled non-drinker reference
## 21 srhFair 5.1026 3.9309 6.6236 Arm 2: pooled non-drinker reference
## 22 srhPoor 8.1754 5.2630 12.6994 Arm 2: pooled non-drinker reference
## 23 physactYes 0.8027 0.7014 0.9187 Arm 2: pooled non-drinker reference
## 24 1|2 19.5463 13.2750 28.7803 Arm 2: pooled non-drinker reference
## 25 2|3 597.2912 402.9775 885.3020 Arm 2: pooled non-drinker reference
## nrows
## 1 8842
## 2 8842
## 3 8842
## 4 8842
## 5 8842
## 6 8842
## 7 8842
## 8 8842
## 9 8842
## 10 8842
## 11 8842
## 12 8842
## 13 8842
## 14 8842
## 15 8842
## 16 8842
## 17 8842
## 18 8842
## 19 8842
## 20 8842
## 21 8842
## 22 8842
## 23 8842
## 24 8842
## 25 8842## term estimate conf.low conf.high
## 1 drinkreFormer, not past heavy 1.2177 0.9949 1.4903
## 2 drinkreFormer, past heavy 0.7673 0.4532 1.2993
## 3 drinkreOccasional 0.8630 0.7191 1.0357
## 4 drinkreLight 0.6988 0.5759 0.8479
## 5 drinkreModerate 1.0092 0.7563 1.3466
## 6 drinkreHeavy 0.8014 0.5987 1.0726
## 7 RIDAGEYR 1.0791 1.0743 1.0839
## 8 sexFemale 0.5738 0.4990 0.6598
## 9 racecatMexican American 0.8669 0.7221 1.0407
## 10 racecatOther Hispanic 0.8270 0.6592 1.0375
## 11 racecatNon-Hispanic Black 1.2827 1.1115 1.4803
## 12 racecatNon-Hispanic Asian 0.8641 0.7121 1.0486
## 13 racecatOther or multiracial 1.2246 0.8660 1.7316
## 14 educat9th to 11th grade 1.4302 1.0690 1.9134
## 15 educatHigh school or GED 1.2190 0.9580 1.5510
## 16 educatSome college 1.3592 1.0710 1.7248
## 17 educatCollege graduate 0.9958 0.7651 1.2959
## 18 pir 0.9092 0.8651 0.9555
## 19 smokecatFormer 1.1047 0.9375 1.3019
## 20 smokecatCurrent 1.0794 0.9233 1.2618
## 21 srhVery good 2.0288 1.6029 2.5678
## 22 srhGood 3.1256 2.4938 3.9176
## 23 srhFair 5.1078 3.9400 6.6217
## 24 srhPoor 8.2848 5.4238 12.6549
## 25 physactYes 0.8040 0.7026 0.9200
## 26 1|2 20.0784 13.5589 29.7327
## 27 2|3 616.4028 409.7471 927.2854
## model nrows
## 1 Arm 3: former drinkers split on past heavy drinking 8841
## 2 Arm 3: former drinkers split on past heavy drinking 8841
## 3 Arm 3: former drinkers split on past heavy drinking 8841
## 4 Arm 3: former drinkers split on past heavy drinking 8841
## 5 Arm 3: former drinkers split on past heavy drinking 8841
## 6 Arm 3: former drinkers split on past heavy drinking 8841
## 7 Arm 3: former drinkers split on past heavy drinking 8841
## 8 Arm 3: former drinkers split on past heavy drinking 8841
## 9 Arm 3: former drinkers split on past heavy drinking 8841
## 10 Arm 3: former drinkers split on past heavy drinking 8841
## 11 Arm 3: former drinkers split on past heavy drinking 8841
## 12 Arm 3: former drinkers split on past heavy drinking 8841
## 13 Arm 3: former drinkers split on past heavy drinking 8841
## 14 Arm 3: former drinkers split on past heavy drinking 8841
## 15 Arm 3: former drinkers split on past heavy drinking 8841
## 16 Arm 3: former drinkers split on past heavy drinking 8841
## 17 Arm 3: former drinkers split on past heavy drinking 8841
## 18 Arm 3: former drinkers split on past heavy drinking 8841
## 19 Arm 3: former drinkers split on past heavy drinking 8841
## 20 Arm 3: former drinkers split on past heavy drinking 8841
## 21 Arm 3: former drinkers split on past heavy drinking 8841
## 22 Arm 3: former drinkers split on past heavy drinking 8841
## 23 Arm 3: former drinkers split on past heavy drinking 8841
## 24 Arm 3: former drinkers split on past heavy drinking 8841
## 25 Arm 3: former drinkers split on past heavy drinking 8841
## 26 Arm 3: former drinkers split on past heavy drinking 8841
## 27 Arm 3: former drinkers split on past heavy drinking 8841Part 7: Tables
Table 2: the weighted stage distribution by drinking category and era
Research question for this block: what does the stage distribution look like written out with its confidence limits?
t2_p <- data.frame( # the pooled block, one row per category and stage
era = rep("Pooled 2011-2020", nrow(sb_p) * 5), # which rows this block describes
category = rep(as.character(sb_p$drinkcat), times = 5), # the drinking category
stage = rep(0:4, each = nrow(sb_p)), # the CKM stage, 0 to 4
pct = round(100 * as.vector(as.matrix(sb_p[, 2:6])), 2), # weighted percentage at that stage
lcl = round(100 * ci_p[, 1], 2), # lower 95% confidence limit
ucl = round(100 * ci_p[, 2], 2)) # upper 95% confidence limit
t2_1 <- data.frame( # the era 1 block
era = rep("Era 1: 2011-2016", nrow(sb_1) * 5), # which rows this block describes
category = rep(as.character(sb_1$drinkcat), times = 5), # the drinking category
stage = rep(0:4, each = nrow(sb_1)), # the CKM stage, 0 to 4
pct = round(100 * as.vector(as.matrix(sb_1[, 2:6])), 2), # weighted percentage at that stage
lcl = round(100 * ci_1[, 1], 2), # lower 95% confidence limit
ucl = round(100 * ci_1[, 2], 2)) # upper 95% confidence limit
t2_2 <- data.frame( # the era 2 block
era = rep("Era 2: 2017-2020", nrow(sb_2) * 5), # which rows this block describes
category = rep(as.character(sb_2$drinkcat), times = 5), # the drinking category
stage = rep(0:4, each = nrow(sb_2)), # the CKM stage, 0 to 4
pct = round(100 * as.vector(as.matrix(sb_2[, 2:6])), 2), # weighted percentage at that stage
lcl = round(100 * ci_2[, 1], 2), # lower 95% confidence limit
ucl = round(100 * ci_2[, 2], 2)) # upper 95% confidence limit
n_1 <- as.data.frame(table(des_x1$variables$drinkcat, des_x1$variables$ckmstage)) # era 1 unweighted cell counts
n_1$era <- rep("Era 1: 2011-2016", nrow(n_1)) # label the block
n_2 <- as.data.frame(table(des_x2$variables$drinkcat, des_x2$variables$ckmstage)) # era 2 unweighted cell counts
n_2$era <- rep("Era 2: 2017-2020", nrow(n_2)) # label the block
n_0 <- as.data.frame(table(des_exp$variables$drinkcat, des_exp$variables$ckmstage)) # pooled unweighted cell counts
n_0$era <- rep("Pooled 2011-2020", nrow(n_0)) # label the block
n_cell <- rbind(n_1, n_2, n_0) # the three blocks of counts
names(n_cell) <- c("category", "stage", "n_unwt", "era") # name them to match the join
n_cell$category <- as.character(n_cell$category) # character, so the join keys have one type
n_cell$stage <- as.integer(as.character(n_cell$stage)) # integer, so the join keys have one type
table2 <- rbind(t2_1, t2_2, t2_p) # era first, then pooled, as every specification asks
table2 <- dplyr::left_join(table2, n_cell, by = dplyr::join_by(era, category, stage)) # the unweighted count behind each percentage
table2 # the table as it will be written## era category stage pct lcl ucl n_unwt
## 1 Era 1: 2011-2016 Lifetime abstainer 0 9.61 6.82 12.41 76
## 2 Era 1: 2011-2016 Former drinker 0 4.18 2.18 6.19 31
## 3 Era 1: 2011-2016 Occasional 0 8.60 6.63 10.57 176
## 4 Era 1: 2011-2016 Light 0 14.83 12.52 17.14 165
## 5 Era 1: 2011-2016 Moderate 0 8.54 4.86 12.23 32
## 6 Era 1: 2011-2016 Heavy 0 8.04 4.12 11.97 20
## 7 Era 1: 2011-2016 Lifetime abstainer 1 23.03 20.05 26.02 189
## 8 Era 1: 2011-2016 Former drinker 1 15.70 12.70 18.71 142
## 9 Era 1: 2011-2016 Occasional 1 27.15 24.53 29.76 546
## 10 Era 1: 2011-2016 Light 1 30.64 27.27 34.02 392
## 11 Era 1: 2011-2016 Moderate 1 26.35 20.20 32.50 100
## 12 Era 1: 2011-2016 Heavy 1 14.64 9.49 19.78 46
## 13 Era 1: 2011-2016 Lifetime abstainer 2 47.24 43.13 51.35 476
## 14 Era 1: 2011-2016 Former drinker 2 49.73 45.59 53.86 511
## 15 Era 1: 2011-2016 Occasional 2 50.64 48.02 53.25 1051
## 16 Era 1: 2011-2016 Light 2 47.14 43.42 50.85 680
## 17 Era 1: 2011-2016 Moderate 2 53.09 47.69 58.49 232
## 18 Era 1: 2011-2016 Heavy 2 66.81 60.70 72.92 198
## 19 Era 1: 2011-2016 Lifetime abstainer 3 9.22 6.74 11.70 109
## 20 Era 1: 2011-2016 Former drinker 3 10.08 8.31 11.84 162
## 21 Era 1: 2011-2016 Occasional 3 5.02 3.93 6.11 147
## 22 Era 1: 2011-2016 Light 3 2.77 1.64 3.90 65
## 23 Era 1: 2011-2016 Moderate 3 2.84 1.60 4.07 24
## 24 Era 1: 2011-2016 Heavy 3 2.89 0.35 5.42 10
## 25 Era 1: 2011-2016 Lifetime abstainer 4 10.89 8.30 13.48 120
## 26 Era 1: 2011-2016 Former drinker 4 20.31 17.12 23.50 222
## 27 Era 1: 2011-2016 Occasional 4 8.60 7.19 10.01 214
## 28 Era 1: 2011-2016 Light 4 4.62 3.26 5.98 81
## 29 Era 1: 2011-2016 Moderate 4 9.18 5.89 12.46 42
## 30 Era 1: 2011-2016 Heavy 4 7.62 3.84 11.40 26
## 31 Era 2: 2017-2020 Lifetime abstainer 0 8.19 3.64 12.75 13
## 32 Era 2: 2017-2020 Former drinker 0 4.12 1.97 6.27 24
## 33 Era 2: 2017-2020 Occasional 0 10.58 6.64 14.51 101
## 34 Era 2: 2017-2020 Light 0 10.59 6.87 14.30 63
## 35 Era 2: 2017-2020 Moderate 0 5.75 0.25 11.25 8
## 36 Era 2: 2017-2020 Heavy 0 10.45 2.90 18.00 11
## 37 Era 2: 2017-2020 Lifetime abstainer 1 18.98 14.06 23.90 69
## 38 Era 2: 2017-2020 Former drinker 1 16.10 12.46 19.73 103
## 39 Era 2: 2017-2020 Occasional 1 25.99 21.41 30.57 292
## 40 Era 2: 2017-2020 Light 1 34.81 29.29 40.34 233
## 41 Era 2: 2017-2020 Moderate 1 19.28 13.17 25.39 37
## 42 Era 2: 2017-2020 Heavy 1 22.59 13.59 31.58 38
## 43 Era 2: 2017-2020 Lifetime abstainer 2 54.74 46.93 62.56 177
## 44 Era 2: 2017-2020 Former drinker 2 49.68 44.59 54.76 356
## 45 Era 2: 2017-2020 Occasional 2 49.71 44.79 54.63 646
## 46 Era 2: 2017-2020 Light 2 45.35 38.60 52.09 423
## 47 Era 2: 2017-2020 Moderate 2 63.97 52.87 75.07 116
## 48 Era 2: 2017-2020 Heavy 2 54.58 43.47 65.69 117
## 49 Era 2: 2017-2020 Lifetime abstainer 3 7.40 4.10 10.69 29
## 50 Era 2: 2017-2020 Former drinker 3 10.20 7.36 13.03 93
## 51 Era 2: 2017-2020 Occasional 3 5.26 3.86 6.66 92
## 52 Era 2: 2017-2020 Light 3 2.43 1.39 3.47 35
## 53 Era 2: 2017-2020 Moderate 3 2.35 -0.19 4.89 8
## 54 Era 2: 2017-2020 Heavy 3 4.01 1.66 6.37 13
## 55 Era 2: 2017-2020 Lifetime abstainer 4 10.69 4.19 17.18 29
## 56 Era 2: 2017-2020 Former drinker 4 19.91 15.63 24.18 165
## 57 Era 2: 2017-2020 Occasional 4 8.46 6.25 10.67 134
## 58 Era 2: 2017-2020 Light 4 6.82 3.81 9.83 65
## 59 Era 2: 2017-2020 Moderate 4 8.65 -2.15 19.46 21
## 60 Era 2: 2017-2020 Heavy 4 8.37 3.79 12.96 27
## 61 Pooled 2011-2020 Lifetime abstainer 0 9.28 6.88 11.68 89
## 62 Pooled 2011-2020 Former drinker 0 4.16 2.68 5.64 55
## 63 Pooled 2011-2020 Occasional 0 9.35 7.41 11.29 277
## 64 Pooled 2011-2020 Light 0 13.21 11.22 15.20 228
## 65 Pooled 2011-2020 Moderate 0 7.60 4.52 10.67 40
## 66 Pooled 2011-2020 Heavy 0 8.99 5.19 12.79 31
## 67 Pooled 2011-2020 Lifetime abstainer 1 22.09 19.50 24.67 258
## 68 Pooled 2011-2020 Former drinker 1 15.86 13.54 18.18 245
## 69 Pooled 2011-2020 Occasional 1 26.71 24.34 29.08 838
## 70 Pooled 2011-2020 Light 1 32.23 29.22 35.25 625
## 71 Pooled 2011-2020 Moderate 1 23.96 19.23 28.69 137
## 72 Pooled 2011-2020 Heavy 1 17.76 13.03 22.49 84
## 73 Pooled 2011-2020 Lifetime abstainer 2 48.99 45.31 52.68 653
## 74 Pooled 2011-2020 Former drinker 2 49.71 46.50 52.92 867
## 75 Pooled 2011-2020 Occasional 2 50.28 47.81 52.76 1697
## 76 Pooled 2011-2020 Light 2 46.46 42.99 49.92 1103
## 77 Pooled 2011-2020 Moderate 2 56.77 51.48 62.07 348
## 78 Pooled 2011-2020 Heavy 2 62.00 56.31 67.69 315
## 79 Pooled 2011-2020 Lifetime abstainer 3 8.79 6.75 10.84 138
## 80 Pooled 2011-2020 Former drinker 3 10.12 8.57 11.68 255
## 81 Pooled 2011-2020 Occasional 3 5.11 4.25 5.97 239
## 82 Pooled 2011-2020 Light 3 2.64 1.83 3.45 100
## 83 Pooled 2011-2020 Moderate 3 2.67 1.48 3.86 32
## 84 Pooled 2011-2020 Heavy 3 3.33 1.53 5.13 23
## 85 Pooled 2011-2020 Lifetime abstainer 4 10.84 8.35 13.34 149
## 86 Pooled 2011-2020 Former drinker 4 20.15 17.58 22.72 387
## 87 Pooled 2011-2020 Occasional 4 8.55 7.34 9.76 348
## 88 Pooled 2011-2020 Light 4 5.46 4.05 6.87 146
## 89 Pooled 2011-2020 Moderate 4 9.00 4.75 13.25 63
## 90 Pooled 2011-2020 Heavy 4 7.92 5.01 10.83 53Table 3: adjusted odds ratios by category and by instrument era
Research question for this block: with the covariate set held fixed, how does each drinking level compare with lifetime abstention?
or_all <- rbind(or_e1, or_e2, or_po, or_hr, or_oc) # the five models of this arm, stacked
table3 <- or_all[startsWith(or_all$term, "drinkcat"), ] # keep the exposure terms only
table3$category <- sub("^drinkcat", "", table3$term) # drop the variable-name prefix
table3$outcome <- rep(outcname, nrow(table3)) # the outcome the cell-check rule chose
table3$outcome[startsWith(table3$model, "Check: high risk")] <- "high risk, stages 3 to 4" # except the binary check
table3$reference <- dplyr::if_else(startsWith(table3$model, "Sensitivity: occasional"), # which reference produced the row
"Occasional drinker", "Lifetime abstainer")
table3$or <- round(table3$estimate, 3) # the odds ratio
table3$lcl <- round(table3$conf.low, 3) # lower 95% confidence limit
table3$ucl <- round(table3$conf.high, 3) # upper 95% confidence limit
table3 <- table3[, c("model", "outcome", "nrows", "category", "reference", "or", "lcl", "ucl")] # the columns that get written
table3 # the table as it will be written## model outcome
## 1 Era 1: 2011-2016 ckmtier (1 low, 2 moderate, 3 high)
## 2 Era 1: 2011-2016 ckmtier (1 low, 2 moderate, 3 high)
## 3 Era 1: 2011-2016 ckmtier (1 low, 2 moderate, 3 high)
## 4 Era 1: 2011-2016 ckmtier (1 low, 2 moderate, 3 high)
## 5 Era 1: 2011-2016 ckmtier (1 low, 2 moderate, 3 high)
## 27 Era 2: 2017-2020 ckmtier (1 low, 2 moderate, 3 high)
## 28 Era 2: 2017-2020 ckmtier (1 low, 2 moderate, 3 high)
## 29 Era 2: 2017-2020 ckmtier (1 low, 2 moderate, 3 high)
## 30 Era 2: 2017-2020 ckmtier (1 low, 2 moderate, 3 high)
## 31 Era 2: 2017-2020 ckmtier (1 low, 2 moderate, 3 high)
## 53 Pooled 2011-2020 ckmtier (1 low, 2 moderate, 3 high)
## 54 Pooled 2011-2020 ckmtier (1 low, 2 moderate, 3 high)
## 55 Pooled 2011-2020 ckmtier (1 low, 2 moderate, 3 high)
## 56 Pooled 2011-2020 ckmtier (1 low, 2 moderate, 3 high)
## 57 Pooled 2011-2020 ckmtier (1 low, 2 moderate, 3 high)
## 80 Check: high risk, stages 3 to 4, pooled high risk, stages 3 to 4
## 81 Check: high risk, stages 3 to 4, pooled high risk, stages 3 to 4
## 82 Check: high risk, stages 3 to 4, pooled high risk, stages 3 to 4
## 83 Check: high risk, stages 3 to 4, pooled high risk, stages 3 to 4
## 84 Check: high risk, stages 3 to 4, pooled high risk, stages 3 to 4
## 104 Sensitivity: occasional drinkers as the reference, pooled ckmtier (1 low, 2 moderate, 3 high)
## 105 Sensitivity: occasional drinkers as the reference, pooled ckmtier (1 low, 2 moderate, 3 high)
## 106 Sensitivity: occasional drinkers as the reference, pooled ckmtier (1 low, 2 moderate, 3 high)
## 107 Sensitivity: occasional drinkers as the reference, pooled ckmtier (1 low, 2 moderate, 3 high)
## 108 Sensitivity: occasional drinkers as the reference, pooled ckmtier (1 low, 2 moderate, 3 high)
## nrows category reference or lcl ucl
## 1 5765 Former drinker Lifetime abstainer 1.257 1.003 1.576
## 2 5765 Occasional Lifetime abstainer 0.895 0.742 1.079
## 3 5765 Light Lifetime abstainer 0.712 0.570 0.890
## 4 5765 Moderate Lifetime abstainer 0.949 0.670 1.346
## 5 5765 Heavy Lifetime abstainer 1.015 0.751 1.373
## 27 3077 Former drinker Lifetime abstainer 0.851 0.503 1.439
## 28 3077 Occasional Lifetime abstainer 0.786 0.492 1.255
## 29 3077 Light Lifetime abstainer 0.665 0.420 1.052
## 30 3077 Moderate Lifetime abstainer 1.135 0.648 1.987
## 31 3077 Heavy Lifetime abstainer 0.567 0.294 1.093
## 53 8842 Former drinker Lifetime abstainer 1.091 0.878 1.356
## 54 8842 Occasional Lifetime abstainer 0.869 0.724 1.043
## 55 8842 Light Lifetime abstainer 0.707 0.583 0.856
## 56 8842 Moderate Lifetime abstainer 1.024 0.762 1.375
## 57 8842 Heavy Lifetime abstainer 0.818 0.611 1.095
## 80 8842 Former drinker Lifetime abstainer 0.929 0.694 1.244
## 81 8842 Occasional Lifetime abstainer 0.821 0.602 1.122
## 82 8842 Light Lifetime abstainer 0.537 0.353 0.819
## 83 8842 Moderate Lifetime abstainer 0.749 0.406 1.382
## 84 8842 Heavy Lifetime abstainer 0.432 0.254 0.735
## 104 8842 Lifetime abstainer Occasional drinker 1.150 0.958 1.380
## 105 8842 Former drinker Occasional drinker 1.256 1.053 1.499
## 106 8842 Light Occasional drinker 0.814 0.700 0.945
## 107 8842 Moderate Occasional drinker 1.179 0.908 1.530
## 108 8842 Heavy Occasional drinker 0.941 0.703 1.260Table 3 on the five-stage outcome: the sensitivity the rule made necessary
Research question for this block: would the arm have said something different if the thinnest cell had held 30 people and the outcome had stayed at five stages?
os_all <- rbind(os_e1, os_e2, os_po) # the three five-stage models
t3s <- os_all[startsWith(os_all$term, "drinkcat"), ] # keep the exposure terms only
t3s$category <- sub("^drinkcat", "", t3s$term) # drop the variable-name prefix
t3s$outcome <- rep("ckmstage (0 to 4)", nrow(t3s)) # the five-stage outcome, named in the table
t3s$reference <- rep("Lifetime abstainer", nrow(t3s)) # the reference group, stated in the table
t3s$or <- round(t3s$estimate, 3) # the odds ratio
t3s$lcl <- round(t3s$conf.low, 3) # lower 95% confidence limit
t3s$ucl <- round(t3s$conf.high, 3) # upper 95% confidence limit
t3s <- t3s[, c("model", "outcome", "nrows", "category", "reference", "or", "lcl", "ucl")] # the columns that get written
t3s # the table as it will be written## model outcome nrows category reference or lcl ucl
## 1 Era 1: 2011-2016 ckmstage (0 to 4) 5765 Former drinker Lifetime abstainer 1.236 0.975 1.567
## 2 Era 1: 2011-2016 ckmstage (0 to 4) 5765 Occasional Lifetime abstainer 0.936 0.787 1.114
## 3 Era 1: 2011-2016 ckmstage (0 to 4) 5765 Light Lifetime abstainer 0.702 0.567 0.870
## 4 Era 1: 2011-2016 ckmstage (0 to 4) 5765 Moderate Lifetime abstainer 0.960 0.706 1.305
## 5 Era 1: 2011-2016 ckmstage (0 to 4) 5765 Heavy Lifetime abstainer 0.943 0.686 1.297
## 29 Era 2: 2017-2020 ckmstage (0 to 4) 3077 Former drinker Lifetime abstainer 0.968 0.590 1.590
## 30 Era 2: 2017-2020 ckmstage (0 to 4) 3077 Occasional Lifetime abstainer 0.864 0.540 1.382
## 31 Era 2: 2017-2020 ckmstage (0 to 4) 3077 Light Lifetime abstainer 0.804 0.548 1.179
## 32 Era 2: 2017-2020 ckmstage (0 to 4) 3077 Moderate Lifetime abstainer 1.343 0.773 2.333
## 33 Era 2: 2017-2020 ckmstage (0 to 4) 3077 Heavy Lifetime abstainer 0.624 0.334 1.169
## 57 Pooled 2011-2020 ckmstage (0 to 4) 8842 Former drinker Lifetime abstainer 1.128 0.908 1.402
## 58 Pooled 2011-2020 ckmstage (0 to 4) 8842 Occasional Lifetime abstainer 0.918 0.767 1.099
## 59 Pooled 2011-2020 ckmstage (0 to 4) 8842 Light Lifetime abstainer 0.746 0.627 0.887
## 60 Pooled 2011-2020 ckmstage (0 to 4) 8842 Moderate Lifetime abstainer 1.083 0.824 1.423
## 61 Pooled 2011-2020 ckmstage (0 to 4) 8842 Heavy Lifetime abstainer 0.805 0.599 1.081write.csv(t3s, file.path(out_dir, "table-3-ckmstage.csv"), row.names = FALSE) # csv first, the file that is compared
flextable::save_as_docx(flextable::flextable(t3s), # then Word
path = file.path(out_dir, "table-3-ckmstage.docx"))Part 8: Figures
Figure 1: the same four categories under three reference groups
Research question for this block: how much of the alcohol-CKM gradient is the reference group rather than the alcohol?
or_fig <- rbind(or_po, or_a2, or_a3) # the three reference constructions, pooled rows
or_fig$category <- sub("^drinkcat|^nondrink|^drinkre", "", or_fig$term) # the category name without its variable prefix
fig1_dat <- or_fig[or_fig$category %in% c("Occasional", "Light", "Moderate", "Heavy"), ] # the four current-drinking categories
fig1_dat$category <- factor(fig1_dat$category, levels = c("Occasional", "Light", "Moderate", "Heavy")) # left to right by volume
fig1_dat$model <- factor(fig1_dat$model, levels = c("Pooled 2011-2020",
"Arm 2: pooled non-drinker reference", "Arm 3: former drinkers split on past heavy drinking"),
labels = c("Arm 1: lifetime abstainers", "Arm 2: pooled non-drinkers", "Arm 3: former drinkers split")) # readable legend
fig1_dat # the twelve points the figure draws## term estimate conf.low conf.high model nrows category
## 2 drinkcatOccasional 0.8690 0.7241 1.0429 Arm 1: lifetime abstainers 8842 Occasional
## 3 drinkcatLight 0.7067 0.5833 0.8560 Arm 1: lifetime abstainers 8842 Light
## 4 drinkcatModerate 1.0239 0.7625 1.3750 Arm 1: lifetime abstainers 8842 Moderate
## 5 drinkcatHeavy 0.8176 0.6106 1.0947 Arm 1: lifetime abstainers 8842 Heavy
## 27 nondrinkOccasional 0.8239 0.7131 0.9521 Arm 2: pooled non-drinkers 8842 Occasional
## 28 nondrinkLight 0.6696 0.5669 0.7909 Arm 2: pooled non-drinkers 8842 Light
## 29 nondrinkModerate 0.9682 0.7202 1.3015 Arm 2: pooled non-drinkers 8842 Moderate
## 30 nondrinkHeavy 0.7721 0.5661 1.0530 Arm 2: pooled non-drinkers 8842 Heavy
## 54 drinkreOccasional 0.8630 0.7191 1.0357 Arm 3: former drinkers split 8841 Occasional
## 55 drinkreLight 0.6988 0.5759 0.8479 Arm 3: former drinkers split 8841 Light
## 56 drinkreModerate 1.0092 0.7563 1.3466 Arm 3: former drinkers split 8841 Moderate
## 57 drinkreHeavy 0.8014 0.5987 1.0726 Arm 3: former drinkers split 8841 Heavyfig1 <- ggplot(fig1_dat, aes(x = category, y = estimate, ymin = conf.low, ymax = conf.high, colour = model)) +
geom_hline(yintercept = 1, linetype = "dashed", colour = "grey40") + # no association
geom_pointrange(position = position_dodge(width = 0.55), size = 0.5) + # estimate and 95% confidence interval
scale_y_log10(breaks = c(0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.5, 2.0)) + # odds ratios read on a log scale
scale_colour_manual(values = c("#003d99", "#990000", "#006633")) + # three colours that survive greyscale printing
labs(title = "Adjusted odds of a higher CKM tier, by drinking category and reference group",
subtitle = "NHANES 2011 to March 2020, survey weighted, adults 20 and older in the fasting subsample",
x = "Drinking category", y = "Adjusted proportional odds ratio (95% CI)", colour = "Reference group",
caption = "Adjusted for age, sex, race and ethnicity, education, income to poverty ratio, smoking, self-rated health, and physical activity.") +
theme_minimal(base_size = 12) + # a plain theme
theme(plot.title = element_text(size = 13, face = "bold"), legend.position = "bottom") # title and legend placement
fig1 # shown in the knitted page
ggsave(file.path(out_dir, "figure-1.png"), fig1, width = 8, height = 5, dpi = 150) # saved at the size the brief fixesPart 9: Save and knit
Research question for this block: can someone else tell exactly which file, which environment, and which decisions produced the output in this directory?
run_log <- c( # the short record that sits beside the tables
"EPI 601 Paper 1, arm 01: lifetime abstainers as the reference group",
paste("Run date:", format(Sys.time(), "%Y-%m-%d %H:%M:%S")), # when this knit happened
paste("R version:", R.version.string), # which R
paste("Input file:", data_file), # which file was read
paste("Input sha256:", tools::sha256sum(data_file)), # and its checksum, so the file can be identified
paste("Teaching sample flag (0 is the real file):", paste(sort(unique(dat$sampflag)), collapse = ", ")),
paste("Exclusions: rows as read", n_read), # the exclusion ladder, matching build/build-log.txt
paste("Exclusions: aged 20 years and older", n_age), # the exclusion ladder
paste("Exclusions: and not pregnant at examination", n_preg), # the exclusion ladder
paste("Exclusions: and in the fasting subsample", n_fast), # the exclusion ladder
paste("Exclusions: and CKM stage could be determined", n_samp), # the exclusion ladder
paste("Arm rows with a drinking category:", n_exp), # what the descriptive tables use
paste("Arm rows with every covariate observed:", n_mod), # what the adjusted models use
"Drinking categories from specs/00-cover.md: occasional 1 or fewer, light above 1 to 7, moderate above 7 to 14, heavy above 14 drinks a week",
paste("Cell check: thinnest exposure-by-stage cell =", min(tab)), # the cell check
paste("Cell check: outcome used =", outcname)) # and its verdict
writeLines(run_log, file.path(out_dir, "run-log.txt")) # written next to the tables and the figure
run_log # printed, so the knitted page carries it too## [1] "EPI 601 Paper 1, arm 01: lifetime abstainers as the reference group"
## [2] "Run date: 2026-09-14 23:11:01"
## [3] "R version: R version 4.5.1 (2025-06-13 ucrt)"
## [4] "Input file: C:/Users/safwa/OneDrive - University at Albany - SUNY/nhanes-ckm-reference-groups/data/nhanes_ckm_2011_2020.rds"
## [5] "Input sha256: NA"
## [6] "Teaching sample flag (0 is the real file): 0"
## [7] "Exclusions: rows as read 26280"
## [8] "Exclusions: aged 20 years and older 26280"
## [9] "Exclusions: and not pregnant at examination 26001"
## [10] "Exclusions: and in the fasting subsample 10679"
## [11] "Exclusions: and CKM stage could be determined 10585"
## [12] "Arm rows with a drinking category: 9823"
## [13] "Arm rows with every covariate observed: 8842"
## [14] "Drinking categories from specs/00-cover.md: occasional 1 or fewer, light above 1 to 7, moderate above 7 to 14, heavy above 14 drinks a week"
## [15] "Cell check: thinnest exposure-by-stage cell = 23"
## [16] "Cell check: outcome used = ckmtier (1 low, 2 moderate, 3 high)"## [1] "arm-01-lifetime-abstainers.html" "figure-1.png"
## [3] "run-log.txt" "table-1.csv"
## [5] "table-1.docx" "table-2.csv"
## [7] "table-2.docx" "table-3-ckmstage.csv"
## [9] "table-3-ckmstage.docx" "table-3.csv"
## [11] "table-3.docx"## R version 4.5.1 (2025-06-13 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_4.0.2 flextable_0.9.11 gtsummary_2.5.0 survey_4.5 survival_3.8-3
## [6] Matrix_1.7-3 dplyr_1.2.0
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.56 bslib_0.10.0 lattice_0.22-7
## [5] vctrs_0.7.1 tools_4.5.1 generics_0.1.4 tibble_3.3.1
## [9] pkgconfig_2.0.3 data.table_1.18.0 RColorBrewer_1.1-3 S7_0.2.1
## [13] gt_1.3.0 uuid_1.2-2 lifecycle_1.0.5 compiler_4.5.1
## [17] farver_2.1.2 stringr_1.6.0 textshaping_1.0.4 mitools_2.4
## [21] litedown_0.9 fontquiver_0.2.1 fontLiberation_0.1.0 htmltools_0.5.9
## [25] sass_0.4.10 yaml_2.3.12 pillar_1.11.1 jquerylib_0.1.4
## [29] tidyr_1.3.2 MASS_7.3-65 openssl_2.3.4 cachem_1.1.0
## [33] fontBitstreamVera_0.1.1 commonmark_2.0.0 tidyselect_1.2.1 zip_2.3.3
## [37] digest_0.6.39 stringi_1.8.7 purrr_1.2.1 splines_4.5.1
## [41] rprojroot_2.1.1 fastmap_1.2.0 here_1.0.2 cli_3.6.5
## [45] magrittr_2.0.4 patchwork_1.3.2 cards_0.7.1 broom_1.0.12
## [49] withr_3.0.2 backports_1.5.0 gdtools_0.5.0 scales_1.4.0
## [53] cardx_0.3.2 rmarkdown_2.30 officer_0.7.3 otel_0.2.0
## [57] askpass_1.2.1 ragg_1.5.0 evaluate_1.0.5 knitr_1.51
## [61] markdown_2.0 rlang_1.1.7 Rcpp_1.1.1 glue_1.8.0
## [65] DBI_1.2.3 xml2_1.5.2 rstudioapi_0.18.0 jsonlite_2.0.0
## [69] R6_2.6.1 fs_1.6.6 systemfonts_1.3.1Last knitted: September 14, 2026