Social Capital, Digital Engagement, and Psychological Well-being

1 Packages

# Only widely-available, stable packages
needed <- c("tidyverse",'rio',"knitr", "broom", "patchwork", "scales", "corrplot")
to_get <- needed[!needed %in% rownames(installed.packages())]
if (length(to_get)) install.packages(to_get, repos = "https://cloud.r-project.org")

library(tidyverse)
library(knitr)
library(rio)
library(broom)
library(patchwork)
library(scales)
library(corrplot)

---

2 Research Objectives

This study uses dd (n = 5,000) to examine psychological well-being across urban and non-urban communities.

Outcome: High life satisfaction = life_satisfaction >= 7 (binary, 0/1)

Objectives:

1. Describe the sample by residence (Table 1)
2. Estimate well-being prevalence across subgroups (Table 2)
3. Identify determinants via logistic regression — urban / non-urban / total (Table 3)
4. Predict well-being probabilities by sex × income × residence (Table 4)
5. Decompose the urban–non-urban gap (Blinder-Oaxaca, Tables 5–6)

---

3 Step 1 — Load and Recode Data

# ── dd must already be in your environment ──────────────────────────────────
dd=import('https://raw.githubusercontent.com/datalake101/hssSampleData/refs/heads/main/healthData.csv')

df <- dd %>%
  mutate(
    # Outcome
    high_wb = as.integer(life_satisfaction >= 7),

    # Residence: 2-level
    res2 = if_else(region == "Urban", "Urban", "Non-urban"),
    res2 = factor(res2, levels = c("Urban", "Non-urban")),

    # Age groups
    age_grp = case_when(
      age <= 24 ~ "18-24",
      age <= 44 ~ "25-44",
      age <= 64 ~ "45-64",
      TRUE      ~ "65+"
    ),
    age_grp = factor(age_grp, levels = c("18-24","25-44","45-64","65+")),

    # Education collapsed
    edu = case_when(
      education %in% c("None","Primary")        ~ "None/Primary",
      education == "Secondary"                  ~ "Secondary",
      education %in% c("Vocational","Graduate") ~ "Tertiary",
      education == "Postgrad"                   ~ "Postgraduate"
    ),
    edu = factor(edu, levels = c("None/Primary","Secondary","Tertiary","Postgraduate")),

    # Income tertile
    inc3 = factor(ntile(monthly_income, 3), labels = c("Low","Medium","High")),

    # Marital
    mar = case_when(
      marital_status == "Married"                  ~ "Married",
      marital_status == "Partnered"                ~ "Partnered",
      marital_status == "Single"                   ~ "Single",
      marital_status %in% c("Divorced","Widowed")  ~ "Div/Widowed"
    ),
    mar = factor(mar, levels = c("Married","Partnered","Single","Div/Widowed")),

    # Sex
    sex = factor(sex, levels = c("Male","Female","Non-binary")),

    # Social capital composite (1-5 scale items)
    sc_index = rowMeans(cbind(community_belonging, neighbor_trust,
                               family_support, friend_support,
                               perceived_cohesion), na.rm = TRUE),
    sc_cat = case_when(
      sc_index <  2   ~ "Low",
      sc_index <  3.5 ~ "Medium",
      TRUE            ~ "High"
    ),
    sc_cat = factor(sc_cat, levels = c("Low","Medium","High")),

    # Loneliness
    lone_cat = case_when(
      loneliness_frequency <= 2 ~ "Low",
      loneliness_frequency <= 3 ~ "Moderate",
      TRUE                      ~ "High"
    ),
    lone_cat = factor(lone_cat, levels = c("Low","Moderate","High")),

    # Screen time
    screen_cat = case_when(
      daily_screen_time <  2 ~ "Low",
      daily_screen_time <  5 ~ "Moderate",
      TRUE                   ~ "High"
    ),
    screen_cat = factor(screen_cat, levels = c("Low","Moderate","High")),

    # Mindfulness
    mind_cat = case_when(
      mindfulness_practice %in% c("Never","Rarely") ~ "Infrequent",
      mindfulness_practice == "Sometimes"            ~ "Sometimes",
      TRUE                                           ~ "Regular"
    ),
    mind_cat = factor(mind_cat, levels = c("Infrequent","Sometimes","Regular")),

    # Simple Yes/No factors
    can_confide   = factor(can_confide,   levels = c("No","Yes")),
    club_member   = factor(club_member,   levels = c("No","Yes")),
    voted         = factor(voted_last_election, levels = c("No","Yes")),
    donated       = factor(donated_charity,     levels = c("No","Yes")),
    uses_sm       = factor(uses_social_media,   levels = c("No","Yes")),
    shares_news   = factor(shares_news,         levels = c("No","Yes")),
    climate_sm    = factor(climate_info_on_sm,  levels = c("No","Yes")),
    therapy       = factor(recent_therapy,      levels = c("No","Yes")),
    discrim       = factor(experienced_discrimination, levels = c("No","Yes")),
    foreign_born  = factor(foreign_born,        levels = c("No","Yes")),
    ethnicity     = factor(ethnicity,
                           levels = c("White","Black","Latino",
                                      "Asian","Indigenous","Other")),

    # Volunteerism
    vol_cat = case_when(
      volunteerism == "Never"     ~ "Never",
      volunteerism == "Rarely"    ~ "Rarely",
      volunteerism == "Sometimes" ~ "Sometimes",
      TRUE                        ~ "Often"
    ),
    vol_cat = factor(vol_cat, levels = c("Never","Rarely","Sometimes","Often"))
  )

cat("Rows:", nrow(df), "\n")

## Rows: 5000

cat("High well-being:", sum(df$high_wb, na.rm=TRUE),
    "(", round(mean(df$high_wb, na.rm=TRUE)*100, 1), "%)\n")

## High well-being: 2020 ( 40.4 %)

cat("Urban:", sum(df$res2=="Urban"),
    "| Non-urban:", sum(df$res2=="Non-urban"), "\n")

## Urban: 1242 | Non-urban: 3758

---

4 Step 2 — Table 1: Sample Description by Residence

# Simple frequency table function
freq_by_res <- function(var_name) {
  data.frame(
    Category  = levels(df[[var_name]]),
    Overall   = paste0(table(df[[var_name]]),
                        " (", round(prop.table(table(df[[var_name]]))*100, 1), "%)"),
    Urban     = paste0(table(df[[var_name]][df$res2=="Urban"]),
                        " (", round(prop.table(table(df[[var_name]][df$res2=="Urban"]))*100, 1), "%)"),
    Non_urban = paste0(table(df[[var_name]][df$res2=="Non-urban"]),
                        " (", round(prop.table(table(df[[var_name]][df$res2=="Non-urban"]))*100, 1), "%)")
  )
}

vars_desc <- c("age_grp","sex","edu","inc3","mar",
               "sc_cat","lone_cat","screen_cat","mind_cat",
               "can_confide","club_member","voted","donated",
               "uses_sm","shares_news","discrim","therapy")

# Build table with section headers
t1_list <- lapply(vars_desc, function(v) {
  header <- data.frame(Category  = paste0("-- ", v, " --"),
                       Overall   = "", Urban = "", Non_urban = "")
  body   <- freq_by_res(v)
  rbind(header, body)
})

t1 <- do.call(rbind, t1_list)
rownames(t1) <- NULL

# Continuous summaries at top
cont_sum <- data.frame(
  Category  = c("-- Age (mean, SD) --",
                 "Age",
                 "-- Social Capital Index (mean, SD) --",
                 "Social Capital Index"),
  Overall   = c("", paste0(round(mean(df$age,      na.rm=TRUE),1),
                             " (", round(sd(df$age,      na.rm=TRUE),1), ")"),
                "", paste0(round(mean(df$sc_index,  na.rm=TRUE),1),
                             " (", round(sd(df$sc_index,  na.rm=TRUE),1), ")")),
  Urban     = c("", paste0(round(mean(df$age[df$res2=="Urban"],     na.rm=TRUE),1),
                             " (", round(sd(df$age[df$res2=="Urban"],     na.rm=TRUE),1), ")"),
                "", paste0(round(mean(df$sc_index[df$res2=="Urban"], na.rm=TRUE),1),
                             " (", round(sd(df$sc_index[df$res2=="Urban"], na.rm=TRUE),1), ")")),
  Non_urban = c("", paste0(round(mean(df$age[df$res2=="Non-urban"],     na.rm=TRUE),1),
                             " (", round(sd(df$age[df$res2=="Non-urban"],     na.rm=TRUE),1), ")"),
                "", paste0(round(mean(df$sc_index[df$res2=="Non-urban"], na.rm=TRUE),1),
                             " (", round(sd(df$sc_index[df$res2=="Non-urban"], na.rm=TRUE),1), ")"))
)

t1_full <- rbind(cont_sum, t1)

kable(t1_full,
      caption = "Table 1. Sample Characteristics by Residence",
      col.names = c("Variable / Category", "Overall", "Urban", "Non-urban"),
      row.names = FALSE)

Table 1. Sample Characteristics by Residence

Variable / Category	Overall	Urban	Non-urban
– Age (mean, SD) –
Age	53.8 (20.7)	54.8 (20.3)	53.5 (20.8)
– Social Capital Index (mean, SD) –
Social Capital Index	3 (0.6)	3 (0.6)	3 (0.6)
– age_grp –
18-24	436 (8.7%)	85 (6.8%)	351 (9.3%)
25-44	1420 (28.4%)	346 (27.9%)	1074 (28.6%)
45-64	1436 (28.7%)	381 (30.7%)	1055 (28.1%)
65+	1708 (34.2%)	430 (34.6%)	1278 (34%)
– sex –
Male	2337 (46.7%)	587 (47.3%)	1750 (46.6%)
Female	2603 (52.1%)	642 (51.7%)	1961 (52.2%)
Non-binary	60 (1.2%)	13 (1%)	47 (1.3%)
– edu –
None/Primary	1666 (33.3%)	409 (32.9%)	1257 (33.4%)
Secondary	817 (16.3%)	198 (15.9%)	619 (16.5%)
Tertiary	1675 (33.5%)	428 (34.5%)	1247 (33.2%)
Postgraduate	842 (16.8%)	207 (16.7%)	635 (16.9%)
– inc3 –
Low	1667 (33.3%)	415 (33.4%)	1252 (33.3%)
Medium	1667 (33.3%)	409 (32.9%)	1258 (33.5%)
High	1666 (33.3%)	418 (33.7%)	1248 (33.2%)
– mar –
Married	998 (20%)	254 (20.5%)	744 (19.8%)
Partnered	1019 (20.4%)	241 (19.4%)	778 (20.7%)
Single	988 (19.8%)	247 (19.9%)	741 (19.7%)
Div/Widowed	1995 (39.9%)	500 (40.3%)	1495 (39.8%)
– sc_cat –
Low	216 (4.3%)	57 (4.6%)	159 (4.2%)
Medium	3667 (73.3%)	896 (72.1%)	2771 (73.7%)
High	1117 (22.3%)	289 (23.3%)	828 (22%)
– lone_cat –
Low	2069 (41.4%)	486 (39.1%)	1583 (42.1%)
Moderate	960 (19.2%)	246 (19.8%)	714 (19%)
High	1971 (39.4%)	510 (41.1%)	1461 (38.9%)
– screen_cat –
Low	1095 (21.9%)	256 (20.6%)	839 (22.3%)
Moderate	1628 (32.6%)	397 (32%)	1231 (32.8%)
High	2277 (45.5%)	589 (47.4%)	1688 (44.9%)
– mind_cat –
Infrequent	2030 (40.6%)	488 (39.3%)	1542 (41%)
Sometimes	943 (18.9%)	237 (19.1%)	706 (18.8%)
Regular	2027 (40.5%)	517 (41.6%)	1510 (40.2%)
– can_confide –
No	1270 (25.4%)	314 (25.3%)	956 (25.4%)
Yes	3730 (74.6%)	928 (74.7%)	2802 (74.6%)
– club_member –
No	3835 (76.7%)	923 (74.3%)	2912 (77.5%)
Yes	1165 (23.3%)	319 (25.7%)	846 (22.5%)
– voted –
No	1704 (34.1%)	424 (34.1%)	1280 (34.1%)
Yes	3296 (65.9%)	818 (65.9%)	2478 (65.9%)
– donated –
No	3575 (71.5%)	882 (71%)	2693 (71.7%)
Yes	1425 (28.5%)	360 (29%)	1065 (28.3%)
– uses_sm –
No	735 (14.7%)	168 (13.5%)	567 (15.1%)
Yes	4265 (85.3%)	1074 (86.5%)	3191 (84.9%)
– shares_news –
No	4053 (81.1%)	1023 (82.4%)	3030 (80.6%)
Yes	947 (18.9%)	219 (17.6%)	728 (19.4%)
– discrim –
No	4608 (92.2%)	1153 (92.8%)	3455 (91.9%)
Yes	392 (7.8%)	89 (7.2%)	303 (8.1%)
– therapy –
No	4092 (81.8%)	1004 (80.8%)	3088 (82.2%)
Yes	908 (18.2%)	238 (19.2%)	670 (17.8%)

---

5 Step 3 — Table 2: Prevalence of High Well-being by Subgroup

# Function: prevalence (%) by one variable and by residence, with chi-sq p
prev_by_var <- function(var_name) {
  d <- df[!is.na(df[[var_name]]) & !is.na(df$high_wb), ]

  # Overall prevalence per category
  ov <- tapply(d$high_wb, d[[var_name]], function(x) round(mean(x)*100, 1))

  # Urban
  du <- d[d$res2 == "Urban",]
  ur <- tapply(du$high_wb, du[[var_name]], function(x) round(mean(x)*100, 1))

  # Non-urban
  dn <- d[d$res2 == "Non-urban",]
  nu <- tapply(dn$high_wb, dn[[var_name]], function(x) round(mean(x)*100, 1))

  # Chi-square p (variable vs high_wb, ignore residence)
  p <- round(chisq.test(table(d[[var_name]], d$high_wb))$p.value, 4)
  p_fmt <- if (p < 0.001) "<0.001" else as.character(p)

  cats <- levels(d[[var_name]])
  out  <- data.frame(
    Category  = cats,
    Overall   = as.numeric(ov[cats]),
    Urban     = as.numeric(ur[cats]),
    Non_urban = as.numeric(nu[cats]),
    p_value   = c(p_fmt, rep("", length(cats)-1)),
    stringsAsFactors = FALSE
  )

  header <- data.frame(Category = paste0("-- ", var_name, " --"),
                       Overall  = NA, Urban = NA, Non_urban = NA,
                       p_value  = "", stringsAsFactors = FALSE)
  rbind(header, out)
}

vars_t2 <- c("age_grp","sex","edu","inc3","mar",
             "sc_cat","lone_cat","screen_cat","mind_cat",
             "can_confide","club_member","voted","donated",
             "uses_sm","discrim","therapy")

t2 <- do.call(rbind, lapply(vars_t2, prev_by_var))
rownames(t2) <- NULL

kable(t2,
      caption = "Table 2. Prevalence of High Well-being (%) by Subgroup and Residence",
      col.names = c("Variable / Category", "Overall (%)", "Urban (%)",
                    "Non-urban (%)", "p-value"),
      na.string = "")

Table 2. Prevalence of High Well-being (%) by Subgroup and Residence

Variable / Category	Overall (%)	Urban (%)	Non-urban (%)	p-value
– age_grp –	NA	NA	NA
18-24	42.2	43.5	41.9	0.8519
25-44	40.6	41.6	40.3
45-64	40.1	37.0	41.2
65+	40.0	38.1	40.6
– sex –	NA	NA	NA
Male	40.1	37.1	41.0	0.4477
Female	40.9	40.8	40.9
Non-binary	33.3	46.2	29.8
– edu –	NA	NA	NA
None/Primary	39.4	36.9	40.2	0.3206
Secondary	39.9	41.9	39.3
Tertiary	42.2	40.2	42.9
Postgraduate	39.3	38.6	39.5
– inc3 –	NA	NA	NA
Low	40.9	37.8	41.9	0.8159
Medium	39.8	36.2	41.0
High	40.5	43.3	39.5
– mar –	NA	NA	NA
Married	40.7	40.6	40.7	0.6529
Partnered	42.0	40.7	42.4
Single	39.8	34.4	41.6
Div/Widowed	39.7	40.0	39.7
– sc_cat –	NA	NA	NA
Low	38.9	26.3	43.4	0.5975
Medium	40.1	40.0	40.2
High	41.6	39.1	42.5
– lone_cat –	NA	NA	NA
Low	40.5	39.7	40.7	0.9357
Moderate	39.9	38.2	40.5
High	40.6	39.0	41.1
– screen_cat –	NA	NA	NA
Low	39.8	37.5	40.5	0.9022
Moderate	40.5	40.6	40.5
High	40.6	38.9	41.2
– mind_cat –	NA	NA	NA
Infrequent	39.5	34.8	40.9	0.0961
Sometimes	38.6	38.8	38.5
Regular	42.2	43.3	41.8
– can_confide –	NA	NA	NA
No	39.5	38.9	39.7	0.4836
Yes	40.7	39.2	41.2
– club_member –	NA	NA	NA
No	39.9	39.2	40.1	0.1762
Yes	42.1	38.9	43.4
– voted –	NA	NA	NA
No	39.1	35.1	40.5	0.2034
Yes	41.0	41.2	41.0
– donated –	NA	NA	NA
No	40.7	39.8	41.0	0.5149
Yes	39.6	37.5	40.4
– uses_sm –	NA	NA	NA
No	41.5	38.1	42.5	0.5384
Yes	40.2	39.3	40.5
– discrim –	NA	NA	NA
No	40.2	38.9	40.6	0.3833
Yes	42.6	41.6	42.9
– therapy –	NA	NA	NA
No	40.5	38.1	41.2	0.8616
Yes	40.1	43.3	39.0

---

6 Step 4 — Figure 1: Prevalence Charts

plot_prev <- function(var, title_txt) {
  df %>%
    filter(!is.na(.data[[var]])) %>%
    group_by(x = .data[[var]], res2) %>%
    summarise(pct = mean(high_wb, na.rm=TRUE)*100, .groups="drop") %>%
    ggplot(aes(x=x, y=pct, fill=res2)) +
    geom_col(position=position_dodge(0.7), width=0.6) +
    scale_fill_manual(values=c("Urban"="#2980B9","Non-urban"="#E07B54")) +
    scale_y_continuous(limits=c(0,100), labels=function(v) paste0(v,"%")) +
    labs(title=title_txt, x=NULL, y="High Well-being (%)", fill=NULL) +
    theme_minimal(base_size=11) +
    theme(legend.position="top",
          axis.text.x=element_text(angle=20, hjust=1, size=9))
}

p1 <- plot_prev("sc_cat",     "Social Capital")
p2 <- plot_prev("inc3",       "Income Tertile")
p3 <- plot_prev("lone_cat",   "Loneliness")
p4 <- plot_prev("screen_cat", "Screen Time")

(p1 + p2) / (p3 + p4) +
  plot_annotation(
    title   = "Figure 1. Prevalence of High Well-being by Key Domains and Residence",
    caption = "Non-urban = Suburban + Rural + Remote")

---

7 Step 5 — Table 3: Logistic Regression (3 Models)

fml <- high_wb ~ sex + age_grp + edu + inc3 + mar + foreign_born +
                 ethnicity + sc_cat + can_confide + club_member +
                 vol_cat + voted + donated + uses_sm +
                 screen_cat + shares_news + climate_sm +
                 lone_cat + mind_cat + therapy + discrim

# Model I: Non-urban only
m_nu <- glm(fml, data = filter(df, res2=="Non-urban"), family=binomial())

# Model II: Urban only
m_ur <- glm(fml, data = filter(df, res2=="Urban"),     family=binomial())

# Model III: Total, add residence
fml_tot <- update(fml, . ~ . + res2)
m_tot   <- glm(fml_tot, data=df, family=binomial())

cat("N Non-urban:", nobs(m_nu), " | N Urban:", nobs(m_ur),
    " | N Total:", nobs(m_tot), "\n")

## N Non-urban: 3758  | N Urban: 1242  | N Total: 5000

# Helper: OR table with stars
or_table <- function(mod) {
  tidy(mod, conf.int=TRUE, exponentiate=TRUE) %>%
    filter(term != "(Intercept)") %>%
    mutate(
      star = case_when(p.value<0.001~"***", p.value<0.01~"**",
                       p.value<0.05~"*", TRUE~""),
      cell = paste0(sprintf("%.2f",estimate),
                    " [", sprintf("%.2f",conf.low),
                    "-",  sprintf("%.2f",conf.high), "]", star)
    ) %>%
    select(term, cell)
}

t3_nu  <- or_table(m_nu)  %>% rename(`Non-urban` = cell)
t3_ur  <- or_table(m_ur)  %>% rename(Urban       = cell)
t3_tot <- or_table(m_tot) %>% rename(Total       = cell)

t3 <- t3_nu %>%
  full_join(t3_ur,  by="term") %>%
  full_join(t3_tot, by="term") %>%
  # readable term labels
  mutate(term = term %>%
    str_replace("^sex",          "Sex: ") %>%
    str_replace("^age_grp",      "Age: ") %>%
    str_replace("^edu",          "Education: ") %>%
    str_replace("^inc3",         "Income: ") %>%
    str_replace("^mar",          "Marital: ") %>%
    str_replace("^foreign_born", "Foreign Born: ") %>%
    str_replace("^ethnicity",    "Ethnicity: ") %>%
    str_replace("^sc_cat",       "Social Capital: ") %>%
    str_replace("^can_confide",  "Can Confide: ") %>%
    str_replace("^club_member",  "Club Member: ") %>%
    str_replace("^vol_cat",      "Volunteering: ") %>%
    str_replace("^voted",        "Voted: ") %>%
    str_replace("^donated",      "Donated: ") %>%
    str_replace("^uses_sm",      "Uses Social Media: ") %>%
    str_replace("^screen_cat",   "Screen Time: ") %>%
    str_replace("^shares_news",  "Shares News: ") %>%
    str_replace("^climate_sm",   "Climate Info SM: ") %>%
    str_replace("^lone_cat",     "Loneliness: ") %>%
    str_replace("^mind_cat",     "Mindfulness: ") %>%
    str_replace("^therapy",      "Therapy: ") %>%
    str_replace("^discrim",      "Discrimination: ") %>%
    str_replace("^res2",         "Residence: ")
  )

kable(t3,
      caption = "Table 3. Adjusted Odds Ratios [95% CI] for High Well-being",
      col.names = c("Predictor","Non-urban AOR [95%CI]",
                    "Urban AOR [95%CI]","Total AOR [95%CI]"),
      align = "lrrr") %>%
  kableExtra::kable_styling(
    bootstrap_options = c("striped","condensed"), full_width=TRUE, font_size=12)

Table 3. Adjusted Odds Ratios [95% CI] for High Well-being

Predictor	Non-urban AOR [95%CI]	Urban AOR [95%CI]	Total AOR [95%CI]
Sex: Female	0.99 [0.87-1.13]	1.19 [0.94-1.51]	1.04 [0.92-1.16]
Sex: Non-binary	0.59 [0.30-1.09]	1.81 [0.56-5.75]	0.73 [0.41-1.24]
Age: 25-44	0.94 [0.73-1.20]	0.87 [0.53-1.44]	0.94 [0.76-1.17]
Age: 45-64	0.97 [0.76-1.25]	0.71 [0.43-1.16]	0.93 [0.74-1.15]
Age: 65+	0.94 [0.74-1.20]	0.76 [0.47-1.25]	0.91 [0.73-1.13]
Education: Secondary	0.96 [0.79-1.17]	1.24 [0.86-1.77]	1.02 [0.86-1.21]
Education: Tertiary	1.11 [0.95-1.31]	1.17 [0.88-1.57]	1.12 [0.98-1.29]
Education: Postgraduate	0.97 [0.80-1.19]	1.11 [0.78-1.58]	1.00 [0.84-1.19]
Income: Medium	0.96 [0.82-1.13]	0.90 [0.67-1.21]	0.96 [0.83-1.10]
Income: High	0.91 [0.78-1.07]	1.29 [0.97-1.72]	0.99 [0.86-1.14]
Marital: Partnered	1.08 [0.88-1.33]	0.98 [0.68-1.43]	1.07 [0.89-1.28]
Marital: Single	1.03 [0.83-1.27]	0.82 [0.57-1.20]	0.97 [0.81-1.16]
Marital: Div/Widowed	0.95 [0.80-1.14]	1.02 [0.74-1.40]	0.97 [0.83-1.13]
Foreign Born: Yes	0.90 [0.74-1.09]	1.02 [0.72-1.44]	0.94 [0.79-1.10]
Ethnicity: Black	1.01 [0.80-1.27]	0.92 [0.61-1.41]	0.98 [0.80-1.19]
Ethnicity: Latino	0.96 [0.77-1.21]	0.99 [0.65-1.50]	0.97 [0.79-1.18]
Ethnicity: Asian	1.04 [0.83-1.31]	0.91 [0.60-1.37]	0.99 [0.81-1.21]
Ethnicity: Indigenous	0.95 [0.75-1.19]	1.18 [0.80-1.75]	1.01 [0.83-1.23]
Ethnicity: Other	1.14 [0.91-1.44]	0.95 [0.63-1.43]	1.08 [0.89-1.32]
Social Capital: Medium	0.89 [0.64-1.23]	1.99 [1.09-3.83]*	1.05 [0.80-1.40]
Social Capital: High	0.97 [0.69-1.38]	1.94 [1.03-3.84]*	1.12 [0.83-1.52]
Can Confide: Yes	1.07 [0.92-1.25]	1.03 [0.78-1.35]	1.05 [0.93-1.20]
Club Member: Yes	1.15 [0.99-1.35]	0.98 [0.75-1.29]	1.11 [0.97-1.26]
Volunteering: Rarely	1.04 [0.88-1.22]	1.22 [0.91-1.63]	1.07 [0.93-1.24]
Volunteering: Sometimes	0.98 [0.83-1.17]	0.89 [0.66-1.20]	0.96 [0.83-1.11]
Volunteering: Often	0.94 [0.75-1.18]	0.91 [0.58-1.40]	0.93 [0.76-1.13]
Voted: Yes	1.03 [0.89-1.18]	1.37 [1.07-1.77]*	1.09 [0.97-1.23]
Donated: Yes	0.98 [0.85-1.14]	0.90 [0.69-1.17]	0.97 [0.85-1.10]
Uses Social Media: Yes	0.92 [0.76-1.10]	1.08 [0.77-1.53]	0.95 [0.81-1.11]
Screen Time: Moderate	0.99 [0.83-1.18]	1.17 [0.84-1.63]	1.03 [0.88-1.20]
Screen Time: High	1.03 [0.87-1.21]	1.07 [0.79-1.47]	1.04 [0.89-1.20]
Shares News: Yes	0.95 [0.80-1.12]	0.95 [0.70-1.29]	0.95 [0.82-1.10]
Climate Info SM: Yes	1.12 [0.95-1.31]	1.03 [0.78-1.38]	1.09 [0.95-1.26]
Loneliness: Moderate	0.99 [0.83-1.19]	0.90 [0.65-1.25]	0.97 [0.83-1.14]
Loneliness: High	1.03 [0.89-1.19]	0.97 [0.74-1.26]	1.01 [0.89-1.15]
Mindfulness: Sometimes	0.91 [0.76-1.09]	1.24 [0.89-1.73]	0.97 [0.83-1.14]
Mindfulness: Regular	1.04 [0.90-1.20]	1.50 [1.16-1.95]**	1.13 [1.00-1.28]
Therapy: Yes	0.91 [0.77-1.08]	1.25 [0.92-1.68]	0.99 [0.85-1.15]
Discrimination: Yes	1.09 [0.86-1.39]	1.08 [0.68-1.69]	1.10 [0.89-1.35]
Residence: Non-urban	NA	NA	1.08 [0.95-1.23]

7.0.1 Model fit

kable(data.frame(
  Model    = c("Non-urban","Urban","Total"),
  N        = c(nobs(m_nu), nobs(m_ur), nobs(m_tot)),
  AIC      = round(c(AIC(m_nu), AIC(m_ur), AIC(m_tot)), 1),
  R2_McFad = round(1 - c(m_nu$deviance/m_nu$null.deviance,
                           m_ur$deviance/m_ur$null.deviance,
                           m_tot$deviance/m_tot$null.deviance), 3)
), caption="Model Fit", col.names=c("Model","N","AIC","McFadden R2"))

Model Fit

Model	N	AIC	McFadden R2
Non-urban	3758	5134.5	0.005
Urban	1242	1698.3	0.027
Total	5000	6799.7	0.004

---

8 Step 6 — Figure 2: Forest Plot (Total Model)

tidy(m_tot, conf.int=TRUE, exponentiate=TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    sig  = p.value < 0.05,
    term = term %>%
      str_replace("^sex","Sex: ") %>%
      str_replace("^age_grp","Age: ") %>%
      str_replace("^edu","Education: ") %>%
      str_replace("^inc3","Income: ") %>%
      str_replace("^mar","Marital: ") %>%
      str_replace("^foreign_born","Foreign Born: ") %>%
      str_replace("^ethnicity","Ethnicity: ") %>%
      str_replace("^sc_cat","Social Capital: ") %>%
      str_replace("^can_confide","Can Confide: ") %>%
      str_replace("^club_member","Club Member: ") %>%
      str_replace("^vol_cat","Volunteering: ") %>%
      str_replace("^voted","Voted: ") %>%
      str_replace("^donated","Donated: ") %>%
      str_replace("^uses_sm","Uses Social Media: ") %>%
      str_replace("^screen_cat","Screen Time: ") %>%
      str_replace("^shares_news","Shares News: ") %>%
      str_replace("^climate_sm","Climate Info SM: ") %>%
      str_replace("^lone_cat","Loneliness: ") %>%
      str_replace("^mind_cat","Mindfulness: ") %>%
      str_replace("^therapy","Therapy: ") %>%
      str_replace("^discrim","Discrimination: ") %>%
      str_replace("^res2","Residence: ")
  ) %>%
  ggplot(aes(x=estimate, y=reorder(term,estimate), colour=sig)) +
  geom_vline(xintercept=1, linetype="dashed", colour="grey50") +
  geom_errorbarh(aes(xmin=conf.low, xmax=conf.high), height=0.3) +
  geom_point(aes(shape=sig), size=3) +
  scale_x_log10() +
  scale_colour_manual(values=c("TRUE"="#2C3E50","FALSE"="#BDC3C7"),
                      labels=c("TRUE"="p<0.05","FALSE"="ns"), name=NULL) +
  scale_shape_manual(values=c("TRUE"=16,"FALSE"=1),
                     labels=c("TRUE"="p<0.05","FALSE"="ns"), name=NULL) +
  labs(title="Figure 2. Adjusted Odds Ratios — Total Model",
       x="AOR (log scale)", y=NULL,
       caption="Error bars = 95% CI. Dashed line = OR of 1.") +
  theme_minimal(base_size=11) +
  theme(legend.position="top", panel.grid.minor=element_blank())

---

9 Step 7 — Table 4: Predicted Probabilities (Sex × Income × Residence)

# Mode helper
mode_val <- function(x) {
  ux <- unique(na.omit(x))
  ux[which.max(tabulate(match(x, ux)))]
}

# Build reference grid
grid <- expand.grid(
  sex      = c("Male","Female"),
  inc3     = c("Low","High"),
  res2     = c("Urban","Non-urban"),
  stringsAsFactors = FALSE
) %>%
  mutate(
    age_grp      = mode_val(df$age_grp) %>% as.character(),
    edu          = mode_val(df$edu)     %>% as.character(),
    mar          = mode_val(df$mar)     %>% as.character(),
    foreign_born = mode_val(df$foreign_born) %>% as.character(),
    ethnicity    = mode_val(df$ethnicity)    %>% as.character(),
    sc_cat       = mode_val(df$sc_cat)       %>% as.character(),
    can_confide  = mode_val(df$can_confide)  %>% as.character(),
    club_member  = mode_val(df$club_member)  %>% as.character(),
    vol_cat      = mode_val(df$vol_cat)      %>% as.character(),
    voted        = mode_val(df$voted)        %>% as.character(),
    donated      = mode_val(df$donated)      %>% as.character(),
    uses_sm      = mode_val(df$uses_sm)      %>% as.character(),
    screen_cat   = mode_val(df$screen_cat)   %>% as.character(),
    shares_news  = mode_val(df$shares_news)  %>% as.character(),
    climate_sm   = mode_val(df$climate_sm)   %>% as.character(),
    lone_cat     = mode_val(df$lone_cat)     %>% as.character(),
    mind_cat     = mode_val(df$mind_cat)     %>% as.character(),
    therapy      = mode_val(df$therapy)      %>% as.character(),
    discrim      = mode_val(df$discrim)      %>% as.character()
  )

# Match factor levels to df
for (col in names(grid)) {
  if (col %in% names(df) && is.factor(df[[col]])) {
    grid[[col]] <- factor(grid[[col]], levels = levels(df[[col]]))
  }
}

p_fit <- predict(m_tot, newdata=grid, type="response", se.fit=TRUE)

t4 <- grid %>%
  select(sex, inc3, res2) %>%
  mutate(
    Probability = round(p_fit$fit, 3),
    Percent     = paste0(round(p_fit$fit * 100, 1), "%"),
    CI_95       = paste0("[",
                         round((p_fit$fit - 1.96*p_fit$se.fit)*100, 1),
                         "% - ",
                         round((p_fit$fit + 1.96*p_fit$se.fit)*100, 1),
                         "%]")
  )

kable(t4,
      caption = "Table 4. Predicted Probability of High Well-being",
      col.names = c("Sex","Income","Residence","Probability","Percent","95% CI"))

Table 4. Predicted Probability of High Well-being

Sex	Income	Residence	Probability	Percent	95% CI
Male	Low	Urban	0.388	38.8%	[32% - 45.7%]
Female	Low	Urban	0.397	39.7%	[32.7% - 46.6%]
Male	High	Urban	0.386	38.6%	[31.7% - 45.4%]
Female	High	Urban	0.394	39.4%	[32.5% - 46.3%]
Male	Low	Non-urban	0.407	40.7%	[34.1% - 47.3%]
Female	Low	Non-urban	0.415	41.5%	[34.9% - 48.1%]
Male	High	Non-urban	0.404	40.4%	[33.9% - 47%]
Female	High	Non-urban	0.413	41.3%	[34.7% - 47.8%]

9.0.1 Figure 3: Predicted Probabilities Plot

grid %>%
  select(sex, inc3, res2) %>%
  mutate(fit   = p_fit$fit,
         lo    = p_fit$fit - 1.96*p_fit$se.fit,
         hi    = p_fit$fit + 1.96*p_fit$se.fit) %>%
  ggplot(aes(x=inc3, y=fit*100, colour=sex, group=sex)) +
  geom_line(linewidth=1) +
  geom_point(size=3) +
  geom_errorbar(aes(ymin=lo*100, ymax=hi*100), width=0.15) +
  facet_wrap(~res2) +
  scale_colour_manual(values=c("Male"="#2980B9","Female"="#C0392B",
                                "Non-binary"="#27AE60")) +
  scale_y_continuous(limits=c(0,100), labels=function(x) paste0(x,"%")) +
  labs(title="Figure 3. Predicted Well-being by Sex, Income, and Residence",
       x="Income Tertile", y="Predicted Probability (%)",
       colour="Sex", caption="Error bars = 95% CI.") +
  theme_minimal(base_size=12) +
  theme(legend.position="top", strip.text=element_text(face="bold"))

---

10 Step 8 — Blinder-Oaxaca Decomposition (Tables 5 & 6)

We implement the threefold decomposition manually using linear probability models (LPM), which avoids any dependency on the oaxaca package.

# Dummy-code all predictors
df_d <- df %>%
  filter(!is.na(high_wb)) %>%
  mutate(
    urban        = as.integer(res2 == "Urban"),
    d_female     = as.integer(sex == "Female"),
    d_nonbin     = as.integer(sex == "Non-binary"),
    d_a25        = as.integer(age_grp == "25-44"),
    d_a45        = as.integer(age_grp == "45-64"),
    d_a65        = as.integer(age_grp == "65+"),
    d_edu_sec    = as.integer(edu == "Secondary"),
    d_edu_ter    = as.integer(edu == "Tertiary"),
    d_edu_post   = as.integer(edu == "Postgraduate"),
    d_inc_med    = as.integer(inc3 == "Medium"),
    d_inc_hi     = as.integer(inc3 == "High"),
    d_mar_part   = as.integer(mar == "Partnered"),
    d_mar_sing   = as.integer(mar == "Single"),
    d_mar_dw     = as.integer(mar == "Div/Widowed"),
    d_foreign    = as.integer(foreign_born == "Yes"),
    d_sc_med     = as.integer(sc_cat == "Medium"),
    d_sc_hi      = as.integer(sc_cat == "High"),
    d_confide    = as.integer(can_confide == "Yes"),
    d_club       = as.integer(club_member == "Yes"),
    d_vol_rare   = as.integer(vol_cat == "Rarely"),
    d_vol_some   = as.integer(vol_cat == "Sometimes"),
    d_vol_often  = as.integer(vol_cat == "Often"),
    d_voted      = as.integer(voted == "Yes"),
    d_donated    = as.integer(donated == "Yes"),
    d_sm         = as.integer(uses_sm == "Yes"),
    d_scr_mod    = as.integer(screen_cat == "Moderate"),
    d_scr_hi     = as.integer(screen_cat == "High"),
    d_news       = as.integer(shares_news == "Yes"),
    d_climate    = as.integer(climate_sm == "Yes"),
    d_lone_mod   = as.integer(lone_cat == "Moderate"),
    d_lone_hi    = as.integer(lone_cat == "High"),
    d_mind_some  = as.integer(mind_cat == "Sometimes"),
    d_mind_reg   = as.integer(mind_cat == "Regular"),
    d_therapy    = as.integer(therapy == "Yes"),
    d_discrim    = as.integer(discrim == "Yes")
  )

# LPM formula (same predictors)
lpm_fml <- high_wb ~ d_female + d_nonbin +
  d_a25 + d_a45 + d_a65 +
  d_edu_sec + d_edu_ter + d_edu_post +
  d_inc_med + d_inc_hi +
  d_mar_part + d_mar_sing + d_mar_dw +
  d_foreign +
  d_sc_med + d_sc_hi + d_confide + d_club +
  d_vol_rare + d_vol_some + d_vol_often +
  d_voted + d_donated +
  d_sm + d_scr_mod + d_scr_hi + d_news + d_climate +
  d_lone_mod + d_lone_hi +
  d_mind_some + d_mind_reg +
  d_therapy + d_discrim

df_u  <- filter(df_d, urban == 1)
df_nu <- filter(df_d, urban == 0)

lpm_u  <- lm(lpm_fml, data = df_u)
lpm_nu <- lm(lpm_fml, data = df_nu)

# Group means
mean_u  <- mean(df_u$high_wb,  na.rm=TRUE)
mean_nu <- mean(df_nu$high_wb, na.rm=TRUE)
gap     <- mean_u - mean_nu

# Means of predictors in each group (exclude intercept)
X_u  <- colMeans(model.matrix(lpm_fml, data=df_u) [, -1])
X_nu <- colMeans(model.matrix(lpm_fml, data=df_nu)[, -1])

# Coefficients (exclude intercept)
B_u  <- coef(lpm_u) [-1]
B_nu <- coef(lpm_nu)[-1]

# Threefold decomposition
# E  = endowments  = (X_u - X_nu) * B_u
# C  = coefficients = X_nu * (B_u - B_nu)
# CE = interaction  = (X_u - X_nu) * (B_u - B_nu)

dX <- X_u  - X_nu
dB <- B_u  - B_nu

E_total  <- sum(dX * B_u,  na.rm=TRUE)
C_total  <- sum(X_nu * dB, na.rm=TRUE)
CE_total <- sum(dX * dB,   na.rm=TRUE)

cat("Gap:", round(gap,4),
    " | E:", round(E_total,4),
    " | C:", round(C_total,4),
    " | CE:", round(CE_total,4),
    " | Sum:", round(E_total+C_total+CE_total, 4), "\n")

## Gap: -0.0169  | E: -0.0011  | C: 0.2977  | CE: -0.0028  | Sum: 0.2938

10.1 Table 5 — Overall Decomposition

t5 <- data.frame(
  Component = c(
    "Urban — mean probability",
    "Non-urban — mean probability",
    "Total difference (Urban minus Non-urban)",
    "  Endowments effect (explained)",
    "  Coefficients effect (unexplained)",
    "  Interaction effect"
  ),
  Estimate = round(c(mean_u, mean_nu, gap,
                     E_total, C_total, CE_total), 4),
  Pct_of_gap = round(c(mean_u*100, mean_nu*100, 100,
                        E_total/gap*100,
                        C_total/gap*100,
                        CE_total/gap*100), 1),
  stringsAsFactors = FALSE
)

kable(t5,
      caption = "Table 5. Blinder-Oaxaca Decomposition: Urban vs Non-urban Well-being Gap",
      col.names = c("Component","Estimate","% of Total Gap"))

Table 5. Blinder-Oaxaca Decomposition: Urban vs Non-urban Well-being Gap

Component	Estimate	% of Total Gap
Urban — mean probability	0.3913	39.1
Non-urban — mean probability	0.4082	40.8
Total difference (Urban minus Non-urban)	-0.0169	100.0
Endowments effect (explained)	-0.0011	6.5
Coefficients effect (unexplained)	0.2977	-1762.5
Interaction effect	-0.0028	16.8

10.2 Table 6 — Detailed Decomposition by Predictor

# Variable-level contributions to the endowments effect
E_vars <- dX * B_u

# Clean labels
var_labels <- c(
  d_female="Female (vs Male)", d_nonbin="Non-binary (vs Male)",
  d_a25="Age 25-44", d_a45="Age 45-64", d_a65="Age 65+",
  d_edu_sec="Secondary Educ.", d_edu_ter="Tertiary Educ.",
  d_edu_post="Postgraduate Educ.",
  d_inc_med="Income: Medium", d_inc_hi="Income: High",
  d_mar_part="Marital: Partnered", d_mar_sing="Marital: Single",
  d_mar_dw="Marital: Div/Widowed",
  d_foreign="Foreign Born",
  d_sc_med="Social Capital: Medium", d_sc_hi="Social Capital: High",
  d_confide="Can Confide", d_club="Club Member",
  d_vol_rare="Volunteering: Rarely", d_vol_some="Volunteering: Sometimes",
  d_vol_often="Volunteering: Often",
  d_voted="Voted", d_donated="Donated to Charity",
  d_sm="Uses Social Media",
  d_scr_mod="Screen Time: Moderate", d_scr_hi="Screen Time: High",
  d_news="Shares News Online", d_climate="Shares Climate Info",
  d_lone_mod="Loneliness: Moderate", d_lone_hi="Loneliness: High",
  d_mind_some="Mindfulness: Sometimes", d_mind_reg="Mindfulness: Regular",
  d_therapy="Recent Therapy", d_discrim="Experienced Discrimination"
)

t6 <- data.frame(
  Variable   = var_labels[names(E_vars)],
  E_coef     = round(as.numeric(E_vars), 5),
  E_pct      = round(as.numeric(E_vars) / gap * 100, 1),
  C_coef     = round(as.numeric(X_nu * dB), 5),
  C_pct      = round(as.numeric(X_nu * dB) / gap * 100, 1),
  CE_coef    = round(as.numeric(dX  * dB), 5),
  CE_pct     = round(as.numeric(dX  * dB) / gap * 100, 1),
  stringsAsFactors = FALSE
)

# Sort by absolute endowments contribution
t6 <- t6[order(abs(t6$E_pct), decreasing=TRUE), ]
rownames(t6) <- NULL

kable(t6,
      caption = "Table 6. Detailed Decomposition: Contribution of Each Predictor",
      col.names = c("Variable",
                    "Endow. Coef","Endow. %",
                    "Coeff. Coef","Coeff. %",
                    "Inter. Coef","Inter. %"))

Table 6. Detailed Decomposition: Contribution of Each Predictor

Variable	Endow. Coef	Endow. %	Coeff. Coef	Coeff. %	Inter. Coef	Inter. %
Social Capital: Medium	-0.00239	14.1	0.13057	-773.0	-0.00282	16.7
Age 45-64	-0.00207	12.2	-0.02050	121.4	-0.00190	11.3
Social Capital: High	0.00179	-10.6	0.03314	-196.2	0.00186	-11.0
Mindfulness: Regular	0.00135	-8.0	0.03411	-201.9	0.00123	-7.3
Recent Therapy	0.00071	-4.2	0.01342	-79.5	0.00100	-5.9
Tertiary Educ.	0.00047	-2.8	0.00390	-23.1	0.00015	-0.9
Volunteering: Sometimes	-0.00042	2.5	-0.00586	34.7	-0.00037	2.2
Age 65+	-0.00038	2.3	-0.01616	95.6	-0.00029	1.7
Volunteering: Often	0.00039	-2.3	-0.00057	3.4	0.00010	-0.6
Screen Time: High	0.00040	-2.3	0.00455	-27.0	0.00025	-1.5
Volunteering: Rarely	-0.00036	2.2	0.01030	-61.0	-0.00029	1.7
Screen Time: Moderate	-0.00029	1.7	0.01295	-76.7	-0.00031	1.9
Non-binary (vs Male)	-0.00027	1.6	0.00315	-18.7	-0.00051	3.0
Income: High	0.00028	-1.6	0.02772	-164.1	0.00037	-2.2
Uses Social Media	0.00026	-1.6	0.03160	-187.1	0.00058	-3.4
Secondary Educ.	-0.00026	1.5	0.00967	-57.2	-0.00031	1.8
Age 25-44	0.00022	-1.3	-0.00470	27.8	0.00012	-0.7
Shares News Online	0.00022	-1.3	0.00012	-0.7	-0.00001	0.1
Female (vs Male)	-0.00020	1.2	0.02203	-130.4	-0.00021	1.2
Donated to Charity	-0.00017	1.0	-0.00652	38.6	-0.00015	0.9
Loneliness: Moderate	-0.00017	1.0	-0.00370	21.9	-0.00016	0.9
Loneliness: High	-0.00017	1.0	-0.00519	30.7	-0.00029	1.7
Mindfulness: Sometimes	0.00015	-0.9	0.01412	-83.6	0.00022	-1.3
Experienced Discrimination	-0.00013	0.8	-0.00062	3.7	0.00007	-0.4
Income: Medium	0.00013	-0.7	-0.00483	28.6	0.00008	-0.5
Marital: Single	-0.00008	0.5	-0.01022	60.5	-0.00009	0.5
Club Member	-0.00007	0.4	-0.00818	48.4	-0.00115	6.8
Postgraduate Educ.	-0.00005	0.3	0.00461	-27.3	-0.00006	0.4
Marital: Partnered	0.00005	-0.3	-0.00456	27.0	0.00029	-1.7
Voted	-0.00006	0.3	0.04516	-267.4	-0.00005	0.3
Foreign Born	-0.00003	0.2	0.00479	-28.4	-0.00015	0.9
Shares Climate Info	0.00003	-0.2	-0.01501	88.9	-0.00009	0.5
Marital: Div/Widowed	0.00002	-0.1	0.00661	-39.1	0.00008	-0.5
Can Confide	0.00001	-0.1	-0.00820	48.5	-0.00002	0.1

10.2.1 Figure 4: Endowments Waterfall

t6 %>%
  filter(abs(E_pct) > 0.5) %>%
  ggplot(aes(x=reorder(Variable, E_pct), y=E_pct,
             fill=E_pct > 0)) +
  geom_col(width=0.7) +
  geom_hline(yintercept=0) +
  coord_flip() +
  scale_fill_manual(values=c("TRUE"="#27AE60","FALSE"="#E74C3C"),
                    labels=c("TRUE"="Widens gap","FALSE"="Narrows gap")) +
  labs(title="Figure 4. Endowments Contribution to Urban-Non-urban Well-being Gap",
       x=NULL, y="% Contribution", fill=NULL,
       caption="|contribution| > 0.5% shown") +
  theme_minimal(base_size=12) +
  theme(legend.position="top", panel.grid.minor=element_blank())

---

11 Step 9 — Correlation Matrix

cont_cols <- df %>%
  select(life_satisfaction, happiness_today, stress_level,
         loneliness_frequency, social_anxiety,
         community_belonging, neighbor_trust,
         family_support, friend_support, perceived_cohesion,
         sc_index, daily_screen_time, number_of_platforms,
         optimism, gratitude, sense_of_purpose)

corrplot(cor(cont_cols, use="pairwise.complete.obs"),
         method="color", type="upper", order="hclust",
         addCoef.col="black", number.cex=0.6,
         tl.cex=0.75, tl.col="black",
         col=colorRampPalette(c("#C0392B","white","#2980B9"))(200),
         title="Figure 5. Correlation Matrix", mar=c(0,0,2,0))

---

12 Step 10 — Well-being by Region and Age

df %>%
  filter(!is.na(region), !is.na(age_grp)) %>%
  group_by(region, age_grp) %>%
  summarise(pct=mean(high_wb, na.rm=TRUE)*100, .groups="drop") %>%
  ggplot(aes(x=region, y=pct, fill=age_grp)) +
  geom_col(position=position_dodge(0.8), width=0.7) +
  scale_fill_brewer(palette="Set2", name="Age") +
  scale_y_continuous(limits=c(0,100), labels=function(x) paste0(x,"%")) +
  labs(title="Figure 6. High Well-being by Region and Age Group",
       x="Region", y="Prevalence (%)") +
  theme_minimal(base_size=12) +
  theme(legend.position="top")

---

13 Session Info

sessionInfo()

## R version 4.3.2 (2023-10-31)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.2.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Toronto
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] corrplot_0.95   scales_1.4.0    patchwork_1.3.0 broom_1.0.8    
##  [5] rio_1.2.3       knitr_1.50      lubridate_1.9.4 forcats_1.0.0  
##  [9] stringr_1.5.1   dplyr_1.1.4     purrr_1.1.0     readr_2.1.5    
## [13] tidyr_1.3.1     tibble_3.3.0    ggplot2_3.5.2   tidyverse_2.0.0
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10        generics_0.1.4     xml2_1.3.8         stringi_1.8.7     
##  [5] hms_1.1.3          digest_0.6.37      magrittr_2.0.3     evaluate_1.0.4    
##  [9] grid_4.3.2         timechange_0.3.0   RColorBrewer_1.1-3 fastmap_1.2.0     
## [13] R.oo_1.27.0        jsonlite_2.0.0     R.utils_2.12.3     backports_1.5.0   
## [17] viridisLite_0.4.2  textshaping_1.0.1  jquerylib_0.1.4    cli_3.6.5         
## [21] rlang_1.1.6        R.methodsS3_1.8.2  withr_3.0.2        cachem_1.1.0      
## [25] yaml_2.3.10        tools_4.3.2        tzdb_0.5.0         kableExtra_1.4.0  
## [29] curl_6.4.0         vctrs_0.6.5        R6_2.6.1           lifecycle_1.0.4   
## [33] MASS_7.3-60        pkgconfig_2.0.3    pillar_1.11.0      bslib_0.9.0       
## [37] gtable_0.3.6       data.table_1.17.8  glue_1.8.0         systemfonts_1.2.3 
## [41] xfun_0.52          tidyselect_1.2.1   rstudioapi_0.17.1  dichromat_2.0-0.1 
## [45] farver_2.1.2       htmltools_0.5.8.1  svglite_2.2.1      labeling_0.4.3    
## [49] rmarkdown_2.29     compiler_4.3.2