Executive summary

Purpose (North Star). Give the Board of a regional U.S. health plan the evidence it needs to redesign how members are invited to its cardiovascular prevention program: replace single-criterion invitations (an age cut-off or a single diagnosis) with a combined risk score, so that limited cardiology capacity finds the most cases.

Big idea. With the same number of invitations, a combined risk score finds up to 16 percentage points more heart-disease cases than any single-criterion rule. The Board should approve a 12-month staged pilot to confirm this in the plan’s own population.

What the data show

The burden is concentrated. 8.6% of the 319,795 U.S. adults surveyed report coronary heart disease or a heart attack. Adults aged 65+ are 35% of the population but 67% of all cases.
Single-criterion rules leave cases behind. Inviting everyone aged 65+ finds about 67% of cases; inviting the same number of people ranked by the combined score finds about 83%.
The score works with 11 simple questions. Two models were compared on the same test set; the logistic regression (AUC 0.84) outperformed a more complex random forest (AUC 0.78) and is easier to explain, so it was selected. It detects about 8 in 10 cases while inviting about 3 in 10 adults: roughly 23 cases per 100 invitations, versus 9 at random.
Averages under-plan capacity. Staffing screening days for the average number of cases leaves cardiology slots short on about 35% of days; planning at the 95th percentile (13 slots per 100 patients) covers 95% of days.
Statistical significance is not managerial relevance. With 320,000 records every test is “significant”; effect sizes show which differences actually matter.

Decision requested from the Board

Approve a 12-month, three-stage pilot of score-based invitations with go/no-go checkpoints at months 3, 6 and 12. The pilot stops at month 3 if, on the plan’s own data, the score does not beat the current rule by at least 5 percentage points of cases found with the same number of invitations.
Adopt percentile-based capacity planning for cardiology screening days.

What this report does not claim. The data are a national survey, not the plan’s members; the outcome is heart disease already diagnosed, while screening seeks undiagnosed people; and there are no cost data, so the report speaks of efficiency (cases per invitation), not savings. The pilot is designed to resolve exactly these gaps.

Beyond the pilot. The report closes with a roadmap to transfer the method (not the numbers) to Colombia’s EPS system, including the data it would use, how the outcome would be improved and the errors to expect.

The interactive Board dashboard that monitors these KPIs accompanies this report: open the Board dashboard.

Strategic framing: the A³ model

The report follows the A³ model (Attention → Analysis → Action): it first captures attention with a management problem, then builds the evidence, and closes with concrete actions.

Attention: the management problem

Coronary heart disease and heart attack are among the leading causes of death in the United States and among the most resource-intensive conditions: specialist visits, hospital stays, intensive care and catheterization procedures. A health plan cannot send all its adult members to cardiology screening at once, so programs select members one criterion at a time: an age cut-off, or a single diagnosis such as diabetes. The management question is whether combining risk factors into one score finds more cases with the same capacity, and how much capacity to plan.

Audience

The Board of Directors decides; the managers propose and execute. The Board’s main concern is a program that falls short of its targets, which is why the proposal is a staged pilot with go/no-go checkpoints.

Audience	What they care about	What they must decide or do
Board of Directors, regional U.S. health plan	Health outcomes, cost, meeting targets	Approve the staged pilot
VP of Population Health / Care Management	Finding high-risk members early	Propose and run the score-based invitations
VP of Operations / Capacity Planning	Slots, staff and beds that match demand	Staff screening days
Chief Medical Officer	Clinical relevance, no missed patients	Validate criteria and follow-up routes

North Star and objectives

General objective: use the BRFSS 2020 heart disease data to test whether a combined risk score should replace single-criterion invitations in a cardiovascular prevention program, and to size the capacity it requires.

Chapter	Specific objective (Board question)	Methods
1. Exploratory Data Analysis	Where is the burden? Describe the population and locate where heart disease concentrates	Data quality, descriptive indicators, exploratory charts
2. Probability & Applied Statistics	How much capacity do we need, and which differences are real?	Empirical and conditional probability, simulation, hypothesis tests, effect sizes
3. Analytical Models	Can we identify high-risk people with simple questions?	Two models (logistic regression vs random forest), cross-validation, diagnostics, AUC and sensitivity
4. Data Storytelling	What should the Board decide, and how will it be monitored? Can the method travel?	Annotated charts, interactive dashboard, executive recommendations, transfer roadmap to Colombia

Success KPIs for the program: the score beats the current rule by ≥ 5 percentage points of cases found with the same invitations; AUC ≥ 0.75 on the plan’s own data; sensitivity ≥ 70%; ≤ 5% of screening days with unmet demand.

# Analytical roadmap: how each chapter feeds the final decision
DiagrammeR::grViz("
digraph roadmap {
  graph [rankdir = LR, fontname = Helvetica]
  node  [shape = box, style = 'rounded,filled', fillcolor = '#E8EEF7',
         fontname = Helvetica, fontsize = 11]

  north [label = 'NORTH STAR\\nCombined risk score instead of\\nsingle-criterion invitations', fillcolor = '#1F78B4', fontcolor = white]
  c1    [label = '1. Exploration\\nWhere is the burden?']
  c2    [label = '2. Probability & statistics\\nHow much capacity?\\nWhich differences matter?']
  c3    [label = '3. Analytical model\\nWho is at high risk?']
  c4    [label = '4. Storytelling & dashboard\\nWhat should the Board decide?']
  dec   [label = 'BOARD DECISION\\n12-month staged pilot +\\npercentile capacity plan', fillcolor = '#FBE3D0']

  north -> c1 -> c2 -> c3 -> c4 -> dec
}
")

1 Exploratory Data Analysis

1.1 The data

Source: Personal Key Indicators of Heart Disease (Kaggle), built from the CDC’s 2020 Behavioral Risk Factor Surveillance System (BRFSS), a telephone survey of U.S. adults. Unit of analysis: one row = one surveyed adult.

The dataset has 319,795 adults and 18 variables. The outcome is HeartDisease: whether the person reports ever having coronary heart disease or a heart attack. The predictors are four numeric variables (BMI, days of poor physical and mental health in the last month, sleep hours) and thirteen categorical variables (smoking, heavy drinking, stroke, difficulty walking, sex, age group, race, diabetes, physical activity, self-rated general health, asthma, kidney disease and skin cancer).

1.1.1 Structure and data quality

# Structure: variable names, types and first values
glimpse(heart)

## Rows: 319,795
## Columns: 18
## $ HeartDisease     <chr> "No", "No", "No", "No", "No", "Yes", "No", "No", "No"…
## $ BMI              <dbl> 16.60, 20.34, 26.58, 24.21, 23.71, 28.87, 21.63, 31.6…
## $ Smoking          <chr> "Yes", "No", "Yes", "No", "No", "Yes", "No", "Yes", "…
## $ AlcoholDrinking  <chr> "No", "No", "No", "No", "No", "No", "No", "No", "No",…
## $ Stroke           <chr> "No", "Yes", "No", "No", "No", "No", "No", "No", "No"…
## $ PhysicalHealth   <dbl> 3, 0, 20, 0, 28, 6, 15, 5, 0, 0, 30, 0, 0, 7, 0, 1, 5…
## $ MentalHealth     <dbl> 30, 0, 30, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 30, 0, 2,…
## $ DiffWalking      <chr> "No", "No", "No", "No", "Yes", "Yes", "No", "Yes", "N…
## $ Sex              <chr> "Female", "Female", "Male", "Female", "Female", "Fema…
## $ AgeCategory      <chr> "55-59", "80 or older", "65-69", "75-79", "40-44", "7…
## $ Race             <chr> "White", "White", "White", "White", "White", "Black",…
## $ Diabetic         <chr> "Yes", "No", "Yes", "No", "No", "No", "No", "Yes", "N…
## $ PhysicalActivity <chr> "Yes", "Yes", "Yes", "No", "Yes", "No", "Yes", "No", …
## $ GenHealth        <chr> "Very good", "Very good", "Fair", "Good", "Very good"…
## $ SleepTime        <dbl> 5, 7, 8, 6, 8, 12, 4, 9, 5, 10, 15, 5, 8, 7, 5, 6, 10…
## $ Asthma           <chr> "Yes", "No", "Yes", "No", "No", "No", "Yes", "Yes", "…
## $ KidneyDisease    <chr> "No", "No", "No", "No", "No", "No", "No", "No", "Yes"…
## $ SkinCancer       <chr> "Yes", "No", "No", "Yes", "No", "No", "Yes", "No", "N…

# Quality check before any analysis: missing values and identical rows
c(missing_values = sum(is.na(heart)),
  duplicated_rows = sum(duplicated(heart)))

##  missing_values duplicated_rows 
##               0           18078

# Outcome distribution: the classes are strongly imbalanced
count(heart, HeartDisease)

Interpretation. There are no missing values. There are 18,078 identical rows (5.7%); since the data have no person identifier and most questions are yes/no, different people can give exactly the same answers, so they were kept. Sleep hours (1 to 24) and BMI (12 to 94.9) include some implausible extreme values; they were kept but are considered when interpreting. Only 8.6% of adults have the outcome: the classes are imbalanced, which shapes every later decision (for example, why accuracy is not used to judge the model).

1.2 Descriptive indicators

Management situation: before comparing groups, management needs the typical values and the variability of the key health indicators.

# BMI indicators with base R functions
c(mean     = mean(heart$BMI),
  median   = median(heart$BMI),
  sd       = sd(heart$BMI),
  variance = var(heart$BMI))

##     mean   median       sd variance 
## 28.32540 27.34000  6.35610 40.40001

range(heart$BMI)

## [1] 12.02 94.85

# Full summary of the four numeric variables
heart |>
  select(BMI, PhysicalHealth, MentalHealth, SleepTime) |>
  summary()

##       BMI        PhysicalHealth    MentalHealth      SleepTime     
##  Min.   :12.02   Min.   : 0.000   Min.   : 0.000   Min.   : 1.000  
##  1st Qu.:24.03   1st Qu.: 0.000   1st Qu.: 0.000   1st Qu.: 6.000  
##  Median :27.34   Median : 0.000   Median : 0.000   Median : 7.000  
##  Mean   :28.33   Mean   : 3.372   Mean   : 3.898   Mean   : 7.097  
##  3rd Qu.:31.42   3rd Qu.: 2.000   3rd Qu.: 3.000   3rd Qu.: 8.000  
##  Max.   :94.85   Max.   :30.000   Max.   :30.000   Max.   :24.000

# The same indicators, split by heart disease status
heart |>
  group_by(HeartDisease) |>
  summarise(
    bmi_mean             = mean(BMI),
    bmi_median           = median(BMI),
    sleep_mean           = mean(SleepTime),
    poor_physical_days   = mean(PhysicalHealth),
    poor_mental_days     = mean(MentalHealth)
  )

1.2.1 Exploratory chart

# Share of adults with heart disease in each age group
prevalence_age <- heart |>
  group_by(AgeCategory) |>
  summarise(prevalence = mean(HeartDisease == "Yes") * 100)

plot_age <- ggplot(prevalence_age, aes(x = AgeCategory, y = prevalence)) +
  geom_col(fill = "#3B528B") +
  labs(
    title = "Share of adults with heart disease by age group",
    x = "Age group (years)",
    y = "% with heart disease"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplotly(plot_age)

1.2.2 Interpretation

BMI. Mean BMI is 28.3 kg/m², with a median of 27.3, a standard deviation of 6.4 (variance 40.4) and a range of 12.0 to 94.9. The mean above the median indicates a right-skewed distribution, and the average adult is already in the overweight range (25–29.9).

Poor-health days and sleep. For poor physical and mental health days the median is 0: most adults report none, but a group reports many (up to 30), which pulls the means up (3.4 and 3.9 days). In such variables the median describes the typical person better than the mean. Average sleep is 7.1 hours.

By group. Adults with heart disease report 7.8 days of poor physical health per month versus 3.0 for the rest (2.6 times more) and a higher mean BMI (29.4 vs 28.2). Differences in mental health (4.6 vs 3.8 days) and sleep (7.1 hours in both groups) are small. The exploratory chart shows prevalence rising steeply with age, from 0.6% to 22.6%.

For management: poor physical health days are a simple proxy for disease burden: people with heart disease lose almost 5 more days per month, which can translate into absenteeism and higher service use. Limitation: all indicators are self-reported.

2 Probability and Applied Statistics

2.1 Applied probability

2.1.1 Random selection and sample space

Management situation: a cardiovascular screening program receives 100 patients. Each one may or may not have heart disease (sample space: {Yes, No}). Using the probability observed in the survey, we simulate how many cases would arrive.

# The seed makes randomness reproducible: anyone who runs the report gets the same result
set.seed(2026)

sample_space <- c("Yes", "No")

# Observed probability in the data: share of "Yes" answers
p_heart <- mean(heart$HeartDisease == "Yes")

screening <- sample(
  sample_space,
  size = 100,
  replace = TRUE,                     # each patient is a new draw
  prob = c(p_heart, 1 - p_heart)      # same order as sample_space
)

table(screening)

## screening
##  No Yes 
##  91   9

Audit sample from the real data: to audit the information without reviewing all 319,795 records, 1,000 adults are selected at random.

set.seed(2026)

# replace = FALSE (the default): nobody can be selected twice
selected_rows <- sample(nrow(heart), size = 1000)
audit_sample <- heart |> slice(selected_rows)

# If both proportions are close, the sample represents the full data well
c(full_data  = p_heart,
  sample_1000 = mean(audit_sample$HeartDisease == "Yes"))

##   full_data sample_1000 
##  0.08559546  0.09200000

Interpretation

Simulated screening. With the observed probability (8.6%), the simulation of 100 screened patients produced 9 cases, close to the expected 8.6 per 100. The exact number changes in each simulation because of chance, so screening capacity (for example, cardiology slots) should be planned with a range, not a single number.

Audit sample. In the random sample of 1,000 adults, 9.2% have heart disease, versus 8.6% in the full data. The small gap is natural sampling variation: a well-drawn random sample reproduces the population, which lets management audit or estimate indicators faster and at lower cost.

2.1.2 Empirical and conditional probabilities

Management situation: to prioritize prevention programs, management needs to know in which groups heart disease concentrates.

# Simple probability: cases / total
freq <- table(heart$HeartDisease)
freq

## 
##     No    Yes 
## 292422  27373

prop.table(freq)

## 
##         No        Yes 
## 0.91440454 0.08559546

# margin = 1 computes the % within each ROW, i.e., within each group:
# "among smokers, what % has heart disease?"
tab_smoking <- table(Smoker = heart$Smoking, HeartDisease = heart$HeartDisease)
round(prop.table(tab_smoking, margin = 1) * 100, 1)

##       HeartDisease
## Smoker   No  Yes
##    No  94.0  6.0
##    Yes 87.8 12.2

tab_stroke <- table(Stroke = heart$Stroke, HeartDisease = heart$HeartDisease)
round(prop.table(tab_stroke, margin = 1) * 100, 1)

##       HeartDisease
## Stroke   No  Yes
##    No  92.5  7.5
##    Yes 63.6 36.4

tab_age <- table(Age = heart$AgeCategory, HeartDisease = heart$HeartDisease)
round(prop.table(tab_age, margin = 1) * 100, 1)

##              HeartDisease
## Age             No  Yes
##   18-24       99.4  0.6
##   25-29       99.2  0.8
##   30-34       98.8  1.2
##   35-39       98.6  1.4
##   40-44       97.7  2.3
##   45-49       96.6  3.4
##   50-54       94.6  5.4
##   55-59       92.6  7.4
##   60-64       90.1  9.9
##   65-69       88.0 12.0
##   70-74       84.4 15.6
##   75-79       81.2 18.8
##   80 or older 77.4 22.6

# High risk does not always mean many cases: compare each group's weight
# in the population and in the total number of cases
age_65_plus <- c("65-69", "70-74", "75-79", "80 or older")
cases <- heart |> filter(HeartDisease == "Yes")

c(
  population_65_plus = mean(heart$AgeCategory %in% age_65_plus),
  cases_65_plus      = mean(cases$AgeCategory %in% age_65_plus),
  cases_with_stroke  = mean(cases$Stroke == "Yes")
)

## population_65_plus      cases_65_plus  cases_with_stroke 
##          0.3466314          0.6738757          0.1603405

Interpretation

Overall probability. 27,373 of 319,795 adults report heart disease: an empirical probability of 0.086 (8.6 per 100).

Conditional probabilities. Among smokers 12.2% have heart disease versus 6.0% of non-smokers: double (relative risk ≈ 2.0). Among adults with a prior stroke, 36.4% versus 7.5%: almost 5 times more. With age, the probability rises from 0.6% (18–24) to 22.6% (80+), accelerating after 50. Because the probability differs from the overall 8.6% in each group, these events are dependent: knowing whether someone smokes, had a stroke or is older changes their probability of heart disease.

For management. With budget for one program, target adults aged 65+: they are 34.7% of the population but 67.4% of all cases. Stroke survivors have the highest risk (36.4%) but contribute only 16% of cases, so they fit a specific follow-up route linked to stroke care. Limitation: cross-sectional, self-reported data show association, not causation.

2.1.3 Discrete distributions: capacity scenarios

Management situation: management must size three services: (1) a prevention program for adults aged 80+, (2) screening days of 100 patients and (3) emergency beds for heart-attack patients.

set.seed(2026)

# Probability among adults aged 80+, taken from the age table above
p_age_80 <- 0.226

# A random number between 0 and 1 falls below 0.226 in 22.6% of cases:
# this simulates whether each patient has heart disease
random_numbers <- runif(100)
cases_age_80 <- sum(random_numbers < p_age_80)

c(general_program = sum(screening == "Yes"),
  program_80_plus = cases_age_80)

## general_program program_80_plus 
##               9              24

set.seed(2026)

# rbinom() repeats the 100-patient screening on 1,000 days and counts cases per day
cases_per_day <- rbinom(n = 1000, size = 100, prob = p_heart)

summary(cases_per_day)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   7.000   8.000   8.504  10.000  18.000

# Range covering the central 90% of days: the basis for planning slots
quantile(cases_per_day, probs = c(0.05, 0.95))

##  5% 95% 
##   4  13

# Share of days on which 9 slots (the average) would not be enough
mean(cases_per_day > 9)

## [1] 0.353

plot_days <- tibble(cases = cases_per_day) |>
  ggplot(aes(x = cases)) +
  geom_bar(fill = "#3B528B") +
  labs(
    title = "Cases per day across 1,000 simulated screening days",
    x = "Patients with heart disease (out of 100 screened)",
    y = "Number of days"
  ) +
  theme_minimal()

# ggplotly() makes the chart interactive: hovering shows exact values
ggplotly(plot_days)

# ASSUMPTION for the example (not from the data): an emergency unit receives
# on average 3 heart-attack patients per day
lambda_mi <- 3

prob_admissions <- tibble(
  admissions = 0:10,
  probability = dpois(admissions, lambda = lambda_mi)
)
prob_admissions

# With 5 beds for new admissions, how likely is it that they are not enough?
1 - sum(dpois(0:5, lambda = lambda_mi))

## [1] 0.08391794

# And with a sixth bed?
1 - sum(dpois(0:6, lambda = lambda_mi))

## [1] 0.03350854

plot_admissions <- prob_admissions |>
  ggplot(aes(x = admissions, y = probability)) +
  geom_col(fill = "#3B528B") +
  labs(
    title = "Probability of daily heart-attack admissions (assumed mean: 3 per day)",
    x = "Heart-attack patients admitted in a day",
    y = "Probability"
  ) +
  theme_minimal()

ggplotly(plot_admissions)

Interpretation

Program for adults 80+ (runif()). With the probability observed in that group (22.6%), 24 of 100 enrolled patients had heart disease, versus 9 in the general program: almost 3 times more cases for the same enrollment, which requires more cardiology hours, follow-up and diagnostics per 100 patients.

Screening days (rbinom()). Across 1,000 simulated days the mean was 8.5 cases per day, ranging from 1 to 18. On 9 of 10 days there were between 4 and 13 cases. If only 9 slots (the average) were booked, they would fall short on 35.3% of days; booking 13 slots covers about 95% of them.

Emergency beds (dpois()). Under the assumption of 3 heart-attack admissions per day, demand exceeds 5 beds on 8.4% of days (about 31 days a year); a sixth bed lowers this to 3.4% (about 12 days a year). This allows a data-based discussion of whether the cost of an extra bed is justified by the saturation days it avoids.

Limitation: simulations depend on their assumptions; the mean of 3 admissions per day is assumed, not observed, and the Poisson model assumes independent arrivals at a constant rate.

2.1.4 Continuous distributions: BMI, arrival times and visit length

Management situation: management needs to (1) estimate obesity to size weight-control programs, (2) know how often heart-attack patients arrive back-to-back to organize the catheterization lab and (3) compute the professional hours a screening day requires.

# Real BMI parameters
bmi_mean <- mean(heart$BMI)
bmi_sd   <- sd(heart$BMI)

set.seed(2026)
# If BMI followed a normal distribution with these parameters...
bmi_simulated <- rnorm(1000, mean = bmi_mean, sd = bmi_sd)

# Compare simulated and real BMI: do the extremes match?
summary(bmi_simulated)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   8.952  24.019  28.419  28.413  32.283  50.041

summary(heart$BMI)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   12.02   24.03   27.34   28.33   31.42   94.85

# Obesity (BMI >= 30): pnorm() gives P(BMI < 30), so it is subtracted from 1
c(normal_model = 1 - pnorm(30, mean = bmi_mean, sd = bmi_sd),
  observed     = mean(heart$BMI >= 30))

## normal_model     observed 
##    0.3960966    0.3215873

# Severe obesity (BMI >= 40)
c(normal_model = 1 - pnorm(40, mean = bmi_mean, sd = bmi_sd),
  observed     = mean(heart$BMI >= 40))

## normal_model     observed 
##   0.03312301   0.05123595

# Theoretical normal curve computed with dnorm(), overlaid on real data
normal_curve <- tibble(
  bmi = seq(10, 70, by = 0.5),
  density = dnorm(bmi, mean = bmi_mean, sd = bmi_sd)
)

plot_bmi <- ggplot() +
  geom_histogram(
    data = audit_sample,
    aes(x = BMI, y = after_stat(density)),
    binwidth = 2, boundary = 0, fill = "#3B528B", alpha = 0.6
  ) +
  geom_line(data = normal_curve, aes(x = bmi, y = density),
            color = "#D95F02", linewidth = 1) +
  labs(
    title = "Real BMI (sample of 1,000) versus the theoretical normal curve",
    x = "Body mass index (kg/m²)",
    y = "Density"
  ) +
  theme_minimal()

ggplotly(plot_bmi)

# Same ASSUMPTION as above: 3 heart attacks per day = 3/24 per hour.
# If arrivals follow a Poisson process, the time between them is exponential
rate_per_hour <- 3 / 24

# Average time between arrivals, in hours
1 / rate_per_hour

## [1] 8

set.seed(2026)
interarrival_times <- rexp(1000, rate = rate_per_hour)
summary(interarrival_times)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##  0.004507  2.277612  5.849029  8.325302 11.943313 55.446716

# If the catheterization lab needs 4 hours to be ready again,
# how likely is it that another heart attack arrives before then?
pexp(4, rate = rate_per_hour)

## [1] 0.3934693

exp_curve <- tibble(
  hours = seq(0, 40, by = 0.5),
  density = dexp(hours, rate = rate_per_hour)
)

plot_arrivals <- ggplot() +
  geom_histogram(
    data = tibble(hours = interarrival_times),
    aes(x = hours, y = after_stat(density)),
    binwidth = 2, boundary = 0, fill = "#3B528B", alpha = 0.6
  ) +
  geom_line(data = exp_curve, aes(x = hours, y = density),
            color = "#D95F02", linewidth = 1) +
  labs(
    title = "Time between heart-attack arrivals (simulated; assumption: 3 per day)",
    x = "Hours between one arrival and the next",
    y = "Density"
  ) +
  theme_minimal()

ggplotly(plot_arrivals)

set.seed(2026)
# ASSUMPTION: each screening visit lasts 15–25 minutes, any value equally likely
# (continuous uniform distribution)
visit_minutes <- runif(100, min = 15, max = 25)
summary(visit_minutes)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   15.01   17.32   19.60   19.84   22.22   24.78

# Professional hours required by a 100-patient screening day
sum(visit_minutes) / 60

## [1] 33.06733

Interpretation

BMI and the normal distribution (rnorm(), dnorm(), pnorm()). Using the real mean (28.3) and standard deviation (6.4), the normal model predicts 39.6% obesity (BMI ≥ 30) when the observed figure is 32.2%, and 3.3% severe obesity (BMI ≥ 40) when the observed figure is 5.1%. The normal model overestimates obesity and underestimates severe obesity because real BMI is right-skewed: a long tail of very high values that a symmetric curve cannot reproduce. A bariatric program planned on the normal curve would underestimate demand by about a third: proportions should be computed from real data, not from a model that does not fit them.

Catheterization lab arrivals (rexp(), dexp(), pexp()). Under 3 heart attacks per day, the mean time between arrivals is 8 hours, yet in 39.3% of cases the next patient arrives in under 4 hours, the assumed time to reset the lab. In about 4 of 10 arrivals the lab could be busy, which justifies a backup team or a transfer protocol.

Visit length (runif()). If each visit lasts 15–25 minutes, a 100-patient day requires 33.1 professional hours. With 8-hour shifts, 5 professionals are needed (always rounded up: a fraction of a shift still has to be covered).

Limitation: BMI figures come from the data; 3 admissions per day, 4 hours of lab reset, visit length and 8-hour shifts are assumptions of the exercise.

2.2 Inferential statistics

The hypothesis tests follow the same logic as the executive memorandum case discussed in class: a null hypothesis (H0: “no difference” or “no association”) is tested against an alternative (H1), and H0 is rejected when the p-value is below α = 0.05.

Warning: with 319,795 records almost any difference, however small, is statistically significant. That is why confidence intervals and effect sizes (measures of how large a difference is) are reported: they show whether a result matters for management.

2.2.1 t-test: BMI by heart disease status

# Welch t-test: compares two means without assuming equal variances
# H0: mean BMI is equal in both groups; H1: it differs
# (R reports the interval for No - Yes, so it is negative)
t_bmi <- t.test(BMI ~ HeartDisease, data = heart)
t_bmi

## 
##  Welch Two Sample t-test
## 
## data:  BMI by HeartDisease
## t = -28.402, df = 32295, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group No and group Yes is not equal to 0
## 95 percent confidence interval:
##  -1.258156 -1.095712
## sample estimates:
##  mean in group No mean in group Yes 
##          28.22466          29.40159

# Effect size (approximate Cohen's d): the difference in standard deviations;
# below 0.2 is considered small
diff(t_bmi$estimate) / sd(heart$BMI)

## mean in group Yes 
##          0.185166

Interpretation

Mean BMI is 29.40 among adults with heart disease and 28.22 among the rest: a difference of 1.18 kg/m² (95% CI: 1.10 to 1.26). The Welch t-test (t = −28.4; p < 0.001) rejects H0, but the effect size is small (d ≈ 0.19). For management: excess weight is associated with heart disease, but BMI alone does not identify who has it; it must be combined with age, smoking and stroke history to target programs.

2.2.2 One-sample t-test: sleep hours

# Same reasoning as the memorandum (one-sided, lower tail):
# H0: mean sleep >= 7 hours (recommended minimum for adults); H1: mean < 7
t.test(heart$SleepTime, mu = 7, alternative = "less")

## 
##  One Sample t-test
## 
## data:  heart$SleepTime
## t = 38.228, df = 319794, p-value = 1
## alternative hypothesis: true mean is less than 7
## 95 percent confidence interval:
##      -Inf 7.101252
## sample estimates:
## mean of x 
##  7.097075

# Complementary indicator: share sleeping less than 7 hours
mean(heart$SleepTime < 7)

## [1] 0.3032755

Interpretation

Mean sleep is 7.10 hours (t = 38.2; p = 1; upper 95% bound: 7.10 hours), so H0 is not rejected: there is no evidence that the average is below 7 hours. However, 30.3% of adults sleep less than 7 hours. As in the memorandum case, a decision should not rest on a simple reading of the mean: here the mean “complies” but hides that 3 in 10 people sleep too little. The useful indicator for a sleep-hygiene program is the share with insufficient sleep.

2.2.3 ANOVA: BMI by self-rated health

# Order self-rated health from worst to best; otherwise R sorts it alphabetically
heart <- heart |>
  mutate(GenHealth = factor(GenHealth,
                            levels = c("Poor", "Fair", "Good", "Very good", "Excellent")))

# Mean BMI and number of people per group
heart |>
  group_by(GenHealth) |>
  summarise(bmi_mean = mean(BMI), people = n())

# H0: mean BMI is equal in the 5 groups; H1: at least one differs
anova_bmi <- aov(BMI ~ GenHealth, data = heart)
summary(anova_bmi)

##                 Df   Sum Sq Mean Sq F value Pr(>F)    
## GenHealth        4   730416  182604    4791 <2e-16 ***
## Residuals   319790 12189265      38                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Effect size (eta squared): share of BMI variation explained by self-rated health
sum_sq <- summary(anova_bmi)[[1]][["Sum Sq"]]
sum_sq[1] / sum(sum_sq)

## [1] 0.05653513

# ANOVA says there are differences but not between which groups:
# Tukey compares all pairs and adjusts for multiple comparisons
TukeyHSD(anova_bmi)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = BMI ~ GenHealth, data = heart)
## 
## $GenHealth
##                           diff         lwr        upr     p adj
## Fair-Poor            0.2445511  0.06206293  0.4270392 0.0023848
## Good-Poor           -0.9061785 -1.07401337 -0.7383436 0.0000000
## Very good-Poor      -2.6405308 -2.80670576 -2.4743559 0.0000000
## Excellent-Poor      -4.3023031 -4.47366881 -4.1309373 0.0000000
## Good-Fair           -1.1507295 -1.25667374 -1.0447854 0.0000000
## Very good-Fair      -2.8850819 -2.98837627 -2.7817875 0.0000000
## Excellent-Fair      -4.5468541 -4.65830744 -4.4354008 0.0000000
## Very good-Good      -1.7343524 -1.80875906 -1.6599457 0.0000000
## Excellent-Good      -3.3961246 -3.48149712 -3.3107520 0.0000000
## Excellent-Very good -1.6617722 -1.74383333 -1.5797111 0.0000000

Interpretation

Mean BMI changes with self-rated health: 30.4 (poor), 30.6 (fair), 29.5 (good), 27.8 (very good) and 26.1 (excellent). The ANOVA (F = 4,790.7; p < 0.001) rejects H0, and Tukey shows all pairs differ; the largest gap is between fair and excellent (4.5 kg/m²) and the smallest between poor and fair (0.25 kg/m², not relevant in practice). The effect size (η² = 0.057) means self-rated health explains about 6% of BMI variation. For management: the self-rated health question is cheap and easy to ask and can serve as a first filter for weight-control programs.

2.2.4 Chi-square: smoking and heart disease

# H0: smoking and heart disease are independent; H1: they are associated
chi_smoking <- chisq.test(tab_smoking)
chi_smoking

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tab_smoking
## X-squared = 3713, df = 1, p-value < 2.2e-16

# Counts expected if both variables were independent
round(chi_smoking$expected)

##       HeartDisease
## Smoker     No   Yes
##    No  171805 16082
##    Yes 120617 11291

# Strength of association (Cramér's V): 0 = none, 1 = perfect
sqrt(chi_smoking$statistic / sum(tab_smoking))

## X-squared 
## 0.1077528

Interpretation

The chi-square test (X² = 3,713; df = 1; p < 0.001) rejects independence. If both variables were independent we would expect about 11,291 smokers with heart disease; we observe 16,037. Cramér’s V (0.11) indicates a statistically weak association, even though smokers have double the probability (12.2% vs 6.0%). For management: smoking-cessation programs are relevant but not sufficient: 41% of cases occur in non-smokers.

2.2.5 Test of proportions: women and men

# H0: the proportion with heart disease is equal in women and men
tab_sex <- table(Sex = heart$Sex, HeartDisease = heart$HeartDisease)
tab_sex

##         HeartDisease
## Sex          No    Yes
##   Female 156571  11234
##   Male   135851  16139

prop.test(x = tab_sex[, "Yes"], n = rowSums(tab_sex))

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  tab_sex[, "Yes"] out of rowSums(tab_sex)
## X-squared = 1568.3, df = 1, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.04120085 -0.03727487
## sample estimates:
##     prop 1     prop 2 
## 0.06694675 0.10618462

Interpretation

6.7% of women (11,234 of 167,805) and 10.6% of men (16,139 of 151,990) report heart disease (X² = 1,568.3; p < 0.001). With 95% confidence, the proportion in men exceeds that in women by 3.7 to 4.1 percentage points: men have about 1.6 times the probability. For management: screening can prioritize men without neglecting women, who account for 41% of cases. Limitation: not adjusted for age or other factors.

2.2.6 Comparative charts with ggpubr

ggboxplot(
  heart, x = "HeartDisease", y = "BMI",
  fill = "HeartDisease", palette = c("#3B528B", "#D95F02"),
  title = "BMI by heart disease status",
  xlab = "Heart disease",
  ylab = "Body mass index (kg/m²)",
  legend = "none"
) +
  stat_compare_means(method = "t.test")

ggboxplot(
  heart, x = "GenHealth", y = "BMI",
  fill = "#9ECAE1",
  title = "BMI by self-rated health",
  xlab = "Self-rated health",
  ylab = "Body mass index (kg/m²)"
) +
  stat_compare_means(method = "anova")

Interpretation

The BMI boxes of both groups overlap almost completely, although the median is higher with heart disease (28.3 vs 27.3): the test is significant, but the overlap shows visually that the difference is small (statistical significance is not practical relevance). The second chart shows a stepwise decrease in BMI from fair/poor to excellent health, consistent with the ANOVA.

3 Analytical Models: who is at high risk?

3.1 Problem, objectives and success metrics

A health plan cannot send all its adult members to cardiology at once; it needs to prioritize who to invite first. The previous chapters showed that heart disease concentrates in older adults, stroke survivors, smokers and men. This chapter combines all these factors into a classification model that assigns each person a probability of heart disease and labels them as high or low risk. Two analytical models are built and compared on the same test set: a logistic regression (model 1, interpretable) and a random forest (model 2, more flexible). The one that best balances performance and explainability is selected.

Metric	What it measures	Target
AUC (area under the ROC curve)	Ability to rank a person with the disease above one without it; 0.5 = chance, 1 = perfect	≥ 0.80 on the test set
Sensitivity	Out of 100 people with heart disease, how many the model flags as high risk	≥ 70%
Train–test difference	Whether the model “memorized” the data (overfitting)	AUC difference ≤ 0.02

Accuracy is not a primary metric: with only 8.6% positives, a model that says nobody has the disease would be more than 91% accurate while detecting no patients.

3.2 Variable selection

Group	Variables included	Rationale
Demographic	Age (`AgeCategory`), sex (`Sex`)	Prevalence rose from 0.6% to 22.6% with age and was higher in men
Habits and weight	Smoking, physical activity, BMI	Modifiable risk factors targeted by prevention programs
Clinical history	Stroke, diabetes, kidney disease	Conditions sharing mechanisms with heart disease; stroke showed the highest risk
Functional status	Self-rated health, difficulty walking, poor physical health days	Summarize disease burden and are easy to ask

Excluded: race, because a model that allocates resources should not use race as a criterion (equity); heavy drinking, because in these data it is associated with less heart disease (5.2% vs 8.8%), probably because sick people stop drinking, which would send a misleading message; asthma, skin cancer, mental health and sleep, which have a less direct clinical link. A parsimonious model (few variables, easy to apply and explain) was preferred.

Working hypotheses: (1) older age increases the probability; (2) stroke, diabetes and kidney disease increase it; (3) smoking increases it and physical activity decreases it; (4) worse self-rated health is associated with higher probability.

3.3 Conceptual model

# Diagram of the relationships between variables and the chosen technique
DiagrammeR::grViz("
digraph model {
  graph [rankdir = LR, fontname = Helvetica]
  node  [shape = box, style = 'rounded,filled', fillcolor = '#E8EEF7',
         fontname = Helvetica, fontsize = 11]

  demo     [label = 'Demographic\\nAge, sex']
  habits   [label = 'Habits and weight\\nSmoking, physical activity, BMI']
  clinical [label = 'Clinical history\\nStroke, diabetes, kidney disease']
  funct    [label = 'Functional status\\nSelf-rated health, difficulty walking,\\npoor physical health days']

  model    [label = 'Model 1: logistic regression\\n(champion; training: 75%)', fillcolor = '#3B528B', fontcolor = white]
  rf       [label = 'Model 2: random forest\\n(challenger)', fillcolor = '#B8C4DE']
  prob     [label = 'Probability of\\nheart disease']
  decision [label = 'High risk: prioritize screening\\nLow risk: routine follow-up', fillcolor = '#FBE3D0']
  evalu    [label = 'Evaluation (test: 25%)\\nAUC and sensitivity', shape = note]

  {demo habits clinical funct} -> model -> prob -> decision
  {demo habits clinical funct} -> rf [style = dashed]
  rf -> evalu [style = dashed, label = 'compared on\\nthe same test set', fontsize = 9]
  prob -> evalu [style = dashed]
}
")

Technique (model 1): logistic regression, the standard model for a yes/no outcome (model 2, the random forest, is described where it is fitted):

log( p / (1 - p) ) = b0 + b1·age + b2·sex + b3·smoking + ... + b11·poor physical health days

p          = probability of heart disease
p / (1-p)  = odds (e.g., p = 0.20 means odds of 1 to 4)
exp(b)     = odds ratio (OR): how many times the odds are multiplied when the
             variable increases by one unit or changes from "No" to "Yes"

3.4 Implementation in RStudio

3.4.1 Data preparation

model_data <- heart |>
  # Only the variables justified above
  select(HeartDisease, AgeCategory, Sex, Smoking, Stroke, Diabetic, KidneyDisease,
         DiffWalking, PhysicalActivity, GenHealth, BMI, PhysicalHealth) |>
  # No missing values here, but the step keeps the workflow valid for other data
  drop_na() |>
  mutate(
    # "No" first, so the model estimates the probability of "Yes"
    HeartDisease = factor(HeartDisease, levels = c("No", "Yes")),
    # Approximate age (start of the group + 2 years; e.g. "55-59" -> 57) to read the effect per year
    age = as.numeric(str_sub(AgeCategory, 1, 2)) + 2,
    # "Excellent" as reference: ORs show how much risk rises as health worsens
    GenHealth = factor(GenHealth, levels = c("Excellent", "Very good", "Good", "Fair", "Poor"))
  ) |>
  select(-AgeCategory)

# Description of the selected variables: type, missing values, distribution
skimr::skim(model_data)

Data summary
Name	model_data
Number of rows	319795
Number of columns	12
_______________________
Column type frequency:
character	7
factor	2
numeric	3
________________________
Group variables	None

Variable type: character

skim_variable	complete_rate	min	max	n_unique
Sex	1	4	6	2
Smoking	1	2	3	2
Stroke	1	2	3	2
Diabetic	1	2	23	4
KidneyDisease	1	2	3	2
DiffWalking	1	2	3	2
PhysicalActivity	1	2	3	2

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
HeartDisease	0	1	FALSE	2	No: 292422, Yes: 27373
GenHealth	0	1	FALSE	5	Ver: 113858, Goo: 93129, Exc: 66842, Fai: 34677

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
BMI	1	28.33	6.36	12.02	24.03	27.34	31.42	94.85	▇▅▁▁▁
PhysicalHealth	1	3.37	7.95	0.00	0.00	0.00	2.00	30.00	▇▁▁▁▁
age	1	54.44	18.07	20.00	42.00	57.00	67.00	82.00	▆▅▇▇▇

3.4.2 Train/test split

# Seed: the random split is reproducible
set.seed(2026)

# Stratifying by HeartDisease keeps the same 8.6% of cases in both sets
data_split <- initial_split(model_data, prop = 0.75, strata = HeartDisease)
train_data <- training(data_split)
test_data  <- testing(data_split)

tibble(
  set       = c("Training", "Test"),
  people    = c(nrow(train_data), nrow(test_data)),
  pct_cases = c(mean(train_data$HeartDisease == "Yes"),
                mean(test_data$HeartDisease == "Yes")) * 100
)

3.4.3 Model 1: logistic regression (fitting)

# Common formula for caret and glm: outcome ~ predictors
model_formula <- HeartDisease ~ age + Sex + Smoking + Stroke + Diabetic +
  KidneyDisease + DiffWalking + PhysicalActivity + GenHealth + BMI + PhysicalHealth

# 5-fold cross-validation: training data is split into 5 parts and the model is
# fitted 5 times, each time evaluated on the part it did not see.
# This estimates performance without touching the test set yet.
cv_control <- caret::trainControl(
  method = "cv",
  number = 5,
  classProbs = TRUE,                        # needed to compute AUC
  summaryFunction = caret::twoClassSummary  # reports AUC under the name "ROC"
)

set.seed(2026)
cv_model <- caret::train(
  model_formula,
  data = train_data,
  method = "glm",            # logistic regression
  family = "binomial",
  trControl = cv_control,
  metric = "ROC"
)

# Mean AUC across the 5 folds and its variability
cv_model$results |> select(ROC, ROCSD)

# The same model with glm(), the function caret uses internally;
# its output is easier to read and diagnose
logit_model <- glm(model_formula, data = train_data, family = binomial)

# Odds ratios with 95% Wald confidence intervals
or_table <- broom::tidy(logit_model) |>
  filter(term != "(Intercept)") |>
  mutate(
    OR      = exp(estimate),
    CI_low  = exp(estimate - 1.96 * std.error),
    CI_high = exp(estimate + 1.96 * std.error)
  ) |>
  select(term, OR, CI_low, CI_high, p.value) |>
  arrange(desc(OR))

knitr::kable(or_table, digits = 3,
             caption = "Odds ratios of the model (sorted from highest to lowest)")

Odds ratios of the model (sorted from highest to lowest)
term	OR	CI_low	CI_high	p.value
GenHealthPoor	6.660	6.078	7.297	0.000
GenHealthFair	4.525	4.205	4.871	0.000
StrokeYes	2.863	2.721	3.013	0.000
GenHealthGood	2.797	2.617	2.989	0.000
SexMale	1.967	1.904	2.033	0.000
KidneyDiseaseYes	1.774	1.679	1.875	0.000
DiabeticYes	1.596	1.537	1.657	0.000
GenHealthVery good	1.557	1.454	1.667	0.000
SmokingYes	1.446	1.400	1.493	0.000
DiffWalkingYes	1.235	1.186	1.287	0.000
DiabeticYes (during pregnancy)	1.218	0.966	1.534	0.095
DiabeticNo, borderline diabetes	1.171	1.066	1.286	0.001
age	1.055	1.054	1.057	0.000
PhysicalActivityYes	1.039	1.002	1.078	0.038
BMI	1.012	1.009	1.014	0.000
PhysicalHealth	1.007	1.005	1.009	0.000

# Named vectors to quote values inside the text
or_vec      <- setNames(or_table$OR, or_table$term)
ci_low_vec  <- setNames(or_table$CI_low, or_table$term)
ci_high_vec <- setNames(or_table$CI_high, or_table$term)

Interpretation of the coefficients

All associations go in the expected direction except physical activity and gestational diabetes (few cases, no clear association). Holding the other variables constant:

Self-rated health carries the most weight: compared with excellent health, the odds are multiplied by 6.7 with poor health and by 4.5 with fair health.
Prior stroke multiplies the odds by 2.9; being male by 2.0.
Kidney disease and diabetes multiply them by 1.8 and 1.6; smoking and difficulty walking by 1.4 and 1.2.
Age: every 10 additional years multiply the odds by 1.7.
BMI and poor physical health days have small effects per unit (OR 1.012 per BMI point and 1.007 per day).
Physical activity is no longer associated with risk once age and health status are accounted for (OR = 1.04; 95% CI 1.00 to 1.08). The gap seen in chapter 2 (13.8% vs 7.1%) was largely explained by those other variables: older and sicker people exercise less. This is the difference between a simple and an adjusted association.

3.4.4 Preliminary diagnostics

In linear regression one checks constant variance (homoscedasticity) and approximately normal residuals. Those assumptions do not apply to logistic regression, whose outcome is yes/no. Instead, four aspects are checked: (1) predictors not too correlated with each other (multicollinearity), (2) extreme residuals, (3) influential observations and (4) calibration (predicted probabilities match observed proportions).

# Variance inflation factor (VIF): values above 5 signal problematic
# multicollinearity. For multi-level variables read the last column,
# GVIF^(1/(2*Df)), comparable with the square root of 5 (about 2.2)
vif_model <- car::vif(logit_model)
vif_model

##                      GVIF Df GVIF^(1/(2*Df))
## age              1.099392  1        1.048519
## Sex              1.052128  1        1.025733
## Smoking          1.032113  1        1.015930
## Stroke           1.033226  1        1.016477
## Diabetic         1.139387  3        1.021987
## KidneyDisease    1.048568  1        1.023996
## DiffWalking      1.445711  1        1.202377
## PhysicalActivity 1.175277  1        1.084102
## GenHealth        1.895001  4        1.083181
## BMI              1.162888  1        1.078373
## PhysicalHealth   1.725014  1        1.313398

# augment() adds each person's fitted value, residual and influence
diagnostics <- broom::augment(logit_model)

# The activity mentions geom_residuals(), which does not exist in ggplot2:
# the residual chart is built with geom_point()
c(pct_extreme_residuals = mean(abs(diagnostics$.std.resid) > 3) * 100,
  max_cooks_distance    = max(diagnostics$.cooksd))

## pct_extreme_residuals    max_cooks_distance 
##           0.108402892           0.000872242

# A sample of 3,000 people is plotted: 240,000 points would be unreadable
set.seed(2026)
plot_resid <- diagnostics |>
  slice_sample(n = 3000) |>
  ggplot(aes(x = .fitted, y = .std.resid, color = HeartDisease)) +
  geom_point(alpha = 0.4) +
  geom_hline(yintercept = 0) +
  geom_hline(yintercept = c(-3, 3), linetype = "dashed") +
  scale_color_manual(values = c(No = "#3B528B", Yes = "#D95F02")) +
  labs(
    title = "Standardized residuals versus fitted values (sample of 3,000)",
    x = "Fitted value (logit scale)",
    y = "Standardized deviance residual",
    color = "Heart disease"
  ) +
  theme_minimal()

ggplotly(plot_resid)

# Calibration: test-set people are grouped into 10 deciles of predicted
# probability and compared with the real share of cases in each decile
prob_test <- predict(logit_model, newdata = test_data, type = "response")

calibration <- tibble(
  observed  = test_data$HeartDisease == "Yes",
  predicted = prob_test
) |>
  mutate(decile = ntile(predicted, 10)) |>
  group_by(decile) |>
  summarise(mean_predicted = mean(predicted), observed_rate = mean(observed))

plot_calibration <- ggplot(calibration, aes(x = mean_predicted, y = observed_rate)) +
  geom_abline(linetype = "dashed", color = "grey50") +   # perfect calibration
  geom_line(color = "#3B528B") +
  geom_point(color = "#3B528B", size = 3) +
  labs(
    title = "Calibration on the test set (by risk decile)",
    x = "Mean predicted probability",
    y = "Observed share with heart disease"
  ) +
  theme_minimal()

ggplotly(plot_calibration)

Interpretation of the diagnostics

Multicollinearity: none. The largest adjusted VIF is 1.31, far below 2.2, so each variable contributes its own information.
Extreme residuals: only 0.11% of people have a standardized residual above 3 in absolute value. They are people with heart disease to whom the model assigned a very low probability: atypical cases (for example, young people without risk factors) that a screening based on these questions would miss.
Influential observations: the maximum Cook’s distance is 0.0009, far below 1: no single person changes the model.
Calibration: points follow the diagonal; in the highest-risk decile the model predicts 36.9% and 36.4% is observed. The probabilities can be used as risk estimates, not only to rank people.

3.4.5 Model 1: performance evaluation

# Classification threshold: the share of cases in training (about 8.6%).
# A person is "high risk" if their probability exceeds the average frequency.
# With the usual 0.5 the model would almost never say "Yes", because the disease is uncommon.
threshold <- mean(train_data$HeartDisease == "Yes")

prob_train <- predict(logit_model, type = "response")

# Same metrics for any data set
evaluate <- function(observed, probability, threshold) {
  predicted_class <- factor(if_else(probability >= threshold, "Yes", "No"), levels = c("No", "Yes"))
  tibble(
    AUC           = yardstick::roc_auc_vec(observed, probability, event_level = "second"),
    Sensitivity   = yardstick::sens_vec(observed, predicted_class, event_level = "second"),
    Specificity   = yardstick::spec_vec(observed, predicted_class, event_level = "second"),
    Accuracy      = yardstick::accuracy_vec(observed, predicted_class),
    Pct_high_risk = mean(predicted_class == "Yes")
  )
}

metrics_train    <- evaluate(train_data$HeartDisease, prob_train, threshold)
metrics_test     <- evaluate(test_data$HeartDisease, prob_test, threshold)
# Baseline for comparison: saying nobody has the disease (probability 0 for all)
metrics_baseline <- evaluate(test_data$HeartDisease, rep(0, nrow(test_data)), 0.5)

metrics_table <- bind_rows(
  metrics_train    |> mutate(Set = "Training"),
  metrics_test     |> mutate(Set = "Test"),
  metrics_baseline |> mutate(Set = "Baseline (all 'No'), test")
) |>
  relocate(Set)

knitr::kable(metrics_table, digits = 3,
             caption = "Performance metrics: training versus test")

Performance metrics: training versus test
Set	AUC	Sensitivity	Specificity	Accuracy	Pct_high_risk
Training	0.839	0.784	0.742	0.746	0.302
Test	0.841	0.786	0.742	0.745	0.305
Baseline (all ‘No’), test	0.500	0.000	1.000	0.912	0.000

# Test-set confusion matrix: hits and errors with the chosen threshold
class_test <- factor(if_else(prob_test >= threshold, "Yes", "No"), levels = c("No", "Yes"))
confusion_test <- table(Predicted = class_test, Observed = test_data$HeartDisease)
confusion_test

##          Observed
## Predicted    No   Yes
##       No  54099  1500
##       Yes 18850  5500

# Yield: share of invited (high-risk) people who actually have heart disease
yield_targeted <- confusion_test["Yes", "Yes"] / sum(confusion_test["Yes", ])
yield_targeted

## [1] 0.2258727

3.4.6 Model 2: random forest (challenger model)

A good modeling practice is to compare the chosen model (champion) with a more flexible alternative (challenger). A random forest combines hundreds of decision trees, each built on a random subset of people and variables, and averages their votes. It can capture non-linear relationships and interactions automatically, but it is harder to explain. It is fitted with caret::train(method = "rf"), as suggested in the activity.

# A random forest on 240,000 people would take too long on a laptop: it is
# trained on a random sample of 20,000 people from the training set and
# evaluated on the same full test set as the logistic regression
set.seed(2026)
train_rf <- train_data |> slice_sample(n = 20000)

rf_model <- caret::train(
  model_formula,
  data = train_rf,
  method = "rf",
  trControl = caret::trainControl(method = "none", classProbs = TRUE),
  tuneGrid = data.frame(mtry = 3),  # variables tried at each split
  ntree = 300,                      # number of trees
  nodesize = 25                     # minimum people per final leaf: avoids memorizing
)

prob_rf <- predict(rf_model, newdata = test_data, type = "prob")[, "Yes"]
metrics_rf <- evaluate(test_data$HeartDisease, prob_rf, threshold)

# Which variables does the forest rely on most?
rf_importance <- caret::varImp(rf_model)$importance |>
  rownames_to_column("variable") |>
  arrange(desc(Overall)) |>
  slice_head(n = 8)

ggplot(rf_importance, aes(x = Overall, y = reorder(variable, Overall))) +
  geom_col(fill = "#1F78B4") +
  labs(
    title = "Random forest: most important variables",
    x = "Relative importance (0-100)", y = NULL
  ) +
  theme_minimal()

3.4.7 Model comparison and selection

# Same test set, same threshold: a fair comparison
bind_rows(
  metrics_test |> mutate(Model = "1. Logistic regression (champion)"),
  metrics_rf   |> mutate(Model = "2. Random forest (challenger)")
) |>
  relocate(Model) |>
  knitr::kable(digits = 3, caption = "Champion versus challenger on the test set")

Champion versus challenger on the test set
Model	AUC	Sensitivity	Specificity	Accuracy	Pct_high_risk
1. Logistic regression (champion)	0.841	0.786	0.742	0.745	0.305
2. Random forest (challenger)	0.776	0.335	0.947	0.894	0.077

# ROC curves: the closer to the top-left corner, the better the model ranks people
roc_data <- bind_rows(
  yardstick::roc_curve(tibble(truth = test_data$HeartDisease, p = prob_test), truth, p,
                       event_level = "second") |> mutate(model = "Logistic regression"),
  yardstick::roc_curve(tibble(truth = test_data$HeartDisease, p = prob_rf), truth, p,
                       event_level = "second") |> mutate(model = "Random forest")
) |>
  # Thousands of thresholds per curve: one point in 40 draws the same line
  # and keeps the interactive chart light
  group_by(model) |>
  slice(unique(c(seq(1, n(), by = 40), n()))) |>
  ungroup()

plot_roc <- ggplot(roc_data, aes(x = 1 - specificity, y = sensitivity, color = model)) +
  geom_abline(linetype = "dashed", color = "grey60") +
  geom_path(linewidth = 1) +
  scale_color_brewer(palette = "Set1") +
  labs(
    title = "ROC curves of both models on the test set",
    x = "False-positive rate (1 - specificity)", y = "Sensitivity", color = NULL
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

ggplotly(plot_roc)

Interpretation

The challenger loses. The fair comparison is the AUC, because it does not depend on the threshold: 0.841 for the logistic regression versus 0.776 for the random forest, a gap of 0.065. In the ROC chart the logistic curve lies above the forest’s over most of its length: for the same share of false positives it finds more cases.

Why sensitivity collapses to 33.5% (vs 78.6%). A random forest’s “probability” is the share of trees that vote “Yes”. Because only 8.6% of people have heart disease, most final leaves contain a majority of “No” and most trees vote “No”, so the forest’s probabilities are squeezed toward zero: with the same 8.5% threshold it flags only 7.7% of people and misses most cases. This is a known weakness of tree ensembles with imbalanced classes, not an error in the code. Its accuracy (89.4%) looks better only because it says “No” more often, the same trap described for the all-“No” baseline.

The forest’s most important variables are age, BMI, StrokeYes (this type of importance tends to favor numeric variables such as age and BMI), so both models agree on the main risk drivers.

Decision: the logistic regression is the champion. It ranks people better, its probabilities are calibrated (section on diagnostics), and it can be explained to the Board and to clinicians factor by factor (odds ratios). For management: a more complex model is not automatically a better one; it has to earn its place with results.

Limitation (and how a fairer challenger would be built): the forest was trained on a sample of 20,000 people, without tuning and without correcting the imbalance. Balanced sampling in each tree (sampsize/strata in randomForest), probability forests (ranger) or gradient boosting would probably narrow the gap; that test is listed under proposed improvements.

3.4.8 3D view of the model

# Probability surface by age and BMI for a reference profile: male, "good"
# health, physically active, no stroke, diabetes or kidney disease.
# Two surfaces are compared: smoker and non-smoker.
ages <- seq(20, 82, by = 2)
bmis <- seq(16, 50, by = 1)

surface_prob <- function(smoker) {
  grid <- expand_grid(
    age = ages, BMI = bmis,
    Sex = "Male", Smoking = smoker, Stroke = "No", Diabetic = "No",
    KidneyDisease = "No", DiffWalking = "No", PhysicalActivity = "Yes",
    GenHealth = "Good", PhysicalHealth = 0
  )
  prob <- predict(logit_model, newdata = grid, type = "response")
  # Rows = BMI and columns = age: the matrix format plotly expects for a surface
  matrix(prob, nrow = length(bmis), ncol = length(ages))
}

z_nonsmoker <- surface_prob("No")
z_smoker    <- surface_prob("Yes")

plot_ly(x = ages, y = bmis) |>
  add_surface(z = z_nonsmoker, name = "Non-smoker", colorscale = "Viridis",
              showscale = FALSE, opacity = 0.9) |>
  add_surface(z = z_smoker, name = "Smoker", colorscale = "Reds",
              showscale = FALSE, opacity = 0.7) |>
  layout(
    title = "Predicted probability by age and BMI: smoker (red) vs non-smoker",
    scene = list(
      xaxis = list(title = "Age (years)"),
      yaxis = list(title = "BMI (kg/m²)"),
      zaxis = list(title = "Probability")
    )
  )

# 3D cloud of 2,000 test-set people: age, BMI and predicted probability,
# colored by whether they actually have heart disease
set.seed(2026)
test_sample <- test_data |>
  mutate(probability = prob_test) |>
  slice_sample(n = 2000)

plot_ly(
  test_sample,
  x = ~age, y = ~BMI, z = ~probability,
  color = ~HeartDisease, colors = c("#3B528B", "#D95F02"),
  type = "scatter3d", mode = "markers", marker = list(size = 3)
) |>
  layout(
    title = "Predicted probability for 2,000 test-set adults (orange = heart disease)",
    scene = list(
      xaxis = list(title = "Age (years)"),
      yaxis = list(title = "BMI (kg/m²)"),
      zaxis = list(title = "Predicted probability")
    )
  )

3.4.9 Interpretation of results

Discrimination: test AUC is 0.841 (cross-validation: 0.839), above the 0.80 target: for a random pair of a person with and without heart disease, the model ranks the first higher in about 8 of 10 cases.
No overfitting: training AUC (0.839) and test AUC (0.841) are practically equal; the model should perform similarly on new people from a similar population.
Sensitivity: with a threshold of 8.5%, the model detects 78.6% of people with heart disease (target ≥ 70%) by flagging 30.5% of the population as high risk. It detects 5,500 cases and misses 1,500.
Yield per invitation: 22.6% of the people invited have heart disease, versus 8.8% if invitations were random: about 2.6 times more cases per invitation.
Cost of that result: specificity is 74.2%, so 18,850 people without the disease would be invited (false positives).
Why accuracy misleads: the model’s accuracy (74.5%) is lower than the baseline that says “No” to everyone (91.2%), yet the baseline detects no patients (sensitivity 0).
3D charts: the surface shows age as the strongest driver; BMI raises the probability only slightly; the smoker surface lies always above the non-smoker one, and the gap widens with age. In the point cloud, people with heart disease concentrate in the high-probability region, although some appear low: the cases the model misses.

All three success metrics are met (AUC ≥ 0.80; sensitivity ≥ 70%; train–test difference ≤ 0.02).

3.4.10 Proposed improvements

Feature engineering (interactions). Following the multiple linear regression reference reviewed in class (where the effect of smoking changes with BMI), test whether smoking weighs more in people with diabetes or whether the age effect differs by sex (Smoking:Diabetic, age:Sex).
Additional data. Add key clinical variables absent here (blood pressure, cholesterol, family history, medication) and retrain with the health plan’s own claims and electronic health records before using it for resource allocation.
A fairer challenger. Retrain the random forest on more data, correcting the class imbalance (balanced sampling per tree) and tuning mtry, or try gradient boosting; adopt it only if it beats the logistic regression’s AUC by enough to justify the loss of interpretability.
Cost-based threshold. If cardiology capacity is limited, raise the threshold to invite fewer people, accepting lower sensitivity; the choice should weigh the cost of a missed case against a false positive.

4 Data Storytelling: from evidence to decision

This chapter turns the evidence into a narrative for the Board, following the A³ model and the principles of the grammar of graphics: the right geometry for each message, visual hierarchy with theme_minimal() and scale_fill_brewer(), explicit annotations with annotate() and readable labels with ggrepel::geom_text_repel(). Each chart answers three questions: What does it show? Why does it matter? What should we do?

The exploratory charts of chapters 1–3 were built to find patterns; the charts below are designed to explain them to decision-makers.

4.1 Attention: the burden is concentrated

# Prevalence and cases by age group, highlighting the priority group (65+)
age_burden <- heart |>
  group_by(AgeCategory) |>
  summarise(
    prevalence = mean(HeartDisease == "Yes") * 100,
    cases      = sum(HeartDisease == "Yes")
  ) |>
  mutate(age_band = factor(if_else(AgeCategory %in% age_65_plus, "65 and older", "Under 65"),
                           levels = c("Under 65", "65 and older")))

share_cases_65 <- sum(age_burden$cases[age_burden$age_band == "65 and older"]) / sum(age_burden$cases)
share_pop_65   <- mean(heart$AgeCategory %in% age_65_plus)

ggplot(age_burden, aes(x = AgeCategory, y = prevalence, fill = age_band)) +
  geom_col(width = 0.75) +
  # Muted color for the rest, strong color for the priority group (visual hierarchy)
  scale_fill_brewer(palette = "Paired") +
  annotate("text", x = "80 or older", y = max(age_burden$prevalence) + 1.5,
           label = paste0(num(max(age_burden$prevalence)), "%"), fontface = "bold") +
  annotate("text", x = 1, y = 20, hjust = 0, size = 4,
           label = paste0("Adults 65+ are ", pct(share_pop_65, 0), " of the population\n",
                          "but ", pct(share_cases_65, 0), " of all heart-disease cases")) +
  labs(
    title = "Heart disease concentrates after age 65",
    subtitle = "Share of adults reporting coronary heart disease or a heart attack, by age group",
    x = "Age group (years)", y = "Prevalence (%)", fill = NULL,
    caption = "Source: CDC BRFSS 2020 (Kaggle cleaned version), n = 319,795 adults"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top",
        panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1))

What it shows. Prevalence rises from 0.6% at 18–24 to 22.6% at 80+. Adults aged 65+ are 35% of the population but concentrate 67% of all cases.

Why it matters. Invitations that ignore age spend most capacity where the disease is rare; two out of three cases are in one clearly identifiable group.

What to do. Age is a strong criterion, but it leaves one in three cases outside the list. The next charts show how combining risk factors finds more of them with the same number of invitations.

4.2 Analysis 1: risk versus volume

# For each risk group: probability within the group (risk), share of all
# cases (volume) and share of the population (bubble size)
risk_profile <- function(label, in_group) {
  tibble(
    group               = label,
    risk                = mean(heart$HeartDisease[in_group] == "Yes") * 100,
    share_of_cases      = sum(in_group & heart$HeartDisease == "Yes") /
                            sum(heart$HeartDisease == "Yes") * 100,
    share_of_population = mean(in_group) * 100
  )
}

risk_volume <- bind_rows(
  risk_profile("Age 65+",             heart$AgeCategory %in% age_65_plus),
  risk_profile("Prior stroke",        heart$Stroke == "Yes"),
  risk_profile("Kidney disease",      heart$KidneyDisease == "Yes"),
  risk_profile("Diabetes",            heart$Diabetic == "Yes"),
  risk_profile("Ever smoker",         heart$Smoking == "Yes"),
  risk_profile("Male",                heart$Sex == "Male"),
  risk_profile("Difficulty walking",  heart$DiffWalking == "Yes"),
  risk_profile("Fair/poor health",    heart$GenHealth %in% c("Fair", "Poor")),
  risk_profile("Physically inactive", heart$PhysicalActivity == "No"),
  risk_profile("Obesity (BMI 30+)",   heart$BMI >= 30)
) |>
  # Strategic role of each group for program design
  mutate(role = case_when(
    share_of_cases >= 50 ~ "Volume driver (> half of cases)",
    risk >= 20           ~ "High-risk niche (risk >= 20%)",
    TRUE                 ~ "Contributing factor"
  ))

overall_prevalence <- p_heart * 100

ggplot(risk_volume, aes(x = risk, y = share_of_cases)) +
  geom_hline(yintercept = 50, linetype = "dashed", color = "grey60") +
  geom_vline(xintercept = overall_prevalence, linetype = "dashed", color = "grey60") +
  geom_point(aes(size = share_of_population, fill = role),
             shape = 21, color = "grey30", alpha = 0.85) +
  ggrepel::geom_text_repel(aes(label = group), size = 3.6, seed = 2026, box.padding = 0.6) +
  scale_fill_brewer(palette = "Set2") +
  scale_size_area(max_size = 16, guide = "none") +
  annotate("text", x = overall_prevalence + 0.4, y = 75, hjust = 0, size = 3.3, color = "grey40",
           label = paste0("Average risk: ", num(overall_prevalence), "%")) +
  annotate("text", x = max(risk_volume$risk), y = 53, hjust = 1, size = 3.3, color = "grey40",
           label = "Above this line: groups holding more than half of all cases") +
  labs(
    title = "Age 65+ holds two-thirds of cases; stroke and kidney disease carry the highest risk",
    subtitle = "Risk within each group versus its share of all heart-disease cases (bubble size = share of population)",
    x = "Risk: % with heart disease within the group",
    y = "Volume: % of all heart-disease cases",
    fill = "Strategic role",
    caption = "Groups overlap (one person can belong to several). Source: CDC BRFSS 2020"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

What it shows. Groups split into two strategic types. Volume drivers (age 65+, men, ever smokers) contain most cases; high-risk niches (prior stroke 36.4%, kidney disease 29.3%, fair/poor health, difficulty walking, diabetes) have 2.5 to 4 times the average risk but are small groups.

Why it matters. A single program cannot serve both: a population-wide program reaches volume, while niches need active case-finding and clinical follow-up.

What to do. Combine a targeted invitation for the high-volume group with automatic flags for niche conditions (stroke, kidney disease) in existing care routes.

4.3 Analysis 2: the combined score beats every single-criterion rule

# Capture curve: people are sorted from highest to lowest predicted risk and we
# track the cumulative share of cases found as more people are invited
case_test <- test_data$HeartDisease == "Yes"

capture_full <- tibble(prob = prob_test, case = case_test) |>
  arrange(desc(prob)) |>
  mutate(
    pct_invited  = row_number() / n() * 100,
    pct_detected = cumsum(case) / sum(case) * 100
  )

# One point every 80 people: same curve, lighter chart
capture <- capture_full |> filter(row_number() %% 80 == 0)

# Single-criterion rules a program could use today, evaluated on the same test set
prior_diagnosis <- test_data$Stroke == "Yes" | test_data$Diabetic == "Yes" |
  test_data$KidneyDisease == "Yes"

rules <- tibble(
  rule = c("Age 65+", "Prior diagnosis\n(stroke, diabetes or kidney)", "Age 65+ or prior diagnosis"),
  flag = list(test_data$age >= 65, prior_diagnosis, test_data$age >= 65 | prior_diagnosis)
) |>
  mutate(
    pct_invited    = map_dbl(flag, \(f) mean(f) * 100),
    pct_detected   = map_dbl(flag, \(f) sum(f & case_test) / sum(case_test) * 100),
    # Cases the score would find inviting exactly the same number of people
    score_detected = map_dbl(pct_invited, \(p) capture_full$pct_detected[ceiling(p / 100 * nrow(capture_full))]),
    gain_points    = score_detected - pct_detected
  ) |>
  select(-flag)

op_invited  <- metrics_test$Pct_high_risk * 100
op_detected <- metrics_test$Sensitivity * 100

plot_capture <- ggplot(capture, aes(x = pct_invited, y = pct_detected)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey50") +
  geom_line(color = "#1F78B4", linewidth = 1.2) +
  geom_segment(data = rules, aes(x = pct_invited, xend = pct_invited,
                                 y = pct_detected, yend = score_detected),
               color = "#FF7F00", linetype = "dotted") +
  geom_point(data = rules, aes(x = pct_invited, y = pct_detected, text = rule),
             shape = 21, fill = "#FF7F00", color = "grey20", size = 4) +
  annotate("point", x = op_invited, y = op_detected, size = 4, color = "#E31A1C") +
  annotate("text", x = op_invited - 2, y = op_detected + 7, hjust = 1, size = 3.6,
           label = paste0("Proposed: top ", num(op_invited, 0), "% by score\n-> ",
                          num(op_detected, 0), "% of cases")) +
  annotate("text", x = 80, y = 66, size = 3.4, color = "grey40", label = "Random invitations") +
  scale_x_continuous(limits = c(0, 100)) +
  scale_y_continuous(limits = c(0, 100)) +
  labs(
    title = "The combined risk score finds more cases than any single-criterion rule",
    subtitle = "Blue: cases found by the score as more members are invited. Orange: single-criterion rules (dotted = cases they miss)",
    x = "% of population invited",
    y = "% of heart-disease cases found",
    caption = "Logistic regression with 11 predictors; test set (25% of the data)"
  ) +
  theme_minimal(base_size = 12)

# Static version with readable, non-overlapping labels for the rules
plot_capture +
  ggrepel::geom_text_repel(data = rules, aes(x = pct_invited, y = pct_detected, label = rule),
                           size = 3.3, seed = 2026, nudge_y = -9, box.padding = 0.5)

rules |>
  mutate(rule = str_replace_all(rule, "\n", " ")) |>
  knitr::kable(
    digits = 1,
    col.names = c("Invitation rule", "% invited", "% of cases found (rule)",
                  "% of cases found (score, same invitations)", "Gain (percentage points)"),
    caption = "Single-criterion rules versus the combined score with the same number of invitations"
  )

Single-criterion rules versus the combined score with the same number of invitations
Invitation rule	% invited	% of cases found (rule)	% of cases found (score, same invitations)	Gain (percentage points)
Age 65+	34.9	67.0	83.2	16.3
Prior diagnosis (stroke, diabetes or kidney)	17.6	47.3	58.8	11.5
Age 65+ or prior diagnosis	42.6	82.2	89.0	6.8

Interactive version. The Board may prefer a different capacity level than the one proposed. The interactive chart (tooltips, zoom and pan) lets each manager read how many cases the score would find for any share of invitations, and compare it with each rule, without asking the analyst.

ggplotly(plot_capture, tooltip = c("text", "x", "y"))

What it shows. Every single-criterion rule (orange) lies below the score’s curve: with the same number of invitations, the score finds between 7 and 16 percentage points more cases. Inviting everyone aged 65+ finds 67% of cases; the score, inviting the same number of people, finds 83%. At the proposed operating point, the top 30% by score finds 79% of cases: 23 cases per 100 invitations versus 9 at random.

Why it matters. Programs already use risk factors, but one at a time. Combining them is what makes the difference: the same cardiology capacity reaches more of the people who need it.

What to do. Build the invitation list from the combined score, and let available capacity decide how far down the list the program goes.

4.4 Action: plan capacity on percentiles, not averages

# Cases per screening day, highlighting days that 9 slots (the average) would not cover
days <- tibble(cases = cases_per_day) |>
  count(cases) |>
  mutate(plan = factor(
    if_else(cases > 9, "Shortfall if planned on the average (9 slots)", "Covered by 9 slots"),
    levels = c("Shortfall if planned on the average (9 slots)", "Covered by 9 slots")
  ))

p95_slots <- unname(quantile(cases_per_day, 0.95))

ggplot(days, aes(x = cases, y = n, fill = plan)) +
  geom_col() +
  geom_vline(xintercept = p95_slots + 0.5, linetype = "dashed") +
  annotate("text", x = p95_slots + 0.8, y = max(days$n) * 0.85, hjust = 0, size = 3.8,
           label = paste0("Plan ", p95_slots, " slots:\ncovers 95% of days")) +
  scale_fill_brewer(palette = "Set1") +
  labs(
    title = paste0("Planning on the average leaves ", pct(mean(cases_per_day > 9), 0),
                   " of screening days short"),
    subtitle = "Heart-disease cases per screening day of 100 patients (1,000 simulated days)",
    x = "Cases per screening day", y = "Number of days", fill = NULL,
    caption = "Simulation with rbinom(); probability of heart disease = 8.6%"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")

What it shows. With 9 slots per day (the average), demand exceeds capacity on 35% of days; with 13 slots, on about 5%.

Why it matters. Unmet demand on screening days means patients sent home without assessment, re-scheduling costs and loss of trust in the program.

What to do. Staff each 100-patient screening day for 13 cardiology assessments and monitor the share of days with unmet demand.

4.5 Executive recommendations

Recommendation 1 — A 12-month, three-stage pilot of score-based invitations. Proposed by the VP of Population Health; decided by the Board at each checkpoint.

Stage	Months	What happens	Board checkpoint (go / adjust / stop)
0. Validation	1–3	Re-estimate and recalibrate the score with the plan’s claims and EHR data (HIPAA-compliant); measure the current rule’s baseline (invitations, attendance, cases found)	Month 3: continue only if, on the plan’s data, the score beats the current rule by ≥ 5 percentage points of cases found with the same invitations and AUC ≥ 0.75
1. Pilot	4–9	One region: invite members by score; automatic flags for stroke and kidney-disease patients; slots planned at the 95th percentile	Month 6: interim review of attendance, cases found and unmet-demand days
2. Evaluation	10–12	Compare with the baseline; estimate cost per case found	Month 12: scale up, adjust or stop

Recommendation 2 — Percentile-based capacity planning. Proposed by the VP of Operations: staff each 100-member screening day for 13 cardiology assessments (95th percentile) instead of the average, starting with the pilot region.

Realistic targets. The targets are deliberately below what the model achieved in this dataset, because models usually lose performance in a new population and because screening looks for undiagnosed people, who are harder to find:

KPI	In this dataset	Pilot target
Gain over the current rule (cases found, same invitations)	7–16 points	≥ 5 points
AUC on the plan’s own data	0.84	≥ 0.75
Sensitivity (cases in the high-risk list)	79%	≥ 70%
Attendance to the screening visit	Not measured	Set from the month-3 baseline
Screening days with unmet demand	35% if planned on the average	≤ 5%

These KPIs are monitored in the Board dashboard: https://rpubs.com/JuanSilva11/1470162

4.6 Limitations and how the pilot addresses them

Limitation	Consequence	How the pilot addresses it
National survey (BRFSS 2020), not the plan’s members; no sampling weights	Figures describe respondents, not the plan’s population	Stage 0 re-estimates the score with the plan’s data
Outcome is heart disease already diagnosed; screening seeks undiagnosed people	Real yield will be lower than in this dataset	Conservative targets, measured against the current rule
Cross-sectional, self-reported data	Associations, not causes	The score is used to prioritize, not to explain causes
No cost data	Savings cannot be claimed	The report speaks of efficiency (cases per invitation); cost per case is measured in stage 2
Capacity scenarios rely on assumptions	Figures are illustrative	Replaced with real attendance and admission data during the pilot

4.7 Beyond the pilot: transferring the method to Colombia

Why this section. What the Board is buying is not a number from a U.S. survey; it is a method: build a risk score, compare it with the current rule at equal invitations, plan capacity on percentiles and decide at checkpoints. A method that only works on one dataset is weak. This section tests whether it can travel to a very different system, Colombia’s. The method transfers; the numbers do not. Every figure in this report would have to be estimated again with Colombian data.

What changes in the Colombian system.

Component	United States (this report)	Colombia (transfer)
Who invites	Regional health plan	EPS (insurer), through its IPS network (providers)
Population	Plan members	Members of the contributory and subsidized regimes, with very different access conditions
Main data source	Self-reported survey (BRFSS)	RIPS (individual service records, now linked to electronic invoicing, Resolution 2275 of 2023), electronic health records, and the Cuenta de Alto Costo registry of hypertension, diabetes and chronic kidney disease
How risk factors are recorded	Yes/No answers	ICD-10 diagnosis codes: stroke I60–I69, diabetes E10–E14, chronic kidney disease N18, hypertension I10–I15; outcome: ischemic heart disease I20–I25
Where screening fits	Plan’s prevention program	The health promotion and maintenance pathway (RIAS, Resolution 3280 of 2018) and the EPS cardiovascular-risk programs
Data protection	HIPAA	Law 1581 of 2012 (personal data protection)

Transfer roadmap (proposal, not part of the 12-month U.S. decision).

Map the variables. Translate the 11 survey questions into Colombian records. Age, sex, stroke, diabetes and kidney disease can be taken from RIPS codes. BMI and smoking exist only in clinical records, often incomplete. Self-rated health and difficulty walking, two of the strongest predictors, are not routinely recorded and would need a short questionnaire in the promotion and maintenance consultation.
Redefine the outcome, and improve on this report. Instead of “already diagnosed” (the main limitation here), predict a new ischemic heart disease diagnosis or hospitalization in the next 12 months among members without one today. Claims data allow this; a one-time survey does not.
Re-estimate and validate retrospectively. Train the same two models (logistic regression and random forest) on historical data from one EPS, and compare them with the EPS’s current selection rule at equal invitations, as in Analysis 2.
Pilot within the existing pathway in one department, with the same go / adjust / stop logic, before scaling.

How much would yield change? Colombia’s prevalence of diagnosed heart disease in an EPS population is not known from this report and is likely lower than 8.6% (younger population, more underdiagnosis). The chunk below shows, as scenarios only, the yield per 100 invitations if the score kept the sensitivity and specificity measured here:

# Yield (positive predictive value) depends on prevalence, not only on the model:
# PPV = sens * p / (sens * p + (1 - spec) * (1 - p)).
# ASSUMPTION: the score keeps this report's sensitivity and specificity,
# which is optimistic; the prevalences are illustrative, not Colombian estimates
sens_score <- metrics_test$Sensitivity
spec_score <- metrics_test$Specificity

colombia_scenarios <- tibble(prevalence = c(0.02, 0.04, 0.06, mean(test_data$HeartDisease == "Yes"))) |>
  mutate(
    yield_score  = sens_score * prevalence /
      (sens_score * prevalence + (1 - spec_score) * (1 - prevalence)),
    yield_random = prevalence,
    lift         = yield_score / yield_random,
    Scenario     = if_else(row_number() == n(), "This dataset (U.S.)", "Colombian scenario")
  ) |>
  relocate(Scenario)

colombia_scenarios |>
  mutate(
    prevalence   = pct(prevalence),
    yield_score  = num(100 * yield_score),
    yield_random = num(100 * yield_random),
    lift         = paste0(num(lift), "x")
  ) |>
  knitr::kable(
    col.names = c("Scenario", "Prevalence", "Cases per 100 invitations (score)",
                  "Cases per 100 invitations (random)", "Lift"),
    caption = "Illustrative yield of the score under lower prevalences"
  )

Illustrative yield of the score under lower prevalences
Scenario	Prevalence	Cases per 100 invitations (score)	Cases per 100 invitations (random)	Lift
Colombian scenario	2.0%	5.8	2.0	2.9x
Colombian scenario	4.0%	11.2	4.0	2.8x
Colombian scenario	6.0%	16.3	6.0	2.7x
This dataset (U.S.)	8.8%	22.6	8.8	2.6x

What it tells us. With a prevalence of 2%, the score would find only about 5.8 cases per 100 invitations, far fewer than here, but still about 2.9 times what random invitations would find. For management: in Colombia the absolute targets must be set lower; the argument for the score is the relative gain over the current rule, not the U.S. yield.

Errors and limitations to expect.

Risk	Why it matters	Mitigation
Diagnosis codes recorded for billing, not for research	Missing or vague codes (symptom codes instead of diagnoses) weaken the score	Audit a sample of records against the clinical history before training
Underdiagnosis and unequal access (rural areas, subsidized regime)	The model learns who gets diagnosed, not who is sick, and may rank under-served groups lower	Evaluate sensitivity separately by regime and region; do not deploy if a group is clearly worse served
Missing BMI, smoking and self-rated health in administrative data	The strongest predictors would be unavailable	Two versions: a records-only score and a records-plus-questionnaire score; keep the one that beats the current rule
Members switching EPS	Incomplete history for part of the population	Require a minimum enrollment period in the training data
Younger age structure	An age-only rule covers fewer people, and calibration from U.S. data would be wrong	Recalibrate on local data; never reuse U.S. coefficients
Data protection and consent	Using health data to invite people requires a legal basis and clear purpose	Legal review under Law 1581 of 2012 before stage 0

The regulatory references above should be confirmed with the EPS’s legal team; they are cited as the starting point, not as legal advice.

References

Centers for Disease Control and Prevention (CDC). Behavioral Risk Factor Surveillance System (BRFSS), 2020 Survey Data and Documentation. https://www.cdc.gov/brfss/annual_data/annual_2020.html
Personal Key Indicators of Heart Disease [dataset]. Kaggle. https://www.kaggle.com/datasets/kamilpytlak/personal-key-indicators-of-heart-disease
García Díaz, C. M. (2025). Probabilidad y estadística aplicada; Modelos de analítica de datos; Storytelling con RStudio [course materials]. Fundación Universitaria de Ciencias de la Salud (FUCS).
García Sabaté, D. Regresión lineal múltiple: predicción de gastos médicos. RPubs. https://rpubs.com/DavidGS/rlm
Martínez Bautista, H. Explorando medicaldata en R. RPubs. https://rpubs.com/Humberto_Mtz_Bta/1229095
Data4HR. Visualización de datos con ggplot2. RPubs. https://rpubs.com/Data4HR/visualizacion-humanos-reales
James, G., Witten, D., Hastie, T. & Tibshirani, R. (2021). An Introduction to Statistical Learning with Applications in R (2nd ed.), ch. 4. Springer.
Nussbaumer Knaflic, C. (2019). Storytelling con datos: ejemplos prácticos. Anaya Multimedia.
Breiman, L. (2001). Random forests. Machine Learning, 45(1), 5–32.
Congreso de la República de Colombia. Ley 1581 de 2012, por la cual se dictan disposiciones generales para la protección de datos personales.
Ministerio de Salud y Protección Social de Colombia. Resolución 3280 de 2018 (Rutas Integrales de Atención en Salud: promoción y mantenimiento de la salud). https://normograma.supersalud.gov.co/compilacion/docs/resolucion_minsaludps_3280_2018.htm
Ministerio de Salud y Protección Social de Colombia. Resolución 2275 de 2023 (Registro Individual de Prestación de Servicios de Salud, RIPS, como soporte de la factura electrónica). https://www.saludcapital.gov.co/DPYS/Normatividad/Resoluciones/Res_2275-2023.pdf

From Data to Decision: A Risk-Targeted Cardiovascular Screening Strategy

Executive report for the Board — Heart disease indicators (CDC BRFSS 2020)

Juan José Silva Camargo · Paula Andrea López Rendón

October 24, 2026

Executive summary

Strategic framing: the A³ model

Attention: the management problem

Audience

North Star and objectives

1 Exploratory Data Analysis

1.1 The data

1.1.1 Structure and data quality

1.2 Descriptive indicators

1.2.1 Exploratory chart

1.2.2 Interpretation

2 Probability and Applied Statistics

2.1 Applied probability

2.1.1 Random selection and sample space

Interpretation

2.1.2 Empirical and conditional probabilities

Interpretation

2.1.3 Discrete distributions: capacity scenarios

Interpretation

2.1.4 Continuous distributions: BMI, arrival times and visit length

Interpretation

2.2 Inferential statistics

2.2.1 t-test: BMI by heart disease status

Interpretation

2.2.2 One-sample t-test: sleep hours

Interpretation

2.2.3 ANOVA: BMI by self-rated health

Interpretation

2.2.4 Chi-square: smoking and heart disease

Interpretation

2.2.5 Test of proportions: women and men

Interpretation

2.2.6 Comparative charts with ggpubr

Interpretation

3 Analytical Models: who is at high risk?

3.1 Problem, objectives and success metrics

3.2 Variable selection

3.3 Conceptual model

3.4 Implementation in RStudio

3.4.1 Data preparation

3.4.2 Train/test split

3.4.3 Model 1: logistic regression (fitting)

Interpretation of the coefficients

3.4.4 Preliminary diagnostics

Interpretation of the diagnostics

3.4.5 Model 1: performance evaluation

3.4.6 Model 2: random forest (challenger model)

3.4.7 Model comparison and selection

Interpretation

3.4.8 3D view of the model

3.4.9 Interpretation of results

3.4.10 Proposed improvements

4 Data Storytelling: from evidence to decision

4.1 Attention: the burden is concentrated

4.2 Analysis 1: risk versus volume

4.3 Analysis 2: the combined score beats every single-criterion rule

4.4 Action: plan capacity on percentiles, not averages

4.5 Executive recommendations

4.6 Limitations and how the pilot addresses them

4.7 Beyond the pilot: transferring the method to Colombia

References