1 What this analysis does

Design. Cross-sectional analysis of individual survey responses from the CDC Behavioral Risk Factor Surveillance System (BRFSS), Arkansas, 2022 and 2024 combined.

Population. Arkansas adults aged 45-75, the ages for which colorectal cancer screening is recommended.

Primary outcome. Up to date with colorectal cancer screening (yes or no), using CDC’s calculated variable for meeting U.S. Preventive Services Task Force recommendations.

Model. Because the outcome is yes/no, the model is logistic regression. Because several predictors are entered together, it is multivariable logistic regression. Because BRFSS is a weighted complex survey, it is fitted with survey weights, strata and sampling units (svyglm in the survey package). Results are adjusted odds ratios with 95% confidence intervals.

Effect modification. Three interactions chosen before looking at results: age group x income, age group x sex, and age group x survey year.

2 Data

The first knit downloads two public files from CDC and keeps only Arkansas. This takes several minutes. Later knits reuse the small saved files.

keep <- paste0("^(_STATE|_PSU|_STSTR|_LLCPWT|_SEX|_IMPRACE|_AGEG5YR|_AGE80|",
               "_EDUCAG|_INCOMG|_HLTHPL|PERSDOC|MEDCOST|CHECKUP|_METSTAT|",
               "_URBSTAT|_SMOKER|_BMI5CAT|_CRCREC|_CLNSCP|_SGMSCP|_RFBLDS|",
               "_STOLDN|_VIRCOL|_SBONTI)")

get_brfss_ar <- function(year) {
  cache <- sprintf("brfss_ar_%d.rds", year)
  if (file.exists(cache)) return(readRDS(cache))

  zipf <- sprintf("LLCP%dXPT.zip", year)
  if (!file.exists(zipf)) {
    download.file(
      sprintf("https://www.cdc.gov/brfss/annual_data/%d/files/LLCP%dXPT.zip",
              year, year),
      zipf, mode = "wb")
  }
  exdir <- sprintf("brfss_%d_unzipped", year)
  dir.create(exdir, showWarnings = FALSE)
  unzip(zipf, exdir = exdir)
  xpt <- list.files(exdir, pattern = "XPT", full.names = TRUE,
                    ignore.case = TRUE)[1]

  # Read only the columns needed; fall back to reading everything
  d <- tryCatch(
    haven::read_xpt(xpt, col_select = tidyselect::matches(keep)),
    error = function(e) haven::read_xpt(xpt))
  d <- as.data.frame(d)
  if (!"_STATE" %in% names(d)) stop("State variable not found in ", xpt)

  d <- d[d[["_STATE"]] == 5, grepl(keep, names(d)), drop = FALSE]  # 5 = Arkansas
  d$year <- year
  saveRDS(d, cache)
  unlink(exdir, recursive = TRUE)   # remove the 1 GB unzipped file
  d
}

ar22 <- get_brfss_ar(2022)
ar24 <- get_brfss_ar(2024)

c(Arkansas_2022 = nrow(ar22), Arkansas_2024 = nrow(ar24))

## Arkansas_2022 Arkansas_2024 
##          5309          5344

CDC’s codebooks list 5,309 Arkansas respondents in 2022 and 5,344 in 2024. The two numbers above should match those.

2.1 Harmonise variable names across years

Some variable names change suffix between years (for example _CRCREC2 in 2022 and _CRCREC3 in 2024), so each is picked up by its prefix.

pick <- function(d, prefix) {
  # An exact name wins (so "_LLCPWT", the final weight, is never confused
  # with "_LLCPWT2", the design weight before raking).
  if (prefix %in% names(d)) return(as.numeric(d[[prefix]]))
  nm <- grep(paste0("^", prefix), names(d), value = TRUE)
  if (length(nm) == 0) stop("No variable starting with ", prefix)
  if (length(nm) > 1) stop("More than one variable starts with ", prefix, ": ",
                           paste(nm, collapse = ", "))
  as.numeric(d[[nm[1]]])
}

harmonise <- function(d) {
  data.frame(
    year   = d$year,
    psu    = pick(d, "_PSU"),
    strat  = pick(d, "_STSTR"),
    wt     = pick(d, "_LLCPWT"),
    age    = pick(d, "_AGE80"),
    sex    = pick(d, "_SEX"),
    race   = pick(d, "_IMPRACE"),
    educ   = pick(d, "_EDUCAG"),
    income = pick(d, "_INCOMG"),
    ins    = pick(d, "_HLTHPL"),
    doc    = pick(d, "PERSDOC"),
    metro  = pick(d, "_METSTAT"),
    smoker = pick(d, "_SMOKER"),
    crc    = pick(d, "_CRCREC"),   # 1 up to date, 2 overdue, 3 never screened
    colo   = pick(d, "_CLNSCP"),   # 1 colonoscopy within 10 years
    stool  = pick(d, "_RFBLDS"),   # 1 stool blood test within 1 year
    sdna   = pick(d, "_STOLDN")    # 1 stool DNA test within 3 years
  )
}

# These must explicitly refer to the FINAL raked BRFSS weight.
stopifnot("_LLCPWT" %in% names(ar22), "_LLCPWT" %in% names(ar24))
d <- rbind(harmonise(ar22), harmonise(ar24))
stopifnot(
  isTRUE(all.equal(d$wt[d$year == 2022], as.numeric(ar22[["_LLCPWT"]]))),
  isTRUE(all.equal(d$wt[d$year == 2024], as.numeric(ar24[["_LLCPWT"]])))
)

2.2 Recode

yes_no <- function(x) ifelse(x == 1, 1, ifelse(x %in% c(2, 3), 0, NA))

d$uptodate <- yes_no(d$crc)
d$never    <- ifelse(d$crc == 3, 1, ifelse(d$crc %in% c(1, 2), 0, NA))
d$colo10   <- yes_no(d$colo)
d$stool1   <- yes_no(d$stool)
d$sdna3    <- yes_no(d$sdna)

d$agegrp <- cut(d$age, breaks = c(44, 54, 64, 75),
                labels = c("45-54", "55-64", "65-75"))
d$agegrp <- relevel(d$agegrp, ref = "65-75")

d$sex <- factor(d$sex, levels = c(2, 1), labels = c("Female", "Male"))

race_lab <- rep("Other", nrow(d))
race_lab[which(d$race == 1)] <- "White"
race_lab[which(d$race == 2)] <- "Black"
race_lab[which(d$race == 5)] <- "Hispanic"
race_lab[is.na(d$race)] <- NA
d$race <- factor(race_lab, levels = c("White", "Black", "Hispanic", "Other"))

d$educ <- factor(ifelse(d$educ %in% 1:4, d$educ, NA), levels = c(4, 3, 2, 1),
                 labels = c("College graduate", "Some college",
                            "High school", "Less than high school"))

inc_lab <- rep(NA_character_, nrow(d))
inc_lab[d$income %in% 5:7] <- "$50,000 or more"
inc_lab[d$income %in% 3:4] <- "$25,000-49,999"
inc_lab[d$income %in% 1:2] <- "Under $25,000"
inc_lab[d$income %in% 9]   <- "Not reported"
d$income3 <- factor(inc_lab, levels = c("$50,000 or more", "$25,000-49,999",
                                        "Under $25,000", "Not reported"))

d$insured  <- factor(ifelse(d$ins == 1, "Insured",
                     ifelse(d$ins == 2, "Uninsured", NA)),
                     levels = c("Insured", "Uninsured"))
d$provider <- factor(ifelse(d$doc %in% 1:2, "Has provider",
                     ifelse(d$doc == 3, "No provider", NA)),
                     levels = c("Has provider", "No provider"))
d$metro    <- factor(ifelse(d$metro == 1, "Metropolitan",
                     ifelse(d$metro == 2, "Nonmetropolitan", NA)),
                     levels = c("Metropolitan", "Nonmetropolitan"))
d$smoke    <- factor(ifelse(d$smoker %in% 1:2, "Current",
                     ifelse(d$smoker == 3, "Former",
                     ifelse(d$smoker == 4, "Never", NA))),
                     levels = c("Never", "Former", "Current"))
d$yearf <- factor(d$year)

# Save the Arkansas analysis file (small; this is the file to share)
write.csv(d, "brfss_arkansas_crc_2022_2024.csv", row.names = FALSE)

2.3 Survey design

d$strat2 <- interaction(d$year, d$strat, drop = TRUE)  # strata within year
d$wt2    <- d$wt / 2                                   # two years pooled

des_all <- svydesign(ids = ~psu, strata = ~strat2, weights = ~wt2,
                     data = d, nest = TRUE)
des <- subset(des_all, !is.na(uptodate))

cat("Original Arkansas records by year:\n")

## Original Arkansas records by year:

print(table(d$year))

## 
## 2022 2024 
## 5309 5344

cat("Screening-valid records by year:\n")

## Screening-valid records by year:

print(table(Year = d$year[!is.na(d$uptodate)]))

## Year
## 2022 2024 
## 2654 2688

cat("Missing screening outcome among age 45–75 by year:\n")

## Missing screening outcome among age 45–75 by year:

print(with(d[!is.na(d$agegrp), ], table(year, is.na(uptodate))))

##       
## year   FALSE TRUE
##   2022  2654  454
##   2024  2688  250

3 Check against CDC’s published figures

The first column is calculated with the final BRFSS survey weight (_LLCPWT). The comparison values below are independently reproduced from the original 2022 and 2024 CDC XPT data; they are not a substitute for verifying your code. The script stops if the reproduced estimates differ by more than 0.15 percentage points.

chk <- svyby(~uptodate, ~yearf, des, svymean, na.rm = TRUE)
ref_rates <- c(`2022` = 64.0279745939, `2024` = 68.0846445592)
actual_rates <- setNames(100 * chk$uptodate, as.character(chk$yearf))
stopifnot(all(names(ref_rates) %in% names(actual_rates)))
if (any(abs(actual_rates[names(ref_rates)] - ref_rates) > 0.15)) {
  stop("Screening prevalence does not match the independent CDC XPT check. Check _LLCPWT selection and original data files.")
}
kable(data.frame(
  Year = chk$yearf,
  `Weighted screening (%)` = round(100 * chk$uptodate, 1),
  `Independent XPT verification (%)` = round(ref_rates[as.character(chk$yearf)], 1),
  check.names = FALSE
))

	Year	Weighted screening (%)	Independent XPT verification (%)
2022	2022	64.0	64.0
2024	2024	68.1	68.1

4 Who is up to date? Weighted percentages

prev_by <- function(var, outcome = "uptodate", design = des) {
  s <- svyby(as.formula(paste0("~", outcome)), as.formula(paste0("~", var)),
             design, svymean, vartype = "ci", na.rm = TRUE)
  n <- table(d[[var]][!is.na(d[[outcome]])])
  data.frame(Characteristic = var,
             Group   = as.character(s[[var]]),
             n       = as.integer(n[as.character(s[[var]])]),
             Percent = round(100 * s[[outcome]], 1),
             Lower   = round(100 * s$ci_l, 1),
             Upper   = round(100 * s$ci_u, 1))
}
vars <- c("yearf", "agegrp", "sex", "race", "educ", "income3",
          "insured", "provider", "metro", "smoke")
kable(do.call(rbind, lapply(vars, prev_by)), row.names = FALSE,
      caption = "Up to date with colorectal cancer screening, Arkansas adults aged 45-75")

Up to date with colorectal cancer screening, Arkansas adults aged 45-75
Characteristic	Group	n	Percent	Lower	Upper
yearf	2022	2654	64.0	61.8	66.3
yearf	2024	2688	68.1	65.9	70.2
agegrp	65-75	2305	79.7	77.7	81.7
agegrp	45-54	1313	45.4	42.3	48.6
agegrp	55-64	1724	72.3	69.8	74.8
sex	Female	2943	68.4	66.3	70.5
sex	Male	2399	63.8	61.4	66.1
race	White	4312	68.1	66.4	69.7
race	Black	662	66.1	61.5	70.7
race	Hispanic	118	37.4	27.4	47.4
race	Other	250	60.6	52.8	68.4
educ	College graduate	1810	71.3	68.8	73.8
educ	Some college	1481	69.4	66.6	72.2
educ	High school	1669	62.7	59.9	65.4
educ	Less than high school	371	57.2	51.3	63.0
income3	$50,000 or more	2216	68.4	66.0	70.7
income3	$25,000-49,999	1242	67.0	63.7	70.3
income3	Under $25,000	1068	62.2	58.6	65.8
income3	Not reported	816	63.8	59.7	67.9
insured	Insured	4911	68.4	66.8	70.0
insured	Uninsured	192	29.4	21.7	37.2
provider	Has provider	4992	69.2	67.7	70.8
provider	No provider	312	27.8	22.0	33.6
metro	Metropolitan	2843	68.1	66.0	70.2
metro	Nonmetropolitan	2499	63.5	61.1	65.8
smoke	Never	2730	67.3	65.2	69.5
smoke	Former	1567	69.5	66.6	72.4
smoke	Current	974	57.9	54.2	61.7

5 Which tests, by income

test_by_income <- function(v, label) {
  ok  <- !is.na(des_all$variables[[v]]) & !is.na(des_all$variables$income3)
  dsg <- des_all[ok, ]
  s <- svyby(as.formula(paste0("~", v)), ~income3, dsg, svymean,
             vartype = "ci", na.rm = TRUE)
  data.frame(Test = label, Income = as.character(s$income3),
             Percent = round(100 * s[[v]], 1),
             Lower = round(100 * s$ci_l, 1), Upper = round(100 * s$ci_u, 1))
}
kable(rbind(test_by_income("colo10", "Colonoscopy within 10 years"),
            test_by_income("stool1", "Stool blood test within 1 year"),
            test_by_income("sdna3",  "Stool DNA test within 3 years")),
      row.names = FALSE)

Test	Income	Percent	Lower	Upper
Colonoscopy within 10 years	$50,000 or more	64.2	61.9	66.6
Colonoscopy within 10 years	$25,000-49,999	62.5	59.1	66.0
Colonoscopy within 10 years	Under $25,000	54.7	51.1	58.4
Colonoscopy within 10 years	Not reported	58.8	54.6	63.0
Stool blood test within 1 year	$50,000 or more	6.8	5.6	8.0
Stool blood test within 1 year	$25,000-49,999	7.1	5.5	8.8
Stool blood test within 1 year	Under $25,000	9.9	7.6	12.1
Stool blood test within 1 year	Not reported	6.0	4.2	7.7
Stool DNA test within 3 years	$50,000 or more	3.2	2.4	4.1
Stool DNA test within 3 years	$25,000-49,999	4.7	2.8	6.6
Stool DNA test within 3 years	Under $25,000	5.6	3.8	7.3
Stool DNA test within 3 years	Not reported	3.6	2.1	5.1

6 Multivariable logistic regression

Odds ratios below 1 mean lower odds of being up to date than the reference group. Reference groups: age 65-75, female, White, 2022, college graduate, income of $50,000 or more, insured, has a provider, metropolitan, never smoked.

or_table <- function(fit) {
  co <- coef(summary(fit))
  ci <- confint(fit)
  out <- data.frame(Term = rownames(co),
                    OR = exp(co[, "Estimate"]),
                    Lower = exp(ci[, 1]), Upper = exp(ci[, 2]),
                    p = co[, ncol(co)], row.names = NULL)
  out[out$Term != "(Intercept)", ]
}

# Model 1: demographics
m1 <- svyglm(uptodate ~ agegrp + sex + race + yearf,
             design = des, family = quasibinomial())
# Model 2: + socioeconomic position
m2 <- svyglm(uptodate ~ agegrp + sex + race + yearf + educ + income3,
             design = des, family = quasibinomial())
# Model 3: + access to care, residence and smoking (main model)
m3 <- svyglm(uptodate ~ agegrp + sex + race + yearf + educ + income3 +
               insured + provider + metro + smoke,
             design = des, family = quasibinomial())

kable(or_table(m1), digits = 3, row.names = FALSE, caption = "Model 1: demographics")

Model 1: demographics
Term	OR	Lower	Upper	p
agegrp45-54	0.218	0.183	0.261	0.000
agegrp55-64	0.678	0.568	0.810	0.000
sexMale	0.805	0.695	0.933	0.004
raceBlack	0.994	0.794	1.244	0.960
raceHispanic	0.368	0.238	0.570	0.000
raceOther	0.878	0.616	1.251	0.471
yearf2024	1.263	1.091	1.462	0.002

kable(or_table(m2), digits = 3, row.names = FALSE, caption = "Model 2: + education and income")

Model 2: + education and income
Term	OR	Lower	Upper	p
agegrp45-54	0.203	0.168	0.244	0.000
agegrp55-64	0.665	0.556	0.795	0.000
sexMale	0.816	0.702	0.948	0.008
raceBlack	1.058	0.846	1.323	0.623
raceHispanic	0.409	0.262	0.637	0.000
raceOther	0.901	0.629	1.291	0.570
yearf2024	1.226	1.058	1.422	0.007
educSome college	0.890	0.734	1.080	0.237
educHigh school	0.653	0.540	0.790	0.000
educLess than high school	0.589	0.434	0.799	0.001
income3$25,000-49,999	0.858	0.698	1.055	0.146
income3Under $25,000	0.791	0.635	0.984	0.036
income3Not reported	0.723	0.573	0.913	0.006

kable(or_table(m3), digits = 3, row.names = FALSE, caption = "Model 3 (main): + insurance, provider, residence, smoking")

Model 3 (main): + insurance, provider, residence, smoking
Term	OR	Lower	Upper	p
agegrp45-54	0.233	0.191	0.284	0.000
agegrp55-64	0.753	0.624	0.909	0.003
sexMale	0.854	0.729	1.002	0.052
raceBlack	1.070	0.837	1.368	0.589
raceHispanic	0.526	0.317	0.875	0.013
raceOther	0.885	0.600	1.307	0.540
yearf2024	1.228	1.051	1.436	0.010
educSome college	0.960	0.784	1.176	0.696
educHigh school	0.746	0.609	0.915	0.005
educLess than high school	0.682	0.488	0.953	0.025
income3$25,000-49,999	0.976	0.785	1.213	0.826
income3Under $25,000	0.890	0.703	1.127	0.334
income3Not reported	0.781	0.605	1.007	0.057
insuredUninsured	0.440	0.272	0.712	0.001
providerNo provider	0.279	0.197	0.396	0.000
metroNonmetropolitan	0.724	0.618	0.848	0.000
smokeFormer	1.093	0.909	1.314	0.347
smokeCurrent	0.844	0.680	1.047	0.124

7 Effect modification

Each interaction is added to the main model and tested with a Wald test.

base_terms <- "sex + race + yearf + educ + income3 + insured + provider + metro + smoke"
int_test <- function(term, label) {
  f   <- as.formula(paste("uptodate ~ agegrp +", base_terms, "+", term))
  fit <- svyglm(f, design = des, family = quasibinomial())
  tt  <- regTermTest(fit, as.formula(paste("~", term)))
  data.frame(Interaction = label, F = as.numeric(tt$Ftest),
             df = as.numeric(tt$df), p = as.numeric(tt$p))
}
kable(rbind(int_test("agegrp:income3", "Age group x income"),
            int_test("agegrp:sex",     "Age group x sex"),
            int_test("agegrp:yearf",   "Age group x survey year")),
      digits = 3, row.names = FALSE)

Interaction	F	df	p
Age group x income	1.136	6	0.339
Age group x sex	1.921	2	0.147
Age group x survey year	0.180	2	0.835

7.1 Income within each age group

strat_or <- do.call(rbind, lapply(c("45-54", "55-64", "65-75"), function(g) {
  dg  <- subset(des, agegrp == g)
  fit <- svyglm(uptodate ~ income3 + sex + race + yearf + educ + insured,
                design = dg, family = quasibinomial())
  tab <- or_table(fit)
  cbind(`Age group` = g, tab[grepl("^income3", tab$Term), ])
}))
kable(strat_or, digits = 3, row.names = FALSE,
      caption = "Adjusted odds ratios for income, by age group")

Adjusted odds ratios for income, by age group
Age group	Term	OR	Lower	Upper	p
45-54	income3$25,000-49,999	0.837	0.562	1.246	0.381
45-54	income3Under $25,000	0.856	0.585	1.253	0.424
45-54	income3Not reported	0.984	0.604	1.601	0.947
55-64	income3$25,000-49,999	0.937	0.656	1.338	0.722
55-64	income3Under $25,000	0.881	0.607	1.278	0.503
55-64	income3Not reported	0.714	0.476	1.070	0.102
65-75	income3$25,000-49,999	0.717	0.497	1.035	0.076
65-75	income3Under $25,000	0.523	0.343	0.797	0.003
65-75	income3Not reported	0.531	0.361	0.782	0.001

7.2 Change from 2022 to 2024 within each age group

s <- svyby(~uptodate, ~agegrp + yearf, des, svymean, vartype = "ci", na.rm = TRUE)
kable(data.frame(`Age group` = s$agegrp, Year = s$yearf,
                 Percent = round(100 * s$uptodate, 1),
                 Lower = round(100 * s$ci_l, 1), Upper = round(100 * s$ci_u, 1),
                 check.names = FALSE), row.names = FALSE)

Age group	Year	Percent	Lower	Upper
65-75	2022	78.7	75.8	81.6
45-54	2022	41.9	37.4	46.3
55-64	2022	69.3	65.5	73.1
65-75	2024	80.6	77.8	83.4
45-54	2024	48.4	44.0	52.9
55-64	2024	75.2	72.0	78.5

8 Secondary outcome: never screened

des_n <- subset(des_all, !is.na(never))
mn <- svyglm(never ~ agegrp + sex + race + yearf + educ + income3 +
               insured + provider + metro + smoke,
             design = des_n, family = quasibinomial())
kable(or_table(mn), digits = 3, row.names = FALSE,
      caption = "Adjusted odds ratios for never having been screened")

Adjusted odds ratios for never having been screened
Term	OR	Lower	Upper	p
agegrp45-54	6.392	5.123	7.975	0.000
agegrp55-64	1.536	1.228	1.921	0.000
sexMale	1.264	1.060	1.508	0.009
raceBlack	1.056	0.807	1.381	0.691
raceHispanic	2.065	1.224	3.482	0.007
raceOther	0.997	0.650	1.528	0.989
yearf2024	0.744	0.625	0.886	0.001
educSome college	1.037	0.826	1.301	0.755
educHigh school	1.357	1.083	1.701	0.008
educLess than high school	1.508	1.046	2.175	0.028
income3$25,000-49,999	1.013	0.795	1.292	0.914
income3Under $25,000	1.173	0.906	1.518	0.226
income3Not reported	1.222	0.918	1.625	0.169
insuredUninsured	2.231	1.355	3.673	0.002
providerNo provider	3.417	2.390	4.884	0.000
metroNonmetropolitan	1.280	1.073	1.527	0.006
smokeFormer	0.898	0.733	1.102	0.303
smokeCurrent	1.125	0.884	1.433	0.339

9 Figure

s <- svyby(~uptodate, ~agegrp + income3, des, svymean, vartype = "ci", na.rm = TRUE)
s <- s[s$income3 != "Not reported", ]
s$income3 <- droplevels(s$income3)
s$agegrp  <- factor(s$agegrp, levels = c("45-54", "55-64", "65-75"))
cols <- c("#762a83", "#e08214", "#1b7837")
k <- as.numeric(s$income3)
x <- as.numeric(s$agegrp) + (k - 2) * 0.15
plot(x, 100 * s$uptodate, ylim = c(0, 100), xlim = c(0.5, 3.5), xaxt = "n",
     pch = 19, col = cols[k], xlab = "Age group",
     ylab = "Up to date with screening (%)",
     main = "Colorectal cancer screening in Arkansas by age and income, 2022-2024")
axis(1, at = 1:3, labels = c("45-54", "55-64", "65-75"))
segments(x, 100 * s$ci_l, x, 100 * s$ci_u, col = cols[k], lwd = 2)
legend("bottomright", legend = levels(s$income3), col = cols, pch = 19,
       bty = "n", title = "Household income")

10 Primary model sample counts

cat("Valid screening outcomes:", sum(!is.na(d$uptodate)), "\n")

## Valid screening outcomes: 5342

cat("Full adjusted regression sample:", nobs(m3), "\n")

## Full adjusted regression sample: 4999

11 Final validation: age comparability and missing screening data

These are sensitivity analyses. The primary outcome uses the official CDC _CRCREC2 (2022) and _CRCREC3 (2024) screening variables and the final BRFSS raked weight _LLCPWT. CDC documents that 2022 _CRCREC2 uses reported AGE, whereas 2024 _CRCREC3 uses _IMPAGE. Both target adults ages 45–75. No sensitivity analysis can make these two derived variables identical; we report the change as a limitation.

11.1 Screening by year in narrower age ranges

If the estimated time trend is comparable when excluding participants near the lower/upper boundaries, this supports robustness to the chosen age range, but does not prove that the outcome algorithms are identical.

age_ranges <- list(`45-75 (primary)` = c(45,75),
                   `50-75` = c(50,75),
                   `50-74` = c(50,74),
                   `55-74` = c(55,74))

age_res <- do.call(rbind, lapply(names(age_ranges), function(label) {
  bounds <- age_ranges[[label]]
  age_des <- subset(des, age >= bounds[1] & age <= bounds[2])
  x <- svyby(~uptodate, ~yearf, age_des, svymean,
             vartype = "ci", na.rm = TRUE)
  data.frame(Age_range = label, Year = as.character(x$yearf),
             Rate_pct = round(100 * x$uptodate, 2),
             CI_low = round(100 * x$ci_l, 2),
             CI_high = round(100 * x$ci_u, 2),
             stringsAsFactors = FALSE)
}))
kable(age_res, caption = "Age-range sensitivity, with correct final BRFSS survey weights")

Age-range sensitivity, with correct final BRFSS survey weights
Age_range	Year	Rate_pct	CI_low	CI_high
45-75 (primary)	2022	64.03	61.75	66.30
45-75 (primary)	2024	68.08	65.93	70.24
50-75	2022	69.20	66.87	71.54
50-75	2024	74.05	71.93	76.18
50-74	2022	68.78	66.39	71.17
50-74	2024	74.02	71.86	76.19
55-74	2022	73.51	71.00	76.03
55-74	2024	77.99	75.81	80.18

cat("\nPercentage-point changes (2024 minus 2022):\n")

## 
## Percentage-point changes (2024 minus 2022):

for (label in names(age_ranges)) {
  vals <- age_res[age_res$Age_range == label, ]
  cat(label, ": ", round(vals$Rate_pct[vals$Year == "2024"] -
                         vals$Rate_pct[vals$Year == "2022"], 2), "\n", sep="")
}

## 45-75 (primary): 4.05
## 50-75: 4.85
## 50-74: 5.24
## 55-74: 4.48

11.2 Missing-outcome patterns

The counts below are unweighted sample characteristics, not estimated population prevalences. Missing screening outcomes can reflect nonresponse, question skip patterns or incomplete information for the CDC-derived outcome; we do not assume missing means unscreened.

audit <- d[!is.na(d$agegrp), , drop = FALSE]
audit$age_detail <- cut(audit$age, c(44,49,54,64,75),
                        labels = c("45-49", "50-54", "55-64", "65-75"))
audit$observed <- as.integer(!is.na(audit$uptodate))
cat("\nOutcome availability by year (unweighted):\n")

## 
## Outcome availability by year (unweighted):

print(with(audit, table(Year = year, Observed = observed)))

##       Observed
## Year      0    1
##   2022  454 2654
##   2024  250 2688

cat("\nOutcome availability by age and year (unweighted):\n")

## 
## Outcome availability by age and year (unweighted):

print(with(audit, table(Year = year, Age = age_detail, Observed = observed)))

## , , Observed = 0
## 
##       Age
## Year   45-49 50-54 55-64 65-75
##   2022    51   111   131   161
##   2024    25    34    93    98
## 
## , , Observed = 1
## 
##       Age
## Year   45-49 50-54 55-64 65-75
##   2022   287   356   845  1166
##   2024   323   347   879  1139

11.3 Inverse-probability weighting for observed screening outcomes

This adjusts for measured predictors of having a usable screening response. It assumes missingness is conditionally independent of screening after adjusting for these variables and cannot correct for unmeasured nonresponse, changing question content or outcome-definition differences. We estimate observation probabilities using the corrected BRFSS survey weights, not the mistakenly chosen _LLCPWT2 weight.

# Keep all age-eligible records when modeling whether their outcome is known.
# Explicitly account for nonresponse on education as a separate category.
audit$edu_audit <- factor(ifelse(is.na(audit$educ), "Missing",
                                  as.character(audit$educ)))

obs_design <- svydesign(ids = ~psu, strata = ~strat2,
                        weights = ~wt2, data = audit, nest = TRUE)

obs_fit <- svyglm(observed ~ yearf * age_detail + sex + race +
                    metro + income3 + edu_audit,
                  design = obs_design, family = quasibinomial())

audit$p_observed <- as.numeric(predict(obs_fit, newdata = audit,
                                      type = "response"))
if (anyNA(audit$p_observed) ||
    any(!is.finite(audit$p_observed)) ||
    any(audit$p_observed <= 0 | audit$p_observed > 1)) {
  stop("Invalid observation propensity: review model predictors")
}

# Multiply final survey weights by inverse propensity ONLY for observed cases.
audit$ipw_weight <- audit$wt2
is_known <- audit$observed == 1
audit$ipw_weight[is_known] <- audit$wt2[is_known] /
                              audit$p_observed[is_known]

orig_obs <- subset(obs_design, observed == 1)
ipw_full <- svydesign(ids = ~psu, strata = ~strat2,
                      weights = ~ipw_weight, data = audit, nest = TRUE)
ipw_obs <- subset(ipw_full, observed == 1)

orig <- svyby(~uptodate, ~yearf, orig_obs, svymean,
              vartype = "ci", na.rm = TRUE)
adj <- svyby(~uptodate, ~yearf, ipw_obs, svymean,
             vartype = "ci", na.rm = TRUE)

ipw_comparison <- data.frame(
  Year = as.character(orig$yearf),
  Original_pct = round(100 * orig$uptodate, 2),
  Original_LCL = round(100 * orig$ci_l, 2),
  Original_UCL = round(100 * orig$ci_u, 2),
  IPW_pct = round(100 * adj$uptodate, 2),
  IPW_LCL = round(100 * adj$ci_l, 2),
  IPW_UCL = round(100 * adj$ci_u, 2)
)
kable(ipw_comparison, caption = "Missingness sensitivity with final raked weight")

Missingness sensitivity with final raked weight
Year	Original_pct	Original_LCL	Original_UCL	IPW_pct	IPW_LCL	IPW_UCL
2022	64.03	61.76	66.30	63.31	60.89	65.73
2024	68.08	65.93	70.24	67.75	65.57	69.93

cat("\nOriginal percentage-point change:",
    round(100 * (orig$uptodate[2] - orig$uptodate[1]), 2), "\n")

## 
## Original percentage-point change: 4.06

cat("IPW-adjusted percentage-point change:",
    round(100 * (adj$uptodate[2] - adj$uptodate[1]), 2), "\n")

## IPW-adjusted percentage-point change: 4.44

cat("\nObserved-case inverse observation probability (weight multiplier):\n")

## 
## Observed-case inverse observation probability (weight multiplier):

print(summary(1 / audit$p_observed[is_known]))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.023   1.047   1.062   1.128   1.102   4.916

cat("\nNumber of complete screening outcomes:", sum(is_known), "\n")

## 
## Number of complete screening outcomes: 5342

# The CIs above use the estimated propensities as fixed and may understate
# uncertainty because they do not re-estimate the propensity model in each
# survey replicate. Interpret as a secondary sensitivity analysis.

11.4 Interpretation for abstract

Report the survey-weighted screening rates and adjusted odds ratios from the primary analysis above. If the missingness-adjusted and narrower-age analyses show similarly directed trends, describe them briefly as sensitivity analyses, not as proof of causality or equivalence of CDC variable coding. The change between survey years is observational and cannot be attributed to Arkansas’s policy changes on these data alone.

12 Session information

sessionInfo()

## R version 4.6.2 (2026-10-29)
## Platform: aarch64-apple-darwin23
## Running under: macOS Golden Gate 27.0
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Chicago
## tzcode source: internal
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] knitr_1.52      survey_4.5      survival_3.8-12 Matrix_1.7-6   
## [5] haven_2.5.5    
## 
## loaded via a namespace (and not attached):
##  [1] vctrs_0.7.3       cli_3.6.6         rlang_1.3.0       xfun_0.62        
##  [5] DBI_1.3.0         mitools_2.7       forcats_1.0.1     jsonlite_2.0.0   
##  [9] glue_1.8.1        htmltools_0.5.9   sass_0.4.10       hms_1.1.4        
## [13] rmarkdown_2.32    evaluate_1.0.5    jquerylib_0.1.4   tibble_3.3.1     
## [17] fastmap_1.2.0     yaml_2.3.12       lifecycle_1.0.5   compiler_4.6.2   
## [21] Rcpp_1.1.2        pkgconfig_2.0.3   rstudioapi_0.19.0 lattice_0.23-1   
## [25] digest_0.6.39     R6_2.6.1          pillar_1.11.1     splines_4.6.2    
## [29] magrittr_2.0.5    bslib_0.12.0      tools_4.6.2       cachem_1.1.0

Colorectal Cancer Screening in Arkansas, 2022–2024: Corrected Weighted Analysis

Samad Alvi

2026-11-13