# Usual source of care and hypertension control among adults with hypertension, NHANES 2007-2012.
# Dinkler, Sugar, Escarce, Ong, and Mangione (2016), Am J Hypertens, doi:10.1093/ajh/hpw010.
# Headline: OR 3.89 (95% CI 2.15-6.98) for having a usual source of care, Table 2 Model 1
# (survey-weighted logistic regression without imputation).

# Usual source of care from the place question. 2007-2012 ask HUQ030 (is there a place you usually
# go when sick or for advice: 1 yes, 2 no place, 3 more than one place) and, of those with a
# place, what kind of place they go to most often (HUQ040: 1 clinic or health center, 2 doctor's
# office or HMO, 3 hospital emergency room, 4 hospital outpatient department, 5 some other place).
# 2021-2023 asks HUQ042 instead (1 doctor's office or health center, 2 urgent care center or clinic
# in a drug store or grocery store, 3 emergency room, 4 VA medical center or outpatient clinic,
# 5 some other place, 6 doesn't go to one place most often). The paper counts no place and the
# emergency room as no usual source; every other place counts as one, so urgent care and retail
# clinics (which 2007-2012 respondents would have reported as a clinic or health center), the VA,
# and "no one place most often" (said by those with more than one place) are a usual source.
row111_usual_source <- function(cycle) {
  d <- component("HUQ", cycle)
  place <- if (has(d, "HUQ042")) d$HUQ042 else d$HUQ040
  places <- if (has(d, "HUQ042")) c(1, 2, 4, 5, 6) else c(1, 2, 4, 5)
  has_place <- d$HUQ030 %in% c(1, 3)
  usoc <- ifelse(d$HUQ030 %in% 2 | (has_place & place %in% 3), 0, ifelse(has_place & place %in% places, 1, NA))
  data.frame(SEQN = d$SEQN, usoc = usoc, health = ifelse(d$HUQ010 %in% 1:5, d$HUQ010, NA))
}

# Meeting the American College of Sports Medicine guideline (Garber et al. 2011, the paper's
# reference 15): moderate activity for 30 minutes or more on 5 or more days a week, or vigorous
# activity for 20 minutes or more on 3 or more days, in any of the Global Physical Activity
# Questionnaire's domains (vigorous and moderate work, walking or cycling to get places, vigorous
# and moderate leisure). Built this way it gives Table 1's counts within 1% (2,277 and 425 in its
# population against 2,262 and 418 with and without a usual source; leisure alone gives about
# 1,600 and 225). 2021-2023 asks about leisure time only, so it is NA there.
row111_activity <- function(cycle) {
  d <- component("PAQ", cycle)
  if (!has(d, "PAQ605")) return(data.frame(SEQN = d$SEQN, active = NA_real_))
  domain <- function(answer, days, minutes, min_days, min_minutes) {
    n <- ifelse(days %in% 1:7, days, NA)
    m <- ifelse(minutes %in% c(7777, 9999), NA, minutes)
    ifelse(answer %in% 2, FALSE, ifelse(answer %in% 1, n >= min_days & m >= min_minutes, NA))
  }
  meets <- cbind(domain(d$PAQ605, d$PAQ610, d$PAD615, 3, 20), domain(d$PAQ620, d$PAQ625, d$PAD630, 5, 30),
                 domain(d$PAQ635, d$PAQ640, d$PAD645, 5, 30), domain(d$PAQ650, d$PAQ655, d$PAD660, 3, 20),
                 domain(d$PAQ665, d$PAQ670, d$PAD675, 5, 30))
  # Meets it in any domain: yes; in none, with every domain known: no; otherwise unknown.
  data.frame(SEQN = d$SEQN, active = ifelse(apply(meets, 1, function(r) any(r %in% TRUE)), 1,
                                            ifelse(apply(meets, 1, function(r) all(r %in% FALSE)), 0, NA)))
}

association <- list(
  id = "row111", row = 111, doi = "10.1093/ajh/hpw010",
  cycles = c("2007-2008", "2009-2010", "2011-2012"),
  # Serum creatinine (BIOPRO) is in the model.
  weight = "WTMEC2YR", blood_file = "BIOPRO",
  family = "logistic", term = "usoc",
  published = list(measure = "OR", estimate = 3.89, low = 2.15, high = 6.98, n = NA, events = NA,
                   contrast = "a usual source of care (any place other than the emergency room) vs none or the emergency room"),
  left_out = c(
    "meeting the physical activity guideline over work, transport, and leisure (GPAQ, 2007-2012): 2021-2023 asks about leisure time only",
    "the paper's complete-case drop of adults who never had their cholesterol checked (BPQ060 = 2, who were not asked BPQ080, so hyperlipidemia was missing): 2021-2023 asks everyone whether they were told of high cholesterol and has no BPQ060, so the harmonized version counts the never checked as not told and keeps them"
  ),
  build = function(cycle) {
    demo <- demographics(cycle)
    bpq <- component("BPQ", cycle)
    d <- merge_all(demo, blood_pressure(cycle, readings = 2:4)[, c("SEQN", "sbp", "dbp")], bp_questions(cycle),
                   row111_usual_source(cycle), component("HIQ", cycle, "HIQ011"),
                   medical_conditions(cycle)[, c("SEQN", "heart_failure", "heart_attack", "stroke", "copd")],
                   diabetes_questions(cycle)[, c("SEQN", "told_diabetes")], body_measures(cycle)[, c("SEQN", "bmi")],
                   biochemistry(cycle)[, c("SEQN", "creatinine")], smoking(cycle), row111_activity(cycle))
    # Blood pressure: readings averaged after discarding the first (2-4 by auscultation in
    # 2007-2012, 2-3 oscillometric in 2021-2023), diastolic zeros left out.
    d$controlled <- ifelse(is.na(d$sbp) | is.na(d$dbp), NA, as.integer(d$sbp < 140 & d$dbp < 90))
    d$age_group <- relevel(cut(d$age, c(-Inf, 34, 44, 54, 64, 74, Inf), labels = c("18-34", "35-44", "45-54", "55-64", "65-74", "75+")), ref = "65-74")
    d$race4 <- factor(ifelse(d$race %in% 1:2, "hispanic", ifelse(d$race == 3, "white", ifelse(d$race == 4, "black", "other"))),
                      levels = c("white", "hispanic", "black", "other"))
    d$male <- as.integer(d$sex == 1)
    d$married <- ifelse(is.na(d$marital), NA, as.integer(d$marital == 1))
    d$insured <- yes(d$HIQ011)
    # Did not complete high school; high school graduate or some college; college graduate (Table 1).
    d$education3 <- factor(ifelse(d$education %in% 1:2, "less than high school", ifelse(d$education %in% 3:4, "high school or some college",
                                  ifelse(d$education %in% 5, "college graduate", NA))),
                           levels = c("less than high school", "high school or some college", "college graduate"))
    d$income4 <- cut(d$pir, c(-Inf, 1.5, 2.5, 3.5, Inf), right = FALSE, labels = c("<150%", "150-249%", "250-349%", "350%+"))
    d$health <- factor(d$health, levels = 1:5)
    d$smoking <- factor(d$smoking, levels = 1:3)
    d$diabetes <- d$told_diabetes
    # Told of high cholesterol (BPQ080). 2007-2012 ask it only of adults 20 and older whose
    # cholesterol was ever checked (BPQ060 = 1); the paper's model drops the rest as missing (about
    # 10% missing, its Results say). The harmonized version counts the never checked as not told.
    d$hyperlipidemia <- d$told_cholesterol
    never_checked <- if (has(bpq, "BPQ060")) bpq$SEQN[bpq$BPQ060 %in% 2] else numeric()
    d$hyperlipidemia_h <- ifelse(is.na(d$told_cholesterol) & d$SEQN %in% never_checked, 0, d$told_cholesterol)
    d$bmi_group <- relevel(cut(d$bmi, c(-Inf, 18.5, 25, 30, Inf), right = FALSE, labels = c("underweight", "normal", "overweight", "obese")), ref = "normal")
    # Adults taking prescribed blood pressure medication (BPQ050A; BPQ150 in 2021-2023) or with
    # mean pressure of 140/90 or more.
    d$in_population <- d$age >= 18 & (d$bp_medication %in% 1 | (d$sbp >= 140) %in% TRUE | (d$dbp >= 90) %in% TRUE)
    # As the paper's numbers show it was built (Stata compares a missing value as larger than any
    # number): adults without a pressure reading are counted as hypertensive and as uncontrolled.
    d$in_population_missing_bp <- d$age >= 18 & (d$bp_medication %in% 1 | !(d$sbp < 140) %in% TRUE | !(d$dbp < 90) %in% TRUE)
    d$controlled_missing_bp <- ifelse(is.na(d$controlled), 0L, d$controlled)
    d
  },
  formula = controlled ~ usoc + age_group + race4 + male + married + insured + education3 + income4 + health + active +
    smoking + diabetes + heart_failure + heart_attack + stroke + hyperlipidemia + creatinine + bmi_group,
  formula_harmonized = controlled ~ usoc + age_group + race4 + male + married + insured + education3 + income4 + health +
    smoking + diabetes + heart_failure + heart_attack + stroke + hyperlipidemia_h + creatinine + bmi_group,
  variants = list(
    list(label = "missing pressure counted hypertensive, uncontrolled (Stata)", sample = "own",
         derive = function(data, constants) { data$in_population <- data$in_population_missing_bp; data },
         formula = controlled_missing_bp ~ usoc + age_group + race4 + male + married + insured + education3 + income4 + health + active +
           smoking + diabetes + heart_failure + heart_attack + stroke + hyperlipidemia + creatinine + bmi_group),
    list(label = "Table 2's covariates only", sample = "own",
         formula = controlled ~ usoc + age_group + race4 + male + married + insured + heart_failure + diabetes + hyperlipidemia + bmi_group),
    list(label = "Table 2's covariates only, Stata missing pressure", sample = "own",
         derive = function(data, constants) { data$in_population <- data$in_population_missing_bp; data },
         formula = controlled_missing_bp ~ usoc + age_group + race4 + male + married + insured + heart_failure + diabetes + hyperlipidemia + bmi_group),
    list(label = "unweighted", weighted = FALSE),
    list(label = "interview weight", weight = "WTINT2YR"),
    list(label = "COPD added (Table 1)", formula = controlled ~ usoc + age_group + race4 + male + married + insured + education3 + income4 + health + active +
           smoking + diabetes + heart_failure + heart_attack + stroke + hyperlipidemia + creatinine + bmi_group + copd, sample = "own"),
    list(label = "without physical activity, paper's sample", formula = controlled ~ usoc + age_group + race4 + male + married + insured + education3 + income4 + health +
           smoking + diabetes + heart_failure + heart_attack + stroke + hyperlipidemia + creatinine + bmi_group)
  ),
  # What the coding check found this paper's analysis to have computed otherwise than its text
  # says (the plan's section 6.6).
  departures = list(
    list(kind = "sample", affects_headline = FALSE, followed = NA,
         detail = "The Methods restrict the sample to adults taking blood pressure medication or with a systolic pressure of 140 or diastolic of 90 at the examination, but Table 1's population (n = 7,653, with its counts by age group, usual source, and treatment within 1%) is reproduced only by also counting as hypertensive about 1,270 adults with no pressure reading. Table 1's control shares (54.8% and 17%) and Table 2's odds ratios leave them out (counted as uncontrolled, they give 3.05 for the headline, against 3.41 without them).")
  ),
  choices = list(
    list(choice = "model covariates", decision = "every variable the Measures section defines: age group, race and ethnicity, sex, married, insured, education, income, self-rated health, physical activity, smoking, diabetes, heart failure, prior heart attack, prior stroke, high cholesterol, creatinine, BMI group",
         reason = "Table 2 shows nine covariates and says not all are shown; with all the Measures variables the shown odds ratios match Table 2's (age 18-34 0.39 vs 0.38, over 74 0.60 vs 0.58, 35-44 0.85 vs 0.84, Black 0.69 vs 0.70, male 0.73 vs 0.76, heart failure 1.28 vs 1.30, diabetes 1.04 vs 1.02, high cholesterol 1.53 vs 1.55, obese 2.05 vs 2.03), better than with the shown ones alone (over 74 0.66, 35-44 0.95, Black 0.79, heart failure 1.40, diabetes 1.22)"),
    list(choice = "population and missing blood pressure", decision = "adults 18 and older taking prescribed blood pressure medicine or with mean pressure 140/90 or more; control is undefined without a reading, so those without one leave the model",
         reason = "the paper's Table 1 counts (n = 7,653; by age group, usual source, and treatment within 1%) are reproduced only if adults without a pressure reading (mostly those interviewed but not examined, hence its 10% missing BMI) are counted hypertensive, the way Stata compares missing values; its Table 1 control shares (54.8% and 17%) and Table 2 odds ratios are those without them, so control was missing for them (a variant counts them as uncontrolled)"),
    list(choice = "blood pressure readings", decision = "mean of the readings after the first (2-4; diastolic zeros left out, as the library does)", reason = "the paper averages after discarding the first reading"),
    list(choice = "usual source of care", decision = "no place (HUQ030 = 2) or the emergency room (HUQ040 = 3) is none; any other place is one; a place of unknown kind is missing",
         reason = "the paper's definition; 'some other place' and more than one place count as a usual source, since only none and the emergency room are named as none"),
    list(choice = "2021-2023 place types", decision = "urgent care or retail clinic, VA, some other place, and no one place most often count as a usual source; the emergency room as none",
         reason = "only the emergency room and no place are none in the paper; urgent care and retail clinics would have been reported as a clinic or health center in 2007-2012, and 'no one place most often' is said by those with more than one place"),
    list(choice = "weights", decision = "examination weight over the three cycles (WTMEC2YR/3; the phlebotomy weight in 2021-2023)", reason = "pressure and BMI come from the examination; Table 1's weighted shares (90.7% with a usual source, 54.8% and 17% controlled, 70.7% and 20.1% treated) match the examination weight"),
    list(choice = "categorical coding", decision = "married or living with a partner vs not; insured vs not (HIQ011); diabetes told (borderline as no); income <150%, 150-249%, 250-349%, 350% or more of poverty; smoking never, former, current; self-rated health in its five answers; creatinine continuous",
         reason = "the paper's definitions where it gives them; smoking in Table 1's three groups; Table 1's diabetes counts match DIQ010 = 1"),
    list(choice = "physical activity", decision = "30 minutes on 5 days of moderate or 20 minutes on 3 days of vigorous activity in any GPAQ domain", reason = "the ACSM guideline the paper cites, which gives Table 1's counts"),
    list(choice = "missing covariates", decision = "complete case", reason = "Model 1 is the model without imputation")
  )
)