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\\nhanes-ckm-reference-groups\\data\\sample\\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] 300 137##
## 1
## 300##
## 1 2
## 199 101Part 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 (118) clusters.
## svydesign(ids = ~SDMVPSU, strata = ~SDMVSTRA, weights = ~wtsafcmb,
## nest = TRUE, data = dat)## [1] 51Part 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] 300des_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] 300des_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] 300des_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] 300des_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] 300## [1] 51Part 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.0 12.0 54.2 52.0 365.0 54Drinks 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.0192 0.2308 1.0000 3.7081 3.5000 42.1154 110quantile(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.00 3.50 10.50 19.21 29.54The 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 33 4
## Former drinker 34 22
## Occasional 66 36
## Light 33 27
## Moderate 12 3
## Heavy 7 6
## <NA> 14 3The 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 67 26
## Occasional 66 36
## Light 33 27
## Moderate 12 3
## Heavy 7 6
## <NA> 14 3 table(des_sub$variables$frmheavy, des_sub$variables$era, useNA = "ifany") # arm 3's split, among former drinkers only##
## 1 2
## Former, not past heavy 28 19
## Former, past heavy 6 3
## <NA> 165 79 table(des_sub$variables$drinkre, des_sub$variables$era, useNA = "ifany") # arm 3's seven-category exposure by era##
## 1 2
## Lifetime abstainer 33 4
## Former, not past heavy 28 19
## Former, past heavy 6 3
## Occasional 66 36
## Light 33 27
## Moderate 12 3
## Heavy 7 6
## <NA> 14 3Sex
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
## 151 149Race 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
## 159 19 31 53
## Non-Hispanic Asian Other or multiracial
## 29 9Education
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
## 24 31 69 87 89Income 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.29 2.18 2.54 3.87 5.00 22Smoking
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
## 161 83 56Self-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>
## 45 85 100 54 15 1Physical 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 101 43
## Yes 98 58The 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 15 18 11 16
## Stage 1 15 12 11 22
## Stage 2 14 9 15 22
## Stage 3 15 16 10 19
## Stage 4 14 15 9 22Complete 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 24 15
## 1 175 86## [1] 261## [1] 39Part 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 8 6 9 8 6
## Former drinker 4 2 11 23 16
## Occasional 18 27 19 16 22
## Light 18 14 13 8 7
## Moderate 5 5 3 0 2
## Heavy 3 2 3 3 2## [1] 0##
## 0 1 2 3 4
## 60 60 60 60 60# 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 0 people, so the outcome used in every model below is ckmtier (1 low, 2 moderate, 3 high)##
## 1 2 3
## 120 60 120The 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 7 6 6 8 6
## Former drinker 2 1 6 14 11
## Occasional 14 16 14 10 12
## Light 11 8 6 6 2
## Moderate 4 4 2 0 2
## Heavy 2 0 2 2 1table(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 1 0 3 0 0
## Former drinker 2 1 5 9 5
## Occasional 4 11 5 6 10
## Light 7 6 7 2 5
## Moderate 1 1 1 0 0
## Heavy 1 2 1 1 1Table 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] 283##
## Lifetime abstainer Former drinker Occasional Light Moderate
## 37 56 102 60 15
## Heavy
## 13tbl1 <- 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 = 7,416,7161 |
Lifetime abstainer N = 540,8561 |
Former drinker N = 1,143,9121 |
Occasional N = 2,993,8871 |
Light N = 1,624,6841 |
Moderate N = 605,8811 |
Heavy N = 507,4941 |
|---|---|---|---|---|---|---|---|
| Age at screening, years | 48.5 (18.8) | 52.6 (20.9) | 61.3 (15.4) | 49.0 (18.4) | 42.7 (18.2) | 43.8 (15.8) | 35.8 (15.5) |
| Sex | |||||||
| Male | 3,414,786 (46%) | 124,054 (23%) | 360,733 (32%) | 1,351,381 (45%) | 851,919 (52%) | 253,459 (42%) | 473,240 (93%) |
| Female | 4,001,930 (54%) | 416,802 (77%) | 783,179 (68%) | 1,642,506 (55%) | 772,765 (48%) | 352,423 (58%) | 34,254 (6.7%) |
| Race and ethnicity | |||||||
| Non-Hispanic White | 5,928,015 (80%) | 318,927 (59%) | 907,402 (79%) | 2,377,010 (79%) | 1,321,985 (81%) | 552,325 (91%) | 450,365 (89%) |
| Mexican American | 263,902 (3.6%) | 64,085 (12%) | 49,905 (4.4%) | 52,815 (1.8%) | 65,784 (4.0%) | 13,006 (2.1%) | 18,308 (3.6%) |
| Other Hispanic | 393,972 (5.3%) | 34,425 (6.4%) | 48,138 (4.2%) | 275,329 (9.2%) | 36,080 (2.2%) | 0 (0%) | 0 (0%) |
| Non-Hispanic Black | 460,532 (6.2%) | 49,811 (9.2%) | 99,857 (8.7%) | 127,294 (4.3%) | 104,199 (6.4%) | 40,550 (6.7%) | 38,821 (7.6%) |
| Non-Hispanic Asian | 224,224 (3.0%) | 56,071 (10%) | 38,610 (3.4%) | 65,328 (2.2%) | 64,214 (4.0%) | 0 (0%) | 0 (0%) |
| Other or multiracial | 146,071 (2.0%) | 17,537 (3.2%) | 0 (0%) | 96,112 (3.2%) | 32,422 (2.0%) | 0 (0%) | 0 (0%) |
| Education | |||||||
| Less than 9th grade | 223,184 (3.0%) | 18,618 (3.4%) | 132,708 (12%) | 42,886 (1.4%) | 15,968 (1.0%) | 13,006 (2.1%) | 0 (0%) |
| 9th to 11th grade | 397,591 (5.4%) | 36,048 (6.7%) | 79,156 (6.9%) | 193,298 (6.5%) | 34,179 (2.1%) | 5,396 (0.9%) | 49,515 (9.8%) |
| High school or GED | 1,949,499 (26%) | 114,873 (21%) | 407,191 (36%) | 1,166,089 (39%) | 178,741 (11%) | 38,993 (6.4%) | 43,612 (8.6%) |
| Some college | 1,975,808 (27%) | 193,982 (36%) | 168,283 (15%) | 669,169 (22%) | 508,610 (31%) | 218,987 (36%) | 216,777 (43%) |
| College graduate | 2,870,634 (39%) | 177,336 (33%) | 356,575 (31%) | 922,446 (31%) | 887,187 (55%) | 329,499 (54%) | 197,590 (39%) |
| Family income to poverty ratio | 3.3 (1.6) | 2.4 (1.6) | 3.0 (1.5) | 3.2 (1.5) | 3.7 (1.4) | 3.2 (1.8) | 3.5 (1.8) |
| Missing | 464,724 | 48,062 | 86,405 | 135,775 | 194,482 | 0 | 0 |
| Smoking status | |||||||
| Never | 3,732,050 (50%) | 475,115 (88%) | 472,720 (41%) | 1,393,361 (47%) | 962,015 (59%) | 373,935 (62%) | 54,905 (11%) |
| Former | 2,305,528 (31%) | 50,607 (9.4%) | 335,174 (29%) | 1,118,522 (37%) | 419,956 (26%) | 208,919 (34%) | 172,349 (34%) |
| Current | 1,379,138 (19%) | 15,134 (2.8%) | 336,019 (29%) | 482,004 (16%) | 242,714 (15%) | 23,027 (3.8%) | 280,240 (55%) |
| Self-rated health (HUQ010) | |||||||
| Excellent | 1,278,730 (17%) | 110,996 (21%) | 135,300 (12%) | 381,989 (13%) | 515,600 (32%) | 100,590 (17%) | 34,254 (6.7%) |
| Very good | 2,863,859 (39%) | 182,247 (34%) | 310,973 (27%) | 1,174,055 (39%) | 566,322 (35%) | 386,501 (64%) | 243,761 (48%) |
| Good | 2,255,694 (30%) | 129,604 (24%) | 321,723 (28%) | 1,018,616 (34%) | 503,767 (31%) | 102,825 (17%) | 179,158 (35%) |
| Fair | 746,729 (10%) | 75,389 (14%) | 269,407 (24%) | 296,651 (9.9%) | 38,996 (2.4%) | 15,966 (2.6%) | 50,320 (9.9%) |
| Poor | 263,624 (3.6%) | 42,620 (7.9%) | 98,427 (8.7%) | 122,578 (4.1%) | 0 (0%) | 0 (0%) | 0 (0%) |
| Missing | 8,081 | 0 | 8,081 | 0 | 0 | 0 | 0 |
| Any vigorous or moderate physical activity | 4,694,855 (63%) | 198,910 (37%) | 581,010 (51%) | 1,889,195 (63%) | 1,018,660 (63%) | 522,633 (86%) | 484,447 (95%) |
| Alcohol instrument era | |||||||
| 2011-2016 (ALQ101, ALQ110) | 4,336,242 (58%) | 491,867 (91%) | 646,957 (57%) | 1,844,697 (62%) | 777,350 (48%) | 398,540 (66%) | 176,831 (35%) |
| 2017-2020 (ALQ111) | 3,080,474 (42%) | 48,990 (9.1%) | 496,955 (43%) | 1,149,190 (38%) | 847,335 (52%) | 207,341 (34%) | 330,663 (65%) |
| CKM syndrome stage | |||||||
| Stage 0 | 2,081,926 (28%) | 202,848 (38%) | 121,997 (11%) | 533,623 (18%) | 715,252 (44%) | 290,535 (48%) | 217,671 (43%) |
| Stage 1 | 1,720,412 (23%) | 83,221 (15%) | 31,710 (2.8%) | 1,024,721 (34%) | 300,119 (18%) | 119,577 (20%) | 161,064 (32%) |
| Stage 2 | 1,534,726 (21%) | 70,538 (13%) | 324,909 (28%) | 627,962 (21%) | 311,596 (19%) | 126,429 (21%) | 73,292 (14%) |
| Stage 3 | 900,136 (12%) | 112,149 (21%) | 314,511 (27%) | 306,063 (10%) | 123,307 (7.6%) | 0 (0%) | 44,106 (8.7%) |
| Stage 4 | 1,179,516 (16%) | 72,100 (13%) | 350,786 (31%) | 501,518 (17%) | 174,411 (11%) | 69,340 (11%) | 11,361 (2.2%) |
| 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] 283## [1] 17des_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] 261des_e1 <- subset(des_mod, era == 1) # 2011 to 2016, ALQ101 and ALQ110
nrow(des_e1) # era 1 model rows## [1] 175## [1] 86##
## 1 2
## Lifetime abstainer 30 3
## Former drinker 31 18
## Occasional 64 32
## Light 31 24
## Moderate 12 3
## Heavy 7 6table(des_mod$variables$drinkcat, des_mod$variables$outcome) # model rows by category and the outcome the rule chose##
## 1 2 3
## Lifetime abstainer 14 8 11
## Former drinker 6 9 34
## Occasional 43 18 35
## Light 30 12 13
## Moderate 10 3 2
## Heavy 5 3 5The 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.3751 0.15387 0.1304 0.20735
## Former drinker Former drinker 0.1066 0.02772 0.2840 0.27494
## Occasional Occasional 0.1782 0.34227 0.2097 0.10223
## Light Light 0.4402 0.18472 0.1918 0.07590
## Moderate Moderate 0.4795 0.19736 0.2087 0.00000
## Heavy Heavy 0.4289 0.31737 0.1444 0.08691
## stagefStage 4 se.stagefStage 0 se.stagefStage 1 se.stagefStage 2
## Lifetime abstainer 0.13331 0.12836 0.06764 0.05183
## Former drinker 0.30665 0.05422 0.02161 0.09975
## Occasional 0.16751 0.04760 0.08333 0.05568
## Light 0.10735 0.08811 0.05170 0.06221
## Moderate 0.11445 0.15994 0.10203 0.13440
## Heavy 0.02239 0.19797 0.20382 0.09727
## se.stagefStage 3 se.stagefStage 4
## Lifetime abstainer 0.07078 0.06841
## Former drinker 0.07142 0.09255
## Occasional 0.03141 0.06334
## Light 0.02968 0.05619
## Moderate 0.00000 0.07983
## Heavy 0.06549 0.01839ci_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 37.51 15.39 13.04 20.74 13.33
## Former drinker 10.66 2.77 28.40 27.49 30.67
## Occasional 17.82 34.23 20.97 10.22 16.75
## Light 44.02 18.47 19.18 7.59 10.74
## Moderate 47.95 19.74 20.87 0.00 11.44
## Heavy 42.89 31.74 14.44 8.69 2.24The 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] 185sb_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.3527 0.16919 0.1035 0.22801
## Former drinker Former drinker 0.1031 0.01669 0.2797 0.28433
## Occasional Occasional 0.2068 0.30603 0.2595 0.10751
## Light Light 0.5018 0.14576 0.2375 0.09552
## Moderate Moderate 0.4564 0.26741 0.1022 0.00000
## Heavy Heavy 0.4672 0.00000 0.3109 0.18973
## stagefStage 4 se.stagefStage 0 se.stagefStage 1 se.stagefStage 2
## Lifetime abstainer 0.14659 0.13963 0.07504 0.04860
## Former drinker 0.31625 0.07324 0.01697 0.12793
## Occasional 0.12019 0.06162 0.07994 0.07946
## Light 0.01939 0.11555 0.04463 0.09877
## Moderate 0.17399 0.16933 0.13771 0.08696
## Heavy 0.03211 0.22866 0.00000 0.19724
## se.stagefStage 3 se.stagefStage 4
## Lifetime abstainer 0.07813 0.07543
## Former drinker 0.08330 0.09084
## Occasional 0.03790 0.04869
## Light 0.03290 0.01428
## Moderate 0.00000 0.11505
## Heavy 0.14960 0.03425The 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] 98sb_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.5995 0.00000 0.40055 0.00000
## Former drinker Former drinker 0.1113 0.04208 0.28971 0.26273
## Occasional Occasional 0.1324 0.40044 0.12993 0.09375
## Light Light 0.3837 0.22047 0.14986 0.05789
## Moderate Moderate 0.5240 0.06273 0.41327 0.00000
## Heavy Heavy 0.4084 0.48709 0.05537 0.03192
## stagefStage 4 se.stagefStage 0 se.stagefStage 1 se.stagefStage 2
## Lifetime abstainer 0.00000 0.27884 0.00000 0.27884
## Former drinker 0.29416 0.08078 0.04582 0.15841
## Occasional 0.24348 0.07285 0.16603 0.06387
## Light 0.18804 0.13229 0.09090 0.07018
## Moderate 0.00000 0.33195 0.07217 0.32613
## Heavy 0.01719 0.27841 0.27629 0.06303
## se.stagefStage 3 se.stagefStage 4
## Lifetime abstainer 0.00000 0.00000
## Former drinker 0.12295 0.17889
## Occasional 0.05410 0.15307
## Light 0.04929 0.09908
## Moderate 0.00000 0.00000
## Heavy 0.03699 0.02013The 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.373e+01 2.237e+00 8.422e+01 Era 1: 2011-2016 175
## 2 drinkcatOccasional 1.914e+00 4.014e-01 9.131e+00 Era 1: 2011-2016 175
## 3 drinkcatLight 2.835e+00 3.690e-01 2.178e+01 Era 1: 2011-2016 175
## 4 drinkcatModerate 7.167e-01 3.981e-02 1.290e+01 Era 1: 2011-2016 175
## 5 drinkcatHeavy 2.185e+01 1.448e+00 3.298e+02 Era 1: 2011-2016 175
## 6 RIDAGEYR 1.104e+00 1.078e+00 1.131e+00 Era 1: 2011-2016 175
## 7 sexFemale 3.851e-01 1.036e-01 1.431e+00 Era 1: 2011-2016 175
## 8 racecatMexican American 4.089e+00 8.242e-02 2.028e+02 Era 1: 2011-2016 175
## 9 racecatOther Hispanic 5.225e-01 8.330e-02 3.277e+00 Era 1: 2011-2016 175
## 10 racecatNon-Hispanic Black 6.620e-01 1.598e-01 2.742e+00 Era 1: 2011-2016 175
## 11 racecatNon-Hispanic Asian 7.759e-01 1.158e-01 5.198e+00 Era 1: 2011-2016 175
## 12 racecatOther or multiracial 8.009e-02 3.486e-03 1.840e+00 Era 1: 2011-2016 175
## 13 educat9th to 11th grade 1.562e+01 3.961e-01 6.161e+02 Era 1: 2011-2016 175
## 14 educatHigh school or GED 3.630e+00 1.974e-01 6.674e+01 Era 1: 2011-2016 175
## 15 educatSome college 5.955e+00 2.800e-01 1.267e+02 Era 1: 2011-2016 175
## 16 educatCollege graduate 5.730e+00 2.707e-01 1.213e+02 Era 1: 2011-2016 175
## 17 pir 7.596e-01 5.108e-01 1.130e+00 Era 1: 2011-2016 175
## 18 smokecatFormer 2.234e+00 7.879e-01 6.336e+00 Era 1: 2011-2016 175
## 19 smokecatCurrent 3.949e-01 1.130e-01 1.380e+00 Era 1: 2011-2016 175
## 20 srhVery good 1.049e+00 1.129e-01 9.734e+00 Era 1: 2011-2016 175
## 21 srhGood 3.466e+00 3.942e-01 3.047e+01 Era 1: 2011-2016 175
## 22 srhFair 2.501e+01 1.536e+00 4.071e+02 Era 1: 2011-2016 175
## 23 srhPoor 1.391e+00 3.129e-02 6.178e+01 Era 1: 2011-2016 175
## 24 physactYes 9.485e-01 3.123e-01 2.881e+00 Era 1: 2011-2016 175
## 25 1|2 6.900e+02 1.650e+01 2.885e+04 Era 1: 2011-2016 175
## 26 2|3 6.012e+03 1.202e+02 3.006e+05 Era 1: 2011-2016 175The 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## Era 2 model did not fit on this file: attempt to find suitable starting values failedor_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## [1] term estimate conf.low conf.high model nrows
## <0 rows> (or 0-length row.names)The 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 5.709e+00 7.349e-01 4.434e+01 Pooled 2011-2020 261
## 2 drinkcatOccasional 1.705e+00 2.709e-01 1.073e+01 Pooled 2011-2020 261
## 3 drinkcatLight 1.838e+00 2.168e-01 1.559e+01 Pooled 2011-2020 261
## 4 drinkcatModerate 1.583e+00 1.503e-01 1.667e+01 Pooled 2011-2020 261
## 5 drinkcatHeavy 5.747e-01 1.723e-02 1.917e+01 Pooled 2011-2020 261
## 6 RIDAGEYR 1.111e+00 1.083e+00 1.139e+00 Pooled 2011-2020 261
## 7 sexFemale 3.989e-01 1.344e-01 1.183e+00 Pooled 2011-2020 261
## 8 racecatMexican American 9.145e-01 1.610e-01 5.196e+00 Pooled 2011-2020 261
## 9 racecatOther Hispanic 3.545e-01 8.564e-02 1.468e+00 Pooled 2011-2020 261
## 10 racecatNon-Hispanic Black 5.150e-01 1.843e-01 1.439e+00 Pooled 2011-2020 261
## 11 racecatNon-Hispanic Asian 9.649e-01 1.956e-01 4.760e+00 Pooled 2011-2020 261
## 12 racecatOther or multiracial 9.410e-02 5.853e-03 1.513e+00 Pooled 2011-2020 261
## 13 educat9th to 11th grade 5.729e+00 7.574e-01 4.333e+01 Pooled 2011-2020 261
## 14 educatHigh school or GED 5.551e+00 1.479e+00 2.084e+01 Pooled 2011-2020 261
## 15 educatSome college 7.555e+00 1.596e+00 3.576e+01 Pooled 2011-2020 261
## 16 educatCollege graduate 3.521e+00 7.061e-01 1.756e+01 Pooled 2011-2020 261
## 17 pir 9.728e-01 6.819e-01 1.388e+00 Pooled 2011-2020 261
## 18 smokecatFormer 1.240e+00 4.475e-01 3.437e+00 Pooled 2011-2020 261
## 19 smokecatCurrent 9.218e-01 2.653e-01 3.204e+00 Pooled 2011-2020 261
## 20 srhVery good 3.206e+00 4.221e-01 2.435e+01 Pooled 2011-2020 261
## 21 srhGood 6.354e+00 8.011e-01 5.039e+01 Pooled 2011-2020 261
## 22 srhFair 4.641e+01 5.302e+00 4.063e+02 Pooled 2011-2020 261
## 23 srhPoor 1.806e+00 7.398e-02 4.408e+01 Pooled 2011-2020 261
## 24 physactYes 4.919e-01 2.013e-01 1.202e+00 Pooled 2011-2020 261
## 25 1|2 1.803e+03 5.447e+01 5.965e+04 Pooled 2011-2020 261
## 26 2|3 1.518e+04 3.380e+02 6.820e+05 Pooled 2011-2020 261## 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 1.74197 1.04594 1.66546
## drinkcatOccasional 0.53343 0.93854 0.56837
## drinkcatLight 0.60893 1.09075 0.55827
## drinkcatModerate 0.45930 1.20116 0.38238
## drinkcatHeavy -0.55385 1.78939 -0.30952
## RIDAGEYR 0.10500 0.01279 8.20770
## sexFemale -0.91915 0.55483 -1.65662
## racecatMexican American -0.08938 0.88638 -0.10084
## racecatOther Hispanic -1.03692 0.72486 -1.43051
## racecatNon-Hispanic Black -0.66358 0.52421 -1.26585
## racecatNon-Hispanic Asian -0.03576 0.81432 -0.04392
## racecatOther or multiracial -2.36343 1.41709 -1.66780
## educat9th to 11th grade 1.74547 1.03231 1.69084
## educatHigh school or GED 1.71406 0.67489 2.53977
## educatSome college 2.02223 0.79323 2.54937
## educatCollege graduate 1.25880 0.81983 1.53545
## pir -0.02760 0.18123 -0.15231
## smokecatFormer 0.21536 0.52008 0.41410
## smokecatCurrent -0.08139 0.63555 -0.12807
## srhVery good 1.16506 1.03444 1.12628
## srhGood 1.84903 1.05654 1.75008
## srhFair 3.83758 1.10688 3.46704
## srhPoor 0.59107 1.63014 0.36259
## physactYes -0.70957 0.45581 -1.55672
##
## Intercepts:
## Value Std. Error t value
## 1|2 7.497 1.785 4.199
## 2|3 9.628 1.941 4.959Sensitivity: 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 4.8652 0.77225 3.065e+01 Era 1: 2011-2016 175
## 2 drinkcatOccasional 1.7736 0.34489 9.121e+00 Era 1: 2011-2016 175
## 3 drinkcatLight 1.2864 0.15790 1.048e+01 Era 1: 2011-2016 175
## 4 drinkcatModerate 0.6953 0.06339 7.627e+00 Era 1: 2011-2016 175
## 5 drinkcatHeavy 2.1054 0.11029 4.019e+01 Era 1: 2011-2016 175
## 6 RIDAGEYR 1.0834 1.05793 1.110e+00 Era 1: 2011-2016 175
## 7 sexFemale 0.5595 0.21709 1.442e+00 Era 1: 2011-2016 175
## 8 racecatMexican American 0.7317 0.01940 2.760e+01 Era 1: 2011-2016 175
## 9 racecatOther Hispanic 0.7353 0.24087 2.245e+00 Era 1: 2011-2016 175
## 10 racecatNon-Hispanic Black 1.3997 0.55108 3.555e+00 Era 1: 2011-2016 175
## 11 racecatNon-Hispanic Asian 1.0482 0.26104 4.209e+00 Era 1: 2011-2016 175
## 12 racecatOther or multiracial 0.1031 0.01434 7.414e-01 Era 1: 2011-2016 175
## 13 educat9th to 11th grade 1.3589 0.07552 2.445e+01 Era 1: 2011-2016 175
## 14 educatHigh school or GED 1.4125 0.17748 1.124e+01 Era 1: 2011-2016 175
## 15 educatSome college 1.4222 0.19779 1.023e+01 Era 1: 2011-2016 175
## 16 educatCollege graduate 1.4327 0.17646 1.163e+01 Era 1: 2011-2016 175
## 17 pir 0.7248 0.52206 1.006e+00 Era 1: 2011-2016 175
## 18 smokecatFormer 2.2655 0.99992 5.133e+00 Era 1: 2011-2016 175
## 19 smokecatCurrent 0.8803 0.22651 3.421e+00 Era 1: 2011-2016 175
## 20 srhVery good 1.3335 0.23821 7.465e+00 Era 1: 2011-2016 175
## 21 srhGood 4.4866 0.92216 2.183e+01 Era 1: 2011-2016 175
## 22 srhFair 16.7951 2.37944 1.185e+02 Era 1: 2011-2016 175
## 23 srhPoor 5.5219 0.22037 1.384e+02 Era 1: 2011-2016 175
## 24 physactYes 1.5565 0.53593 4.521e+00 Era 1: 2011-2016 175
## 25 0|1 28.3581 1.60563 5.009e+02 Era 1: 2011-2016 175
## 26 1|2 116.0230 5.76550 2.335e+03 Era 1: 2011-2016 175
## 27 2|3 864.4893 33.36154 2.240e+04 Era 1: 2011-2016 175
## 28 3|4 3998.8979 150.30090 1.064e+05 Era 1: 2011-2016 175Sensitivity 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## Era 2 five-stage model did not fit on this file: attempt to find suitable starting values failedos_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## [1] term estimate conf.low conf.high model nrows
## <0 rows> (or 0-length row.names)Sensitivity 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 3.400e+00 0.44512 2.596e+01 Pooled 2011-2020 261
## 2 drinkcatOccasional 1.979e+00 0.28516 1.374e+01 Pooled 2011-2020 261
## 3 drinkcatLight 1.532e+00 0.16277 1.442e+01 Pooled 2011-2020 261
## 4 drinkcatModerate 1.148e+00 0.09278 1.421e+01 Pooled 2011-2020 261
## 5 drinkcatHeavy 6.598e-01 0.02680 1.624e+01 Pooled 2011-2020 261
## 6 RIDAGEYR 1.084e+00 1.06141 1.108e+00 Pooled 2011-2020 261
## 7 sexFemale 4.479e-01 0.21833 9.187e-01 Pooled 2011-2020 261
## 8 racecatMexican American 1.040e+00 0.35651 3.031e+00 Pooled 2011-2020 261
## 9 racecatOther Hispanic 7.043e-01 0.23923 2.074e+00 Pooled 2011-2020 261
## 10 racecatNon-Hispanic Black 1.246e+00 0.58397 2.660e+00 Pooled 2011-2020 261
## 11 racecatNon-Hispanic Asian 1.028e+00 0.33829 3.122e+00 Pooled 2011-2020 261
## 12 racecatOther or multiracial 1.371e-01 0.03209 5.859e-01 Pooled 2011-2020 261
## 13 educat9th to 11th grade 1.080e+00 0.22582 5.169e+00 Pooled 2011-2020 261
## 14 educatHigh school or GED 2.725e+00 0.59486 1.248e+01 Pooled 2011-2020 261
## 15 educatSome college 2.209e+00 0.49033 9.949e+00 Pooled 2011-2020 261
## 16 educatCollege graduate 1.598e+00 0.41383 6.168e+00 Pooled 2011-2020 261
## 17 pir 1.012e+00 0.78087 1.312e+00 Pooled 2011-2020 261
## 18 smokecatFormer 1.325e+00 0.66148 2.654e+00 Pooled 2011-2020 261
## 19 smokecatCurrent 2.058e+00 0.58521 7.238e+00 Pooled 2011-2020 261
## 20 srhVery good 4.178e+00 1.02703 1.700e+01 Pooled 2011-2020 261
## 21 srhGood 7.247e+00 2.04822 2.564e+01 Pooled 2011-2020 261
## 22 srhFair 3.318e+01 6.48422 1.698e+02 Pooled 2011-2020 261
## 23 srhPoor 5.763e+00 0.48108 6.903e+01 Pooled 2011-2020 261
## 24 physactYes 6.333e-01 0.27827 1.441e+00 Pooled 2011-2020 261
## 25 0|1 8.503e+01 5.44018 1.329e+03 Pooled 2011-2020 261
## 26 1|2 4.698e+02 27.37435 8.063e+03 Pooled 2011-2020 261
## 27 2|3 3.162e+03 143.69857 6.957e+04 Pooled 2011-2020 261
## 28 3|4 1.172e+04 477.97858 2.872e+05 Pooled 2011-2020 261Check: 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
## 161 100fit_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) 3.212e-07 7.344e-10 1.405e-04
## 2 drinkcatFormer drinker 2.178e+00 3.623e-01 1.309e+01
## 3 drinkcatOccasional 1.065e+00 2.106e-01 5.384e+00
## 4 drinkcatLight 5.174e-01 4.867e-02 5.501e+00
## 5 drinkcatModerate 2.148e+00 4.107e-02 1.124e+02
## 6 drinkcatHeavy 2.675e-01 1.732e-02 4.131e+00
## 7 RIDAGEYR 1.213e+00 1.114e+00 1.320e+00
## 8 sexFemale 1.091e-01 3.082e-02 3.859e-01
## 9 racecatMexican American 2.270e+00 1.888e-01 2.729e+01
## 10 racecatOther Hispanic 8.964e-01 6.701e-02 1.199e+01
## 11 racecatNon-Hispanic Black 1.549e+00 1.915e-01 1.253e+01
## 12 racecatNon-Hispanic Asian 1.834e+00 1.585e-01 2.123e+01
## 13 racecatOther or multiracial 5.652e-08 3.778e-09 8.456e-07
## 14 educat9th to 11th grade 3.584e+00 1.402e-01 9.160e+01
## 15 educatHigh school or GED 4.436e+01 3.931e+00 5.005e+02
## 16 educatSome college 3.955e+01 1.878e+00 8.326e+02
## 17 educatCollege graduate 6.593e+00 5.371e-01 8.093e+01
## 18 pir 9.152e-01 5.757e-01 1.455e+00
## 19 smokecatFormer 1.140e+00 2.836e-01 4.580e+00
## 20 smokecatCurrent 4.496e+00 8.784e-01 2.301e+01
## 21 srhVery good 1.411e+00 1.248e-01 1.595e+01
## 22 srhGood 5.505e+00 3.606e-01 8.403e+01
## 23 srhFair 7.233e+01 6.597e+00 7.930e+02
## 24 srhPoor 9.798e+00 6.082e-01 1.579e+02
## 25 physactYes 3.519e-01 7.393e-02 1.675e+00
## model nrows
## 1 Check: high risk, stages 3 to 4, pooled 261
## 2 Check: high risk, stages 3 to 4, pooled 261
## 3 Check: high risk, stages 3 to 4, pooled 261
## 4 Check: high risk, stages 3 to 4, pooled 261
## 5 Check: high risk, stages 3 to 4, pooled 261
## 6 Check: high risk, stages 3 to 4, pooled 261
## 7 Check: high risk, stages 3 to 4, pooled 261
## 8 Check: high risk, stages 3 to 4, pooled 261
## 9 Check: high risk, stages 3 to 4, pooled 261
## 10 Check: high risk, stages 3 to 4, pooled 261
## 11 Check: high risk, stages 3 to 4, pooled 261
## 12 Check: high risk, stages 3 to 4, pooled 261
## 13 Check: high risk, stages 3 to 4, pooled 261
## 14 Check: high risk, stages 3 to 4, pooled 261
## 15 Check: high risk, stages 3 to 4, pooled 261
## 16 Check: high risk, stages 3 to 4, pooled 261
## 17 Check: high risk, stages 3 to 4, pooled 261
## 18 Check: high risk, stages 3 to 4, pooled 261
## 19 Check: high risk, stages 3 to 4, pooled 261
## 20 Check: high risk, stages 3 to 4, pooled 261
## 21 Check: high risk, stages 3 to 4, pooled 261
## 22 Check: high risk, stages 3 to 4, pooled 261
## 23 Check: high risk, stages 3 to 4, pooled 261
## 24 Check: high risk, stages 3 to 4, pooled 261
## 25 Check: high risk, stages 3 to 4, pooled 261The 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
## 96 33 49 55 15
## Heavy
## 13fit_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 5.866e-01 9.321e-02 3.692e+00
## 2 drinkcatFormer drinker 3.349e+00 1.011e+00 1.109e+01
## 3 drinkcatLight 1.078e+00 3.291e-01 3.534e+00
## 4 drinkcatModerate 9.285e-01 2.799e-01 3.080e+00
## 5 drinkcatHeavy 3.371e-01 1.584e-02 7.175e+00
## 6 RIDAGEYR 1.111e+00 1.083e+00 1.139e+00
## 7 sexFemale 3.989e-01 1.344e-01 1.183e+00
## 8 racecatMexican American 9.145e-01 1.609e-01 5.196e+00
## 9 racecatOther Hispanic 3.545e-01 8.564e-02 1.468e+00
## 10 racecatNon-Hispanic Black 5.150e-01 1.843e-01 1.439e+00
## 11 racecatNon-Hispanic Asian 9.648e-01 1.956e-01 4.760e+00
## 12 racecatOther or multiracial 9.409e-02 5.852e-03 1.513e+00
## 13 educat9th to 11th grade 5.729e+00 7.575e-01 4.333e+01
## 14 educatHigh school or GED 5.552e+00 1.479e+00 2.084e+01
## 15 educatSome college 7.556e+00 1.596e+00 3.577e+01
## 16 educatCollege graduate 3.521e+00 7.061e-01 1.756e+01
## 17 pir 9.728e-01 6.820e-01 1.388e+00
## 18 smokecatFormer 1.240e+00 4.475e-01 3.437e+00
## 19 smokecatCurrent 9.218e-01 2.653e-01 3.204e+00
## 20 srhVery good 3.206e+00 4.222e-01 2.435e+01
## 21 srhGood 6.354e+00 8.011e-01 5.039e+01
## 22 srhFair 4.642e+01 5.302e+00 4.063e+02
## 23 srhPoor 1.806e+00 7.398e-02 4.408e+01
## 24 physactYes 4.919e-01 2.013e-01 1.202e+00
## 25 1|2 1.057e+03 4.475e+01 2.499e+04
## 26 2|3 8.906e+03 2.616e+02 3.031e+05
## model nrows
## 1 Sensitivity: occasional drinkers as the reference, pooled 261
## 2 Sensitivity: occasional drinkers as the reference, pooled 261
## 3 Sensitivity: occasional drinkers as the reference, pooled 261
## 4 Sensitivity: occasional drinkers as the reference, pooled 261
## 5 Sensitivity: occasional drinkers as the reference, pooled 261
## 6 Sensitivity: occasional drinkers as the reference, pooled 261
## 7 Sensitivity: occasional drinkers as the reference, pooled 261
## 8 Sensitivity: occasional drinkers as the reference, pooled 261
## 9 Sensitivity: occasional drinkers as the reference, pooled 261
## 10 Sensitivity: occasional drinkers as the reference, pooled 261
## 11 Sensitivity: occasional drinkers as the reference, pooled 261
## 12 Sensitivity: occasional drinkers as the reference, pooled 261
## 13 Sensitivity: occasional drinkers as the reference, pooled 261
## 14 Sensitivity: occasional drinkers as the reference, pooled 261
## 15 Sensitivity: occasional drinkers as the reference, pooled 261
## 16 Sensitivity: occasional drinkers as the reference, pooled 261
## 17 Sensitivity: occasional drinkers as the reference, pooled 261
## 18 Sensitivity: occasional drinkers as the reference, pooled 261
## 19 Sensitivity: occasional drinkers as the reference, pooled 261
## 20 Sensitivity: occasional drinkers as the reference, pooled 261
## 21 Sensitivity: occasional drinkers as the reference, pooled 261
## 22 Sensitivity: occasional drinkers as the reference, pooled 261
## 23 Sensitivity: occasional drinkers as the reference, pooled 261
## 24 Sensitivity: occasional drinkers as the reference, pooled 261
## 25 Sensitivity: occasional drinkers as the reference, pooled 261
## 26 Sensitivity: occasional drinkers as the reference, pooled 261The 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] 261 fit_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 4.667e-01 1.467e-01 1.485e+00 Arm 2: pooled non-drinker reference
## 2 nondrinkLight 5.166e-01 1.223e-01 2.182e+00 Arm 2: pooled non-drinker reference
## 3 nondrinkModerate 4.476e-01 8.308e-02 2.411e+00 Arm 2: pooled non-drinker reference
## 4 nondrinkHeavy 1.602e-01 7.719e-03 3.323e+00 Arm 2: pooled non-drinker reference
## 5 RIDAGEYR 1.110e+00 1.083e+00 1.138e+00 Arm 2: pooled non-drinker reference
## 6 sexFemale 4.222e-01 1.367e-01 1.304e+00 Arm 2: pooled non-drinker reference
## 7 racecatMexican American 6.810e-01 1.406e-01 3.298e+00 Arm 2: pooled non-drinker reference
## 8 racecatOther Hispanic 3.673e-01 9.112e-02 1.481e+00 Arm 2: pooled non-drinker reference
## 9 racecatNon-Hispanic Black 5.542e-01 1.964e-01 1.564e+00 Arm 2: pooled non-drinker reference
## 10 racecatNon-Hispanic Asian 6.335e-01 1.109e-01 3.619e+00 Arm 2: pooled non-drinker reference
## 11 racecatOther or multiracial 8.332e-02 6.462e-03 1.074e+00 Arm 2: pooled non-drinker reference
## 12 educat9th to 11th grade 4.051e+00 5.815e-01 2.822e+01 Arm 2: pooled non-drinker reference
## 13 educatHigh school or GED 4.059e+00 1.203e+00 1.369e+01 Arm 2: pooled non-drinker reference
## 14 educatSome college 4.963e+00 1.144e+00 2.154e+01 Arm 2: pooled non-drinker reference
## 15 educatCollege graduate 2.503e+00 4.684e-01 1.338e+01 Arm 2: pooled non-drinker reference
## 16 pir 1.029e+00 7.183e-01 1.474e+00 Arm 2: pooled non-drinker reference
## 17 smokecatFormer 1.324e+00 4.807e-01 3.645e+00 Arm 2: pooled non-drinker reference
## 18 smokecatCurrent 1.147e+00 3.518e-01 3.741e+00 Arm 2: pooled non-drinker reference
## 19 srhVery good 3.062e+00 4.170e-01 2.249e+01 Arm 2: pooled non-drinker reference
## 20 srhGood 6.297e+00 8.354e-01 4.747e+01 Arm 2: pooled non-drinker reference
## 21 srhFair 4.293e+01 5.126e+00 3.595e+02 Arm 2: pooled non-drinker reference
## 22 srhPoor 1.989e+00 8.753e-02 4.519e+01 Arm 2: pooled non-drinker reference
## 23 physactYes 5.480e-01 2.304e-01 1.303e+00 Arm 2: pooled non-drinker reference
## 24 1|2 4.843e+02 2.424e+01 9.677e+03 Arm 2: pooled non-drinker reference
## 25 2|3 3.843e+03 1.344e+02 1.099e+05 Arm 2: pooled non-drinker reference
## nrows
## 1 261
## 2 261
## 3 261
## 4 261
## 5 261
## 6 261
## 7 261
## 8 261
## 9 261
## 10 261
## 11 261
## 12 261
## 13 261
## 14 261
## 15 261
## 16 261
## 17 261
## 18 261
## 19 261
## 20 261
## 21 261
## 22 261
## 23 261
## 24 261
## 25 261## term estimate conf.low conf.high
## 1 drinkreFormer, not past heavy 5.883e+00 7.499e-01 4.615e+01
## 2 drinkreFormer, past heavy 2.530e+00 1.564e-01 4.093e+01
## 3 drinkreOccasional 1.719e+00 2.745e-01 1.077e+01
## 4 drinkreLight 1.846e+00 2.187e-01 1.559e+01
## 5 drinkreModerate 1.591e+00 1.515e-01 1.671e+01
## 6 drinkreHeavy 5.652e-01 1.655e-02 1.930e+01
## 7 RIDAGEYR 1.111e+00 1.083e+00 1.140e+00
## 8 sexFemale 3.941e-01 1.318e-01 1.178e+00
## 9 racecatMexican American 9.893e-01 1.643e-01 5.956e+00
## 10 racecatOther Hispanic 3.502e-01 8.391e-02 1.461e+00
## 11 racecatNon-Hispanic Black 5.204e-01 1.870e-01 1.448e+00
## 12 racecatNon-Hispanic Asian 9.647e-01 1.954e-01 4.762e+00
## 13 racecatOther or multiracial 9.304e-02 5.721e-03 1.513e+00
## 14 educat9th to 11th grade 5.313e+00 7.699e-01 3.667e+01
## 15 educatHigh school or GED 5.116e+00 1.445e+00 1.812e+01
## 16 educatSome college 7.018e+00 1.520e+00 3.239e+01
## 17 educatCollege graduate 3.260e+00 6.781e-01 1.567e+01
## 18 pir 9.755e-01 6.819e-01 1.396e+00
## 19 smokecatFormer 1.255e+00 4.476e-01 3.520e+00
## 20 smokecatCurrent 9.184e-01 2.648e-01 3.185e+00
## 21 srhVery good 3.223e+00 4.213e-01 2.465e+01
## 22 srhGood 6.396e+00 8.015e-01 5.104e+01
## 23 srhFair 4.742e+01 5.315e+00 4.231e+02
## 24 srhPoor 1.758e+00 7.340e-02 4.209e+01
## 25 physactYes 4.924e-01 2.019e-01 1.201e+00
## 26 1|2 1.724e+03 5.291e+01 5.615e+04
## 27 2|3 1.457e+04 3.277e+02 6.477e+05
## model nrows
## 1 Arm 3: former drinkers split on past heavy drinking 261
## 2 Arm 3: former drinkers split on past heavy drinking 261
## 3 Arm 3: former drinkers split on past heavy drinking 261
## 4 Arm 3: former drinkers split on past heavy drinking 261
## 5 Arm 3: former drinkers split on past heavy drinking 261
## 6 Arm 3: former drinkers split on past heavy drinking 261
## 7 Arm 3: former drinkers split on past heavy drinking 261
## 8 Arm 3: former drinkers split on past heavy drinking 261
## 9 Arm 3: former drinkers split on past heavy drinking 261
## 10 Arm 3: former drinkers split on past heavy drinking 261
## 11 Arm 3: former drinkers split on past heavy drinking 261
## 12 Arm 3: former drinkers split on past heavy drinking 261
## 13 Arm 3: former drinkers split on past heavy drinking 261
## 14 Arm 3: former drinkers split on past heavy drinking 261
## 15 Arm 3: former drinkers split on past heavy drinking 261
## 16 Arm 3: former drinkers split on past heavy drinking 261
## 17 Arm 3: former drinkers split on past heavy drinking 261
## 18 Arm 3: former drinkers split on past heavy drinking 261
## 19 Arm 3: former drinkers split on past heavy drinking 261
## 20 Arm 3: former drinkers split on past heavy drinking 261
## 21 Arm 3: former drinkers split on past heavy drinking 261
## 22 Arm 3: former drinkers split on past heavy drinking 261
## 23 Arm 3: former drinkers split on past heavy drinking 261
## 24 Arm 3: former drinkers split on past heavy drinking 261
## 25 Arm 3: former drinkers split on past heavy drinking 261
## 26 Arm 3: former drinkers split on past heavy drinking 261
## 27 Arm 3: former drinkers split on past heavy drinking 261Part 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 35.27 7.90 62.64 7
## 2 Era 1: 2011-2016 Former drinker 0 10.31 -4.05 24.66 2
## 3 Era 1: 2011-2016 Occasional 0 20.68 8.60 32.76 14
## 4 Era 1: 2011-2016 Light 0 50.18 27.54 72.83 11
## 5 Era 1: 2011-2016 Moderate 0 45.64 12.45 78.83 4
## 6 Era 1: 2011-2016 Heavy 0 46.72 1.91 91.54 2
## 7 Era 1: 2011-2016 Lifetime abstainer 1 16.92 2.21 31.63 6
## 8 Era 1: 2011-2016 Former drinker 1 1.67 -1.66 5.00 1
## 9 Era 1: 2011-2016 Occasional 1 30.60 14.94 46.27 16
## 10 Era 1: 2011-2016 Light 1 14.58 5.83 23.32 8
## 11 Era 1: 2011-2016 Moderate 1 26.74 -0.25 53.73 4
## 12 Era 1: 2011-2016 Heavy 1 0.00 0.00 0.00 0
## 13 Era 1: 2011-2016 Lifetime abstainer 2 10.35 0.83 19.88 6
## 14 Era 1: 2011-2016 Former drinker 2 27.97 2.89 53.04 6
## 15 Era 1: 2011-2016 Occasional 2 25.95 10.37 41.52 14
## 16 Era 1: 2011-2016 Light 2 23.75 4.39 43.11 6
## 17 Era 1: 2011-2016 Moderate 2 10.22 -6.82 27.27 2
## 18 Era 1: 2011-2016 Heavy 2 31.09 -7.57 69.75 2
## 19 Era 1: 2011-2016 Lifetime abstainer 3 22.80 7.49 38.11 8
## 20 Era 1: 2011-2016 Former drinker 3 28.43 12.11 44.76 14
## 21 Era 1: 2011-2016 Occasional 3 10.75 3.32 18.18 10
## 22 Era 1: 2011-2016 Light 3 9.55 3.10 16.00 6
## 23 Era 1: 2011-2016 Moderate 3 0.00 0.00 0.00 0
## 24 Era 1: 2011-2016 Heavy 3 18.97 -10.35 48.29 2
## 25 Era 1: 2011-2016 Lifetime abstainer 4 14.66 -0.13 29.44 6
## 26 Era 1: 2011-2016 Former drinker 4 31.63 13.82 49.43 11
## 27 Era 1: 2011-2016 Occasional 4 12.02 2.48 21.56 12
## 28 Era 1: 2011-2016 Light 4 1.94 -0.86 4.74 2
## 29 Era 1: 2011-2016 Moderate 4 17.40 -5.15 39.95 2
## 30 Era 1: 2011-2016 Heavy 4 3.21 -3.50 9.92 1
## 31 Era 2: 2017-2020 Lifetime abstainer 0 59.95 5.29 114.60 1
## 32 Era 2: 2017-2020 Former drinker 0 11.13 -4.70 26.96 2
## 33 Era 2: 2017-2020 Occasional 0 13.24 -1.04 27.52 4
## 34 Era 2: 2017-2020 Light 0 38.37 12.44 64.30 7
## 35 Era 2: 2017-2020 Moderate 0 52.40 -12.66 117.46 1
## 36 Era 2: 2017-2020 Heavy 0 40.84 -13.72 95.41 1
## 37 Era 2: 2017-2020 Lifetime abstainer 1 0.00 0.00 0.00 0
## 38 Era 2: 2017-2020 Former drinker 1 4.21 -4.77 13.19 1
## 39 Era 2: 2017-2020 Occasional 1 40.04 7.50 72.59 11
## 40 Era 2: 2017-2020 Light 1 22.05 4.23 39.86 6
## 41 Era 2: 2017-2020 Moderate 1 6.27 -7.87 20.42 1
## 42 Era 2: 2017-2020 Heavy 1 48.71 -5.44 102.86 2
## 43 Era 2: 2017-2020 Lifetime abstainer 2 40.05 -14.60 94.71 3
## 44 Era 2: 2017-2020 Former drinker 2 28.97 -2.08 60.02 5
## 45 Era 2: 2017-2020 Occasional 2 12.99 0.47 25.51 5
## 46 Era 2: 2017-2020 Light 2 14.99 1.23 28.74 7
## 47 Era 2: 2017-2020 Moderate 2 41.33 -22.59 105.25 1
## 48 Era 2: 2017-2020 Heavy 2 5.54 -6.82 17.89 1
## 49 Era 2: 2017-2020 Lifetime abstainer 3 0.00 0.00 0.00 0
## 50 Era 2: 2017-2020 Former drinker 3 26.27 2.18 50.37 9
## 51 Era 2: 2017-2020 Occasional 3 9.37 -1.23 19.98 6
## 52 Era 2: 2017-2020 Light 3 5.79 -3.87 15.45 2
## 53 Era 2: 2017-2020 Moderate 3 0.00 0.00 0.00 0
## 54 Era 2: 2017-2020 Heavy 3 3.19 -4.06 10.44 1
## 55 Era 2: 2017-2020 Lifetime abstainer 4 0.00 0.00 0.00 0
## 56 Era 2: 2017-2020 Former drinker 4 29.42 -5.64 64.48 5
## 57 Era 2: 2017-2020 Occasional 4 24.35 -5.65 54.35 10
## 58 Era 2: 2017-2020 Light 4 18.80 -0.62 38.22 5
## 59 Era 2: 2017-2020 Moderate 4 0.00 0.00 0.00 0
## 60 Era 2: 2017-2020 Heavy 4 1.72 -2.23 5.66 1
## 61 Pooled 2011-2020 Lifetime abstainer 0 37.51 12.35 62.66 8
## 62 Pooled 2011-2020 Former drinker 0 10.66 0.04 21.29 4
## 63 Pooled 2011-2020 Occasional 0 17.82 8.49 27.15 18
## 64 Pooled 2011-2020 Light 0 44.02 26.75 61.29 18
## 65 Pooled 2011-2020 Moderate 0 47.95 16.60 79.30 5
## 66 Pooled 2011-2020 Heavy 0 42.89 4.09 81.69 3
## 67 Pooled 2011-2020 Lifetime abstainer 1 15.39 2.13 28.64 6
## 68 Pooled 2011-2020 Former drinker 1 2.77 -1.46 7.01 2
## 69 Pooled 2011-2020 Occasional 1 34.23 17.89 50.56 27
## 70 Pooled 2011-2020 Light 1 18.47 8.34 28.61 14
## 71 Pooled 2011-2020 Moderate 1 19.74 -0.26 39.73 5
## 72 Pooled 2011-2020 Heavy 1 31.74 -8.21 71.68 2
## 73 Pooled 2011-2020 Lifetime abstainer 2 13.04 2.88 23.20 9
## 74 Pooled 2011-2020 Former drinker 2 28.40 8.85 47.95 11
## 75 Pooled 2011-2020 Occasional 2 20.97 10.06 31.89 19
## 76 Pooled 2011-2020 Light 2 19.18 6.99 31.37 13
## 77 Pooled 2011-2020 Moderate 2 20.87 -5.47 47.21 3
## 78 Pooled 2011-2020 Heavy 2 14.44 -4.62 33.51 3
## 79 Pooled 2011-2020 Lifetime abstainer 3 20.74 6.86 34.61 8
## 80 Pooled 2011-2020 Former drinker 3 27.49 13.50 41.49 23
## 81 Pooled 2011-2020 Occasional 3 10.22 4.07 16.38 16
## 82 Pooled 2011-2020 Light 3 7.59 1.77 13.41 8
## 83 Pooled 2011-2020 Moderate 3 0.00 0.00 0.00 0
## 84 Pooled 2011-2020 Heavy 3 8.69 -4.14 21.53 3
## 85 Pooled 2011-2020 Lifetime abstainer 4 13.33 -0.08 26.74 6
## 86 Pooled 2011-2020 Former drinker 4 30.67 12.53 48.81 16
## 87 Pooled 2011-2020 Occasional 4 16.75 4.34 29.17 22
## 88 Pooled 2011-2020 Light 4 10.74 -0.28 21.75 7
## 89 Pooled 2011-2020 Moderate 4 11.44 -4.20 27.09 2
## 90 Pooled 2011-2020 Heavy 4 2.24 -1.37 5.84 2Table 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 Pooled 2011-2020 ckmtier (1 low, 2 moderate, 3 high)
## 28 Pooled 2011-2020 ckmtier (1 low, 2 moderate, 3 high)
## 29 Pooled 2011-2020 ckmtier (1 low, 2 moderate, 3 high)
## 30 Pooled 2011-2020 ckmtier (1 low, 2 moderate, 3 high)
## 31 Pooled 2011-2020 ckmtier (1 low, 2 moderate, 3 high)
## 54 Check: high risk, stages 3 to 4, pooled high risk, stages 3 to 4
## 55 Check: high risk, stages 3 to 4, pooled high risk, stages 3 to 4
## 56 Check: high risk, stages 3 to 4, pooled high risk, stages 3 to 4
## 57 Check: high risk, stages 3 to 4, pooled high risk, stages 3 to 4
## 58 Check: high risk, stages 3 to 4, pooled high risk, stages 3 to 4
## 78 Sensitivity: occasional drinkers as the reference, pooled ckmtier (1 low, 2 moderate, 3 high)
## 79 Sensitivity: occasional drinkers as the reference, pooled ckmtier (1 low, 2 moderate, 3 high)
## 80 Sensitivity: occasional drinkers as the reference, pooled ckmtier (1 low, 2 moderate, 3 high)
## 81 Sensitivity: occasional drinkers as the reference, pooled ckmtier (1 low, 2 moderate, 3 high)
## 82 Sensitivity: occasional drinkers as the reference, pooled ckmtier (1 low, 2 moderate, 3 high)
## nrows category reference or lcl ucl
## 1 175 Former drinker Lifetime abstainer 13.727 2.237 84.221
## 2 175 Occasional Lifetime abstainer 1.914 0.401 9.131
## 3 175 Light Lifetime abstainer 2.835 0.369 21.779
## 4 175 Moderate Lifetime abstainer 0.717 0.040 12.903
## 5 175 Heavy Lifetime abstainer 21.853 1.448 329.766
## 27 261 Former drinker Lifetime abstainer 5.709 0.735 44.344
## 28 261 Occasional Lifetime abstainer 1.705 0.271 10.729
## 29 261 Light Lifetime abstainer 1.838 0.217 15.592
## 30 261 Moderate Lifetime abstainer 1.583 0.150 16.669
## 31 261 Heavy Lifetime abstainer 0.575 0.017 19.169
## 54 261 Former drinker Lifetime abstainer 2.178 0.362 13.092
## 55 261 Occasional Lifetime abstainer 1.065 0.211 5.384
## 56 261 Light Lifetime abstainer 0.517 0.049 5.501
## 57 261 Moderate Lifetime abstainer 2.148 0.041 112.355
## 58 261 Heavy Lifetime abstainer 0.267 0.017 4.131
## 78 261 Lifetime abstainer Occasional drinker 0.587 0.093 3.692
## 79 261 Former drinker Occasional drinker 3.349 1.011 11.092
## 80 261 Light Occasional drinker 1.078 0.329 3.534
## 81 261 Moderate Occasional drinker 0.929 0.280 3.080
## 82 261 Heavy Occasional drinker 0.337 0.016 7.175Table 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) 175 Former drinker Lifetime abstainer 4.865 0.772 30.651
## 2 Era 1: 2011-2016 ckmstage (0 to 4) 175 Occasional Lifetime abstainer 1.774 0.345 9.121
## 3 Era 1: 2011-2016 ckmstage (0 to 4) 175 Light Lifetime abstainer 1.286 0.158 10.481
## 4 Era 1: 2011-2016 ckmstage (0 to 4) 175 Moderate Lifetime abstainer 0.695 0.063 7.627
## 5 Era 1: 2011-2016 ckmstage (0 to 4) 175 Heavy Lifetime abstainer 2.105 0.110 40.190
## 29 Pooled 2011-2020 ckmstage (0 to 4) 261 Former drinker Lifetime abstainer 3.400 0.445 25.965
## 30 Pooled 2011-2020 ckmstage (0 to 4) 261 Occasional Lifetime abstainer 1.979 0.285 13.739
## 31 Pooled 2011-2020 ckmstage (0 to 4) 261 Light Lifetime abstainer 1.532 0.163 14.417
## 32 Pooled 2011-2020 ckmstage (0 to 4) 261 Moderate Lifetime abstainer 1.148 0.093 14.214
## 33 Pooled 2011-2020 ckmstage (0 to 4) 261 Heavy Lifetime abstainer 0.660 0.027 16.241write.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 1.7048 0.270884 10.729 Arm 1: lifetime abstainers 261 Occasional
## 3 drinkcatLight 1.8385 0.216777 15.592 Arm 1: lifetime abstainers 261 Light
## 4 drinkcatModerate 1.5830 0.150329 16.669 Arm 1: lifetime abstainers 261 Moderate
## 5 drinkcatHeavy 0.5747 0.017232 19.169 Arm 1: lifetime abstainers 261 Heavy
## 27 nondrinkOccasional 0.4667 0.146662 1.485 Arm 2: pooled non-drinkers 261 Occasional
## 28 nondrinkLight 0.5166 0.122271 2.182 Arm 2: pooled non-drinkers 261 Light
## 29 nondrinkModerate 0.4476 0.083080 2.411 Arm 2: pooled non-drinkers 261 Moderate
## 30 nondrinkHeavy 0.1602 0.007719 3.323 Arm 2: pooled non-drinkers 261 Heavy
## 54 drinkreOccasional 1.7194 0.274488 10.771 Arm 3: former drinkers split 261 Occasional
## 55 drinkreLight 1.8465 0.218684 15.591 Arm 3: former drinkers split 261 Light
## 56 drinkreModerate 1.5910 0.151467 16.712 Arm 3: former drinkers split 261 Moderate
## 57 drinkreHeavy 0.5652 0.016553 19.299 Arm 3: former drinkers split 261 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:08:29"
## [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): 1"
## [7] "Exclusions: rows as read 300"
## [8] "Exclusions: aged 20 years and older 300"
## [9] "Exclusions: and not pregnant at examination 300"
## [10] "Exclusions: and in the fasting subsample 300"
## [11] "Exclusions: and CKM stage could be determined 300"
## [12] "Arm rows with a drinking category: 283"
## [13] "Arm rows with every covariate observed: 261"
## [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 = 0"
## [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