---
title: "Risk-Targeted Cardiovascular Screening · Board Dashboard"
author: "J. J. Silva Camargo · P. A. López Rendón"
output:
  flexdashboard::flex_dashboard:
    orientation: rows
    vertical_layout: fill
    theme: flatly
    source_code: embed
---

```{r setup, include=FALSE}
# Hide package messages: the dashboard should show only results
knitr::opts_chunk$set(message = FALSE, warning = FALSE)

library(flexdashboard)  # dashboard layout, valueBox() and gauge()
library(tidyverse)      # data wrangling and ggplot2 graphics
library(here)           # paths relative to the project folder
library(plotly)         # interactive and 3D charts
library(rsample)        # train/test split (same as the report)
library(ggrepel)        # non-overlapping labels

# Number formatting for values shown in boxes and text
num <- function(x, d = 1) formatC(as.numeric(x), format = "f", digits = d, big.mark = ",")
pct <- function(x, d = 1) paste0(num(100 * x, d), "%")

heart <- read_csv(here("datos", "heart_2020_cleaned.csv"), show_col_types = FALSE)

# ---- Burden KPIs --------------------------------------------------------------
p_heart        <- mean(heart$HeartDisease == "Yes")
age_65_plus    <- c("65-69", "70-74", "75-79", "80 or older")
share_pop_65   <- mean(heart$AgeCategory %in% age_65_plus)
share_cases_65 <- mean(heart$AgeCategory[heart$HeartDisease == "Yes"] %in% age_65_plus)

# ---- Targeting model: same variables, seed and split as the executive report --
model_data <- heart |>
  select(HeartDisease, AgeCategory, Sex, Smoking, Stroke, Diabetic, KidneyDisease,
         DiffWalking, PhysicalActivity, GenHealth, BMI, PhysicalHealth) |>
  drop_na() |>
  mutate(
    HeartDisease = factor(HeartDisease, levels = c("No", "Yes")),
    age = as.numeric(str_sub(AgeCategory, 1, 2)) + 2,
    GenHealth = factor(GenHealth, levels = c("Excellent", "Very good", "Good", "Fair", "Poor"))
  ) |>
  select(-AgeCategory)

set.seed(2026)
data_split <- initial_split(model_data, prop = 0.75, strata = HeartDisease)
train_data <- training(data_split)
test_data  <- testing(data_split)

model_formula <- HeartDisease ~ age + Sex + Smoking + Stroke + Diabetic +
  KidneyDisease + DiffWalking + PhysicalActivity + GenHealth + BMI + PhysicalHealth
logit_model <- glm(model_formula, data = train_data, family = binomial)

# Operating threshold = share of cases in training (prioritizes not missing cases)
threshold <- mean(train_data$HeartDisease == "Yes")
prob_test <- predict(logit_model, newdata = test_data, type = "response")
flag_test <- prob_test >= threshold
case_test <- test_data$HeartDisease == "Yes"

auc_test    <- yardstick::roc_auc_vec(test_data$HeartDisease, prob_test, event_level = "second")
sensitivity <- sum(flag_test & case_test) / sum(case_test)
pct_invited <- mean(flag_test)
yield       <- sum(flag_test & case_test) / sum(flag_test)

# ---- Capacity simulation (same seed as the report) ----------------------------
set.seed(2026)
cases_per_day      <- rbinom(n = 1000, size = 100, prob = p_heart)
p95_slots          <- unname(quantile(cases_per_day, 0.95))
pct_days_short_avg <- mean(cases_per_day > 9)
```

```{r charts, include=FALSE}
# Charts are built once here and displayed in the pages below

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

plot_age <- ggplot(age_burden, aes(x = AgeCategory, y = prevalence, fill = age_band)) +
  geom_col(width = 0.75) +
  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 = 3.8,
           label = paste0("Adults 65+: ", pct(share_pop_65, 0), " of the population,\n",
                          pct(share_cases_65, 0), " of all cases")) +
  labs(title = "Heart disease concentrates after age 65",
       x = "Age group (years)", y = "Prevalence (%)", fill = NULL,
       caption = "Source: CDC BRFSS 2020, 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))

# 2. Analysis 1: risk versus volume by risk group
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)
) |>
  mutate(role = case_when(
    share_of_cases >= 50 ~ "Volume driver",
    risk >= 20           ~ "High-risk niche",
    TRUE                 ~ "Contributing factor"
  ))

plot_risk <- ggplot(risk_volume, aes(x = risk, y = share_of_cases)) +
  geom_hline(yintercept = 50, linetype = "dashed", color = "grey60") +
  geom_vline(xintercept = p_heart * 100, 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.5, seed = 2026, box.padding = 0.6) +
  scale_fill_brewer(palette = "Set2") +
  scale_size_area(max_size = 15, guide = "none") +
  annotate("text", x = p_heart * 100 + 0.4, y = 75, hjust = 0, size = 3.2, color = "grey40",
           label = paste0("Average risk: ", pct(p_heart))) +
  labs(title = "Volume drivers versus high-risk niches",
       x = "Risk: % with heart disease within the group",
       y = "Volume: % of all heart-disease cases",
       fill = "Strategic role",
       caption = "Bubble size = share of population. Groups overlap.") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

# 3. Analysis 2: capture curve of the score versus single-criterion rules
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)
capture <- capture_full |> filter(row_number() %% 80 == 0)

prior_diagnosis <- test_data$Stroke == "Yes" | test_data$Diabetic == "Yes" |
  test_data$KidneyDisease == "Yes"

rules <- tibble(
  rule = c("Age 65+", "Prior diagnosis", "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),
    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)

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 = pct_invited * 100, y = sensitivity * 100, size = 4, color = "#E31A1C") +
  annotate("text", x = pct_invited * 100 - 2, y = sensitivity * 100 + 7, hjust = 1, size = 3.4,
           label = paste0("Proposed: top ", num(pct_invited * 100, 0), "% -> ",
                          num(sensitivity * 100, 0), "% of cases")) +
  annotate("text", x = 80, y = 66, size = 3.3, color = "grey40", label = "Random invitations") +
  scale_x_continuous(limits = c(0, 100)) +
  scale_y_continuous(limits = c(0, 100)) +
  labs(title = "Score (blue) versus single-criterion rules (orange)",
       x = "% of population invited", y = "% of heart-disease cases found") +
  theme_minimal(base_size = 12)

# 4. Action: screening-day capacity
days <- tibble(cases = cases_per_day) |>
  count(cases) |>
  mutate(plan = factor(
    if_else(cases > 9, "Shortfall with 9 slots (average)", "Covered by 9 slots"),
    levels = c("Shortfall with 9 slots (average)", "Covered by 9 slots")
  ))

plot_capacity <- 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.6,
           label = paste0("Plan ", p95_slots, " slots:\n95% of days covered")) +
  scale_fill_brewer(palette = "Set1") +
  labs(title = paste0("Planning on the average leaves ", pct(pct_days_short_avg, 0),
                      " of screening days short"),
       x = "Cases per screening day (100 patients)", y = "Number of days (of 1,000 simulated)",
       fill = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")
```

Executive overview {data-icon="fa-chart-bar"}
=====================================

Row {data-height=170}
-------------------------------------

### Prevalence

```{r}
valueBox(pct(p_heart), caption = "of adults report heart disease or a heart attack",
         icon = "fa-heartbeat", color = "danger")
```

### Concentration

```{r}
valueBox(pct(share_cases_65, 0),
         caption = paste0("of cases are adults 65+ (", pct(share_pop_65, 0), " of the population)"),
         icon = "fa-users", color = "warning")
```

### Model AUC (target ≥ 0.80)

```{r}
gauge(round(auc_test, 2), min = 0.5, max = 1,
      sectors = gaugeSectors(success = c(0.8, 1), warning = c(0.7, 0.8), danger = c(0.5, 0.7)))
```

### Sensitivity (target ≥ 70%)

```{r}
gauge(round(sensitivity * 100), min = 0, max = 100, symbol = "%",
      sectors = gaugeSectors(success = c(70, 100), warning = c(50, 70), danger = c(0, 50)))
```

### Gain over the age-65 rule

```{r}
valueBox(paste0("+", num(rules$gain_points[1], 0), " pts"),
         caption = "more cases found by the score than by 'Age 65+', same invitations",
         icon = "fa-bullseye", color = "success")
```

Row
-------------------------------------

### Where is the burden? (interactive: hover, zoom and pan)

```{r}
ggplotly(plot_age)
```

### Which groups matter, and how?

```{r}
plot_risk
```

The story (A³) {.storyboard data-icon="fa-book"}
=====================================

### **Attention.** Two out of three cases are adults aged 65+

```{r}
plot_age
```

***

**What it shows.** Prevalence rises from `r num(min(age_burden$prevalence))`% at 18–24 to `r num(max(age_burden$prevalence))`% at 80+. Adults 65+ are `r pct(share_pop_65, 0)` of the population but `r pct(share_cases_65, 0)` of all cases.

**Why it matters.** Invitations that ignore age spend most capacity where the disease is rare.

**What to do.** Age is a strong criterion, but it leaves one in three cases off the list: combine it with other risk factors (next frames).

### **Analysis 1.** Volume drivers need a program; high-risk niches need flags

```{r}
plot_risk
```

***

**What it shows.** Age 65+, men and ever-smokers contain most cases (*volume drivers*). Stroke, kidney disease, fair/poor health, difficulty walking and diabetes carry 2.5–4 times the average risk but are small groups (*high-risk niches*).

**Why it matters.** One program cannot serve both: volume needs population outreach; niches need active case-finding in clinical routes.

**What to do.** Targeted invitations for the volume group plus automatic flags for stroke and kidney-disease patients.

### **Analysis 2.** The combined score beats every single-criterion rule

```{r}
ggplotly(plot_capture, tooltip = c("text", "x", "y"))
```

***

**What it shows.** Every single-criterion rule (orange) lies below the score's curve (blue). With the same invitations, the score finds `r num(min(rules$gain_points), 0)`–`r num(max(rules$gain_points), 0)` percentage points more cases. "Age 65+" finds `r num(rules$pct_detected[1], 0)`% of cases; the score, inviting the same number of people, `r num(rules$score_detected[1], 0)`%.

**Why it matters.** Programs already use risk factors, but one at a time. Combining 11 simple questions (AUC = `r num(auc_test, 2)`) reaches more of the people who need care with the same capacity: `r round(yield * 100)` cases per 100 invitations versus `r round(p_heart * 100)` at random.

**What to do.** Build the invitation list from the score; capacity decides how far down the list to go. *Hover to compare any capacity level with each rule.*

### **Action.** Staff screening days at the 95th percentile, not the average

```{r}
plot_capacity
```

***

**What it shows.** With 9 slots (the average), demand exceeds capacity on `r pct(pct_days_short_avg, 0)` of days; with `r p95_slots` slots, on about 5%.

**Why it matters.** Unmet demand means patients sent home unassessed, rescheduling costs and loss of trust.

**What to do.** Staff each 100-patient day for `r p95_slots` cardiology assessments and track days with unmet demand.

3D risk views {data-icon="fa-cube"}
=====================================

Row
-------------------------------------

### Risk surface by age and BMI: smoker (red) vs non-smoker

```{r}
# Probability surface for a reference profile: male, "good" health, physically
# active, no stroke, diabetes or kidney disease
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
  )
  matrix(predict(logit_model, newdata = grid, type = "response"),
         nrow = length(bmis), ncol = length(ages))
}

plot_ly(x = ages, y = bmis) |>
  add_surface(z = surface_prob("No"), name = "Non-smoker", colorscale = "Viridis",
              showscale = FALSE, opacity = 0.9) |>
  add_surface(z = surface_prob("Yes"), name = "Smoker", colorscale = "Reds",
              showscale = FALSE, opacity = 0.7) |>
  layout(scene = list(xaxis = list(title = "Age (years)"),
                      yaxis = list(title = "BMI (kg/m²)"),
                      zaxis = list(title = "Probability")))
```

### Predicted risk for 2,000 test-set adults (orange = heart disease)

```{r}
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(scene = list(xaxis = list(title = "Age (years)"),
                      yaxis = list(title = "BMI (kg/m²)"),
                      zaxis = list(title = "Predicted probability")))
```

Row {data-height=150}
-------------------------------------

### How to read the surface

Age is the main driver: the surface climbs steeply along the age axis and only gently along BMI. The smoker surface (red) lies always above the non-smoker one, and the gap widens with age: **smoking cessation pays off most in older adults**.

### How to read the point cloud

Adults with heart disease (orange) concentrate at the top of the cloud, where predicted risk is highest. The orange points near the floor are the cases the model misses: atypical patients that **clinical flags** must catch.

Action plan {data-icon="fa-tasks"}
=====================================

Row
-------------------------------------

### Recommendation 1 · 12-month staged pilot of score-based invitations

- **Proposed by:** VP of Population Health. **Decided by:** the Board, at each checkpoint.
- **Months 1–3 · Validation:** re-estimate the score with the plan's claims and EHR data; measure the current rule's baseline. **Go/no-go:** the score must beat the current rule by ≥ 5 points of cases found (same invitations) and reach AUC ≥ 0.75.
- **Months 4–9 · Pilot** in one region, with automatic flags for stroke and kidney-disease patients. **Month 6:** interim review.
- **Months 10–12 · Evaluation:** cost per case found. **Month 12:** scale up, adjust or stop.

### Recommendation 2 · Percentile-based capacity plan

- **Proposed by:** VP of Operations.
- **Action:** staff each 100-member screening day for `r p95_slots` cardiology assessments (95th percentile), starting in the pilot region.
- **Target:** ≤ 5% of screening days with unmet demand (planning on the average: `r pct(pct_days_short_avg, 0)`).
- **Prerequisite:** replace the simulation with real attendance data during the pilot.

Row
-------------------------------------

### Why this model · two models compared

- **Model 1 · Logistic regression (champion):** AUC `r num(auc_test, 2)`, explained factor by factor with odds ratios.
- **Model 2 · Random forest (challenger):** lower AUC (0.78 in the report) and much lower sensitivity, because with only 8.6% cases most trees vote "No"; also harder to explain and audit.
- **Decision:** keep the logistic regression; a more complex model has to earn its place with results.
- **Next:** re-run the comparison on the plan's data in months 1–3.

### Beyond the pilot · transfer to Colombia

- **The method transfers; the numbers do not.** Re-estimate with EPS data (RIPS, clinical records, Cuenta de Alto Costo).
- **Better outcome:** predict *new* ischemic heart disease (I20–I25) in 12 months, not prior diagnosis.
- **Expected errors:** billing-oriented coding, underdiagnosis in rural and subsidized populations, missing self-rated health and BMI.
- **Rule:** lower absolute targets; judge the score by its gain over the EPS's current rule.

Row
-------------------------------------

### KPI monitoring

```{r}
# Pilot targets are deliberately below this dataset's results: models usually lose
# performance in a new population, and screening seeks undiagnosed people
kpis <- tibble(
  KPI = c("Gain over the current rule (cases found, same invitations)",
          "Model AUC",
          "Sensitivity (cases in the high-risk list)",
          "Attendance to the screening visit",
          "Screening days with unmet demand",
          "Cases per 100 invitations (score vs random)"),
  `In this dataset` = c(paste0(num(min(rules$gain_points), 0), "-", num(max(rules$gain_points), 0), " pts"),
                        num(auc_test, 2), pct(sensitivity, 0), "Not measured",
                        paste0(pct(pct_days_short_avg, 0), " (planned on average)"),
                        paste0(round(yield * 100), " vs ", round(p_heart * 100))),
  `Pilot target` = c(">= 5 pts", ">= 0.75", ">= 70%", "Set from month-3 baseline", "<= 5%",
                     "Monitor (expected lower: undiagnosed cases)"),
  Owner = c("Population Health", "Analytics", "Population Health", "Population Health",
            "Operations", "Population Health")
)

knitr::kable(kpis)
```