Load libraries.

# R markdown link: https://rpubs.com/jasohosk/1428122

# Load libraries
packages <- c("ideanet", "kableExtra", "tidyverse", "psych", "purrr", "ltm")

installed <- packages %in% rownames(installed.packages())
if (any(!installed)) {
  install.packages(packages[!installed])
}
lapply(packages, library, character.only = TRUE)

# Set working directory
setwd(
  "/Users/jasonhoskin/Library/CloudStorage/OneDrive-IndianaUniversity/study-sgm-caregivers"
)

Load survey/demographics dataset. Test for and remove any rows with NAs.

# Import dataset
url <- "2026-05-18_final-table_sans-sna.csv"
df0 <- read.csv(url)

## Assess for NAs
colSums(is.na(df0))
##                        id            current_gender                      race 
##                         0                         0                         0 
##                 ethnicity        sexual_orientation                 education 
##                         0                         0                         0 
##                  military          living_situation            marital_status 
##                         0                         0                         0 
##       biological_children   not_biological_children         employment_status 
##                         0                         0                         0 
##       work_hours_per_week                 happiness             companionship 
##                         0                         0                         0 
##                  left_out                  isolated                  contacts 
##                         0                         0                         0 
##             vigorous_days          vigorous_minutes              vigorous_met 
##                         0                         0                         0 
##             moderate_days          moderate_minutes              moderate_met 
##                         0                         0                         0 
##              walking_days           walking_minutes               walking_met 
##                         0                         0                         0 
##              sitting_days                   no_food                skip_meals 
##                         0                         0                         0 
##                no_shelter         personal_income_1         personal_income_2 
##                         0                         0                         0 
##         personal_income_3         personal_income_4         personal_income_5 
##                         0                         0                         0 
##         personal_income_6         personal_income_7         personal_income_8 
##                         0                         0                         0 
##        household_income_1        household_income_2        household_income_3 
##                         0                         0                         0 
##        household_income_4        household_income_5        household_income_6 
##                         0                         0                         0 
##        household_income_7        household_income_8            patient_demand 
##                         0                         0                         0 
##                 time_loss          patient_behavior             patient_anger 
##                         0                         0                         0 
##            patient_impact            patient_future        patient_dependency 
##                         0                         0                         0 
##           personal_income           food_insecurity        zarit_burden_score 
##                         0                         0                         1 
##  ipaq_sf_score_continuous ipaq_sf_score_categorical                loneliness 
##                         0                         0                         1 
##             perceived_qol 
##                         0
# Conduct listwise deletions
df0 <- na.omit(df0)

Calculate Cronbach’s for Caregiver burden and IPAQ.

# Caregiver burden
cronbach.alpha(df0[46:54])
## 
## Cronbach's alpha for the 'df0[46:54]' data-set
## 
## Items: 9
## Sample units: 19
## alpha: 0.753
# Food insecurity
cronbach.alpha(df0[18:20])
## 
## Cronbach's alpha for the 'df0[18:20]' data-set
## 
## Items: 3
## Sample units: 19
## alpha: 0.016
# IPAQ
cronbach.alpha(df0[8:17])
## 
## Cronbach's alpha for the 'df0[8:17]' data-set
## 
## Items: 10
## Sample units: 19
## alpha: -0.008

Read in Network Canvas data.

networks <- nc_read(
  path = "networkcanvasexport_merged")

# Exclude ego_id = 1 (missing data, excluded from analysis)
networks$egos             <- networks$egos             |> filter(ego_id != 1)
networks$alters           <- networks$alters           |> filter(ego_id != 1)
networks$alter_edgelists  <- networks$alter_edgelists  |> filter(ego_id != 1)

res <- ego_netwrite(networks$egos, 
            ego_id = "ego_id",
            networks$alters,
            alter_id = "alter_id",
            alter_ego = "ego_id",
            max_alters = Inf,
            alter_alter = networks$alter_edgelists,
            aa_ego = "ego_id",
            i_elements = "from",
            j_elements = "to")

Check results.

# Average across all networks
res$overall_summary |>
  knitr::kable() |>
  kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover")) |> 
  column_spec(1, bold = TRUE)
measure_labels measure_descriptions measures
Number of egos/ego networks Total number of egos providing ego networks in dataset 19
Number of alters Total number of alters nominated by egos across entire dataset 214
Number of isolates Number of egos who did not report any alters in their personal network 0
Number of one-node networks Number of egos who reported only one alter in their personal network 0
Smallest non-isolate network size Smallest number of alters provided by a single ego 4
Largest network size Largest number of alters provided by a single ego 21
Average network size Average number of alters provided by a single ego 11.2631578947368
Average network density The average density of personal networks provided by egos (networks with 0-1 alters excluded from calculation) 0.465546476227591
Average fragmentation The mean fragmentation index score of personal networks provided by egos (networks with 0-1 alters excluded from calculation) 0.215077354705838 
# Summaries for each network (ego)
res$summaries |>
  dplyr::select(ego_id, network_size, mean_degree, density, effective_size) |>
  knitr::kable() |>
  kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover")) |> 
  column_spec(1, bold = TRUE)
ego_id network_size mean_degree density effective_size
1 9 5.1111111 0.6388889 3.888889
2 15 3.4666667 0.2476190 11.533333
3 21 6.5714286 0.3285714 14.428571
4 9 6.2222222 0.7777778 2.777778
5 13 9.6923077 0.8076923 3.307692
6 4 1.5000000 0.5000000 2.500000
7 6 2.6666667 0.5333333 3.333333
8 6 0.6666667 0.1333333 5.333333
9 6 5.0000000 1.0000000 1.000000
10 18 3.2222222 0.1895425 14.777778
11 18 6.8888889 0.4052288 11.111111
12 11 2.0000000 0.2000000 9.000000
13 15 4.4000000 0.3142857 10.600000
14 12 2.5000000 0.2272727 9.500000
15 7 2.2857143 0.3809524 4.714286
16 9 2.8888889 0.3611111 6.111111
17 4 3.0000000 1.0000000 1.000000
18 17 5.4117647 0.3382353 11.588235
19 14 6.0000000 0.4615385 8.000000 
# Find ego IDs for the smallest and largest networks by network size
ego_ids_in_igraph <- sapply(res$igraph_objects, \(x) x$ego)

smallest_ego_id <- res$summaries |> slice_min(network_size, n = 1) |> pull(ego_id)
largest_ego_id  <- res$summaries |> slice_max(network_size, n = 1) |> pull(ego_id)

smallest_ego_id
## [1]  6 17
largest_ego_id
## [1] 3
## Because two egos are tied for smallest, I will manually insert their index number and print both.

# Plot side by side
plot(
  res$igraph_objects[[3]]$igraph_ego,
  main = paste0("Ego 3 — Largest")
)

plot(
  res$igraph_objects[[6]]$igraph_ego,
  main = paste0("Ego 6 — Smallest 1")
)

plot(
  res$igraph_objects[[17]]$igraph_ego,
  main = paste0("Ego 17 — Smallest 2")
)

Incorporate SNA data into df0.

# Extract SNA variables from ego_netwrite results
sna_vars <- res$summaries |>
  dplyr::select(
    ego_id,
    network_size,
    mean_degree,
    density,
    effective_size
  )
# Match IDs
df0$id = 1:19

## Join SNA results into df0
df1 <- df0 |>
  left_join(sna_vars, by = c("id" = "ego_id"))

Provide descriptive statistics for continuous variables.

# Create descriptive table
continuous_vars <- df1 |>
  dplyr::select(
    zarit_burden_score,
    ipaq_sf_score_continuous,
    network_size,
    effective_size,
    mean_degree,
    density,
    work_hours_per_week,
    loneliness,
    perceived_qol,
    biological_children,
    not_biological_children,
    contacts
  )

continuous_descriptives <- psych::describe(continuous_vars) |>
  as.data.frame() |>
  dplyr::select(n, mean, sd, median, min, max, skew, kurtosis)

continuous_descriptives |>
  mutate(across(where(is.numeric), \(x) round(x, 2))) |>
  knitr::kable(
    caption = "Descriptive Statistics for Continuous Variables",
    col.names = c("n", "Mean", "SD", "Median", "Min", "Max", "Skew", "Kurtosis"),
    align = c("r", "r", "r", "r", "r", "r", "r", "r")
  ) |>
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = FALSE,
    position = "left"
  )
Descriptive Statistics for Continuous Variables
n Mean SD Median Min Max Skew Kurtosis
zarit_burden_score 19 19.84 5.30 20.00 11.00 29.00 0.12 -0.84
ipaq_sf_score_continuous 19 3249.95 3794.43 1884.00 0.00 13488.00 1.90 2.45
network_size 19 11.26 5.16 11.00 4.00 21.00 0.21 -1.30
effective_size 19 7.08 4.43 6.11 1.00 14.78 0.23 -1.40
mean_degree 19 4.18 2.25 3.47 0.67 9.69 0.57 -0.36
density 19 0.47 0.26 0.38 0.13 1.00 0.77 -0.68
work_hours_per_week 19 27.37 19.55 35.00 0.00 60.00 -0.17 -1.52
loneliness 19 8.74 3.28 8.00 3.00 14.00 -0.24 -1.07
perceived_qol 19 6.11 2.00 7.00 2.00 9.00 -0.45 -1.09
biological_children 19 0.21 0.54 0.00 0.00 2.00 2.25 4.09
not_biological_children 19 0.42 1.12 0.00 0.00 4.00 2.31 3.96
contacts 19 571.42 1230.54 247.00 30.00 5562.00 3.49 11.25 
# Create histograms
continuous_vars <- c(
  "biological_children",
  "not_biological_children",
  "work_hours_per_week",
  "density",
  "network_size",
  "effective_size",
  "mean_degree",
  "ipaq_sf_score_continuous",
  "loneliness",
  "perceived_qol",
  "contacts"
)

df_hist_long <- df1 |>
  dplyr::select(all_of(continuous_vars)) |>
  pivot_longer(cols = everything(),
               names_to = "variable",
               values_to = "value")

ggplot(df_hist_long, aes(x = value)) +
  geom_histogram(bins = 30, color = "white") +
  facet_wrap(~ variable, scales = "free") +
  theme_bw()

# Create scatterplots of continuous variables by caregiver burden
df_scatter_long <- df1 |>
  dplyr::select(zarit_burden_score, all_of(continuous_vars)) |>
  pivot_longer(
    cols = all_of(continuous_vars),
    names_to = "variable",
    values_to = "value"
  )

ggplot(df_scatter_long, aes(x = value, y = zarit_burden_score)) +
  geom_point(alpha = 0.7) +
  facet_wrap(~ variable, scales = "free_x") +
  theme_bw() + 
  labs(x = "Value", y = "Caregiver Burden")

Provide descriptive statistics grouped by IPAQ category.

cat_vars <- c(
  "current_gender", "race", "ethnicity", "sexual_orientation",
  "education", "military", "living_situation", "marital_status",
  "employment_status", "personal_income", "food_insecurity"
)

cont_vars <- c(
  "zarit_burden_score", "ipaq_sf_score_continuous", "network_size", "effective_size",
  "mean_degree", "density", "work_hours_per_week", "loneliness",
  "perceived_qol", "biological_children", "not_biological_children", "contacts"
)

group_var <- "ipaq_sf_score_categorical"
groups <- c("Low", "Moderate", "High")

cat_rows <- map_dfr(cat_vars, function(var) {
  map_dfr(groups, function(grp) {
    df1 |>
      filter(.data[[group_var]] == grp) |>
      count(level = as.character(.data[[var]])) |>
      mutate(
        variable = var,
        group = grp,
        stat = sprintf("%d (%.1f%%)", n, 100 * n / sum(n))
      ) |>
      dplyr::select(variable, level, group, stat)
  }) |>
    pivot_wider(names_from = group, values_from = stat, values_fill = "0 (0.0%)")
})

cont_rows <- map_dfr(cont_vars, function(var) {
  map_dfr(groups, function(grp) {
    vals <- df1 |> filter(.data[[group_var]] == grp) |> pull(.data[[var]])
    tibble(
      group = grp,
      stat = sprintf("%.2f (%.2f)", median(vals, na.rm = TRUE), mad(vals, na.rm = TRUE))
    )
  }) |>
    pivot_wider(names_from = group, values_from = stat) |>
    mutate(variable = var, level = "Median (MAD)") |>
    dplyr::select(variable, level, everything())
})

desc_table <- bind_rows(cat_rows, cont_rows)

count_row <- df1 |>
  count(ipaq_sf_score_categorical) |>
  pivot_wider(names_from = ipaq_sf_score_categorical, values_from = n) |>
  mutate(variable = "N", level = "") |>
  mutate(across(c(Low, Moderate, High), as.character)) |>
  dplyr::select(variable, level, Low, Moderate, High)

desc_table_with_n <- bind_rows(count_row, desc_table)

desc_table_with_n |>
  knitr::kable(
    col.names = c("Variable", "Level", "Low", "Moderate", "High"),
    caption = "Descriptive Statistics by IPAQ Category",
    align = c("l", "l", "r", "r", "r")
  ) |>
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = FALSE, position = "left"
  ) |>
  pack_rows(index = table(desc_table_with_n$variable)[unique(desc_table_with_n$variable)])
Descriptive Statistics by IPAQ Category
Variable Level Low Moderate High
N
N 3 11 5
current_gender
current_gender 0 1 (33.3%) 8 (72.7%) 0 (0.0%)
current_gender 1 1 (33.3%) 2 (18.2%) 5 (100.0%)
current_gender 2 1 (33.3%) 1 (9.1%) 0 (0.0%)
race
race 0 3 (100.0%) 8 (72.7%) 2 (40.0%)
race 2 0 (0.0%) 1 (9.1%) 0 (0.0%)
race 3 0 (0.0%) 2 (18.2%) 2 (40.0%)
race 1 0 (0.0%) 0 (0.0%) 1 (20.0%)
ethnicity
ethnicity 0 2 (66.7%) 11 (100.0%) 5 (100.0%)
ethnicity 1 1 (33.3%) 0 (0.0%) 0 (0.0%)
sexual_orientation
sexual_orientation 0 2 (66.7%) 10 (90.9%) 3 (60.0%)
sexual_orientation 3 1 (33.3%) 0 (0.0%) 0 (0.0%)
sexual_orientation 2 0 (0.0%) 1 (9.1%) 0 (0.0%)
sexual_orientation 1 0 (0.0%) 0 (0.0%) 2 (40.0%)
education
education 2 2 (66.7%) 4 (36.4%) 2 (40.0%)
education 4 1 (33.3%) 5 (45.5%) 1 (20.0%)
education 1 0 (0.0%) 1 (9.1%) 0 (0.0%)
education 5 0 (0.0%) 1 (9.1%) 0 (0.0%)
education 0 0 (0.0%) 0 (0.0%) 2 (40.0%)
military
military 0 3 (100.0%) 10 (90.9%) 5 (100.0%)
military 1 0 (0.0%) 1 (9.1%) 0 (0.0%)
living_situation
living_situation 0 1 (33.3%) 9 (81.8%) 2 (40.0%)
living_situation 1 1 (33.3%) 2 (18.2%) 0 (0.0%)
living_situation 3 1 (33.3%) 0 (0.0%) 2 (40.0%)
living_situation 2 0 (0.0%) 0 (0.0%) 1 (20.0%)
marital_status
marital_status 0 1 (33.3%) 4 (36.4%) 2 (40.0%)
marital_status 2 2 (66.7%) 5 (45.5%) 2 (40.0%)
marital_status 1 0 (0.0%) 2 (18.2%) 1 (20.0%)
employment_status
employment_status 0 1 (33.3%) 6 (54.5%) 1 (20.0%)
employment_status 3 1 (33.3%) 0 (0.0%) 0 (0.0%)
employment_status 6 1 (33.3%) 0 (0.0%) 2 (40.0%)
employment_status 1 0 (0.0%) 3 (27.3%) 0 (0.0%)
employment_status 2 0 (0.0%) 2 (18.2%) 1 (20.0%)
employment_status 4 0 (0.0%) 0 (0.0%) 1 (20.0%)
personal_income
personal_income 1 2 (66.7%) 2 (18.2%) 0 (0.0%)
personal_income 6 1 (33.3%) 2 (18.2%) 0 (0.0%)
personal_income 4 0 (0.0%) 3 (27.3%) 1 (20.0%)
personal_income 7 0 (0.0%) 1 (9.1%) 1 (20.0%)
personal_income 8 0 (0.0%) 3 (27.3%) 3 (60.0%)
food_insecurity
food_insecurity 1 3 (100.0%) 4 (36.4%) 1 (20.0%)
food_insecurity 0 0 (0.0%) 7 (63.6%) 4 (80.0%)
zarit_burden_score
zarit_burden_score Median (MAD) 23.00 (2.97) 20.00 (4.45) 18.00 (2.97)
ipaq_sf_score_continuous
ipaq_sf_score_continuous Median (MAD) 433.00 (139.36) 1884.00 (531.51) 4653.00 (2154.96)
network_size
network_size Median (MAD) 13.00 (7.41) 12.00 (4.45) 6.00 (2.97)
effective_size
effective_size Median (MAD) 4.71 (2.09) 9.50 (3.10) 5.33 (3.95)
mean_degree
mean_degree Median (MAD) 3.22 (1.39) 4.40 (2.57) 2.89 (3.13)
density
density Median (MAD) 0.38 (0.28) 0.34 (0.16) 0.46 (0.15)
work_hours_per_week
work_hours_per_week Median (MAD) 0.00 (0.00) 36.00 (8.90) 15.00 (14.83)
loneliness
loneliness Median (MAD) 11.00 (2.97) 10.00 (2.97) 6.00 (2.97)
perceived_qol
perceived_qol Median (MAD) 7.00 (0.00) 5.00 (1.48) 8.00 (1.48)
biological_children
biological_children Median (MAD) 0.00 (0.00) 0.00 (0.00) 0.00 (0.00)
not_biological_children
not_biological_children Median (MAD) 0.00 (0.00) 0.00 (0.00) 0.00 (0.00)
contacts
contacts Median (MAD) 300.00 (308.38) 200.00 (74.13) 400.00 (284.66)

Generate per-person table.

df1 |>
  mutate(current_gender = case_when(
    current_gender == 0 ~ "Male",
    current_gender == 1 ~ "Female",
    current_gender == 2 ~ "Non-binary"
  ),
  sexual_orientation = case_when(
   sexual_orientation == 0 ~ "Gay or Lesbian",
   sexual_orientation == 1 ~ "Bisexual",
   sexual_orientation == 2 ~ "Queer",
   sexual_orientation == 3 ~ "Asexual"
  )) |>
  dplyr::select(
    id, current_gender, sexual_orientation,
    zarit_burden_score, ipaq_sf_score_continuous, ipaq_sf_score_categorical,
    contacts, loneliness, perceived_qol,
    network_size, mean_degree, density, effective_size
  ) |>
  arrange(id) |>
  knitr::kable(
    col.names = c(
      "ID", "Gender", "Sexual Orientation",
      "Zarit Burden", "IPAQ (continuous)", "IPAQ (categorical)",
      "Contacts", "Loneliness", "Perceived QoL",
      "Network Size", "Mean Degree", "Density", "Effective Size"
    ),
    caption = "Participant Responses and Social Network Measures",
    align = c("r", "r", "r", "r", "r", "r", "r", "r", "r", "r", "r", "r"),
    digits = 2
  ) |>
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = FALSE,
    position = "left"
  )
Participant Responses and Social Network Measures
ID Gender Sexual Orientation Zarit Burden IPAQ (continuous) IPAQ (categorical) Contacts Loneliness Perceived QoL Network Size Mean Degree Density Effective Size
1 Male Gay or Lesbian 20 1879.5 Moderate 120 12 6 9 5.11 0.64 3.89
2 Male Gay or Lesbian 13 2346.0 Moderate 100 8 4 15 3.47 0.25 11.53
3 Male Gay or Lesbian 21 3306.0 Moderate 150 8 8 21 6.57 0.33 14.43
4 Female Gay or Lesbian 20 13488.0 Moderate 250 10 5 9 6.22 0.78 2.78
5 Male Gay or Lesbian 23 433.0 Low 300 9 7 13 9.69 0.81 3.31
6 Female Bisexual 20 3199.5 High 150 3 9 4 1.50 0.50 2.50
7 Female Gay or Lesbian 14 2646.0 Moderate 1000 14 4 6 2.67 0.53 3.33
8 Female Gay or Lesbian 11 13293.0 High 57 6 7 6 0.67 0.13 5.33
9 Female Gay or Lesbian 18 4833.0 High 400 7 8 6 5.00 1.00 1.00
10 Female Gay or Lesbian 28 0.0 Low 508 13 7 18 3.22 0.19 14.78
11 Male Gay or Lesbian 20 1525.5 Moderate 200 3 8 18 6.89 0.41 11.11
12 Male Gay or Lesbian 21 1272.0 Moderate 247 8 8 11 2.00 0.20 9.00
13 Non-binary Queer 29 1884.0 Moderate 341 12 2 15 4.40 0.31 10.60
14 Male Gay or Lesbian 29 2205.0 Moderate 200 11 6 12 2.50 0.23 9.50
15 Non-binary Asexual 21 527.0 Low 30 11 4 7 2.29 0.38 4.71
16 Female Bisexual 19 1126.5 High 592 4 8 9 2.89 0.36 6.11
17 Male Gay or Lesbian 15 1593.0 Moderate 150 7 5 4 3.00 1.00 1.00
18 Male Gay or Lesbian 23 1539.0 Moderate 5562 12 3 17 5.41 0.34 11.59
19 Female Gay or Lesbian 12 4653.0 High 500 8 7 14 6.00 0.46 8.00

Provide descriptive statistics for categorical variables.

# Create cross tables
categorical_vars <- c(
  "current_gender",
  "race",
  "ethnicity",
  "sexual_orientation",
  "education",
  "military",
  "living_situation",
  "marital_status",
  "employment_status",
  "personal_income",
  "food_insecurity",
  "ipaq_sf_score_categorical"
)

df1 <- df1 |>
  mutate(across(all_of(categorical_vars), as.factor))

cat_tables <- map(categorical_vars, function(v) {
  df1 |>
    count(.data[[v]], name = "n") |>
    mutate(percent = round(100 * n / sum(n), 1), variable = v) |>
    relocate(variable)
})

names(cat_tables) <- categorical_vars

cat_summary <- map_dfr(categorical_vars, function(v) {
  df1 |>
    count(.data[[v]], name = "n") |>
    mutate(percent = round(100 * n / sum(n), 1), variable = v) |>
    rename(category = .data[[v]]) |>
    dplyr::select(variable, category, n, percent)
})

cat_summary |>
  knitr::kable(
    col.names = c("Variable", "Category", "n", "%"),
    caption = "Descriptive Statistics for Categorical Variables",
    align = c("l", "l", "r", "r")
  ) |>
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = FALSE,
    position = "left"
  ) |>
  kableExtra::pack_rows(index = table(cat_summary$variable)[unique(cat_summary$variable)])
Descriptive Statistics for Categorical Variables
Variable Category n %
current_gender
current_gender 0 9 47.4
current_gender 1 8 42.1
current_gender 2 2 10.5
race
race 0 13 68.4
race 1 1 5.3
race 2 1 5.3
race 3 4 21.1
ethnicity
ethnicity 0 18 94.7
ethnicity 1 1 5.3
sexual_orientation
sexual_orientation 0 15 78.9
sexual_orientation 1 2 10.5
sexual_orientation 2 1 5.3
sexual_orientation 3 1 5.3
education
education 0 2 10.5
education 1 1 5.3
education 2 8 42.1
education 4 7 36.8
education 5 1 5.3
military
military 0 18 94.7
military 1 1 5.3
living_situation
living_situation 0 12 63.2
living_situation 1 3 15.8
living_situation 2 1 5.3
living_situation 3 3 15.8
marital_status
marital_status 0 7 36.8
marital_status 1 3 15.8
marital_status 2 9 47.4
employment_status
employment_status 0 8 42.1
employment_status 1 3 15.8
employment_status 2 3 15.8
employment_status 3 1 5.3
employment_status 4 1 5.3
employment_status 6 3 15.8
personal_income
personal_income 1 4 21.1
personal_income 4 4 21.1
personal_income 6 3 15.8
personal_income 7 2 10.5
personal_income 8 6 31.6
food_insecurity
food_insecurity 0 11 57.9
food_insecurity 1 8 42.1
ipaq_sf_score_categorical
ipaq_sf_score_categorical High 5 26.3
ipaq_sf_score_categorical Low 3 15.8
ipaq_sf_score_categorical Moderate 11 57.9 
# Create boxplots
df_cat_long <- df1 |>
  dplyr::select(zarit_burden_score, all_of(categorical_vars)) |>
  pivot_longer(
    cols = all_of(categorical_vars),
    names_to = "variable",
    values_to = "category"
  ) |>
  mutate(category = as.factor(category))

ggplot(df_cat_long, aes(x = category, y = zarit_burden_score)) +
  geom_boxplot() +
  facet_wrap( ~ variable, scales = "free_x") +
  theme_minimal() +
  labs(x = "Category", y = "Caregiver burden") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))