ED_thrive_LMC_Fall_2025_Analysis
K.Evans
2026-01-17
# Path to Fall 2025 data folder
data_dir <- "/Users/krisevans/Dropbox (Personal)/Gradschool/Grad school PC/R Mac working folder/Educators thrive/Labor Management Collaboration (LMC) and Teacher Well-being (TWB)/Fall 2025 Data analysis"
# Excel file name
file_fall <- "Fall 2025 NEA Data.xlsx"
# Full file path
path_fall <- file.path(data_dir, file_fall)
# Load data
survey_raw <- readxl::read_excel(path_fall, sheet = "Survey Data") %>%
janitor::clean_names()
demo_raw <- readxl::read_excel(path_fall, sheet = "Demographics") %>%
janitor::clean_names()
meta_items <- readxl::read_excel(path_fall, sheet = "Metadata") %>%
janitor::clean_names()
# Quick structural checks
glimpse(survey_raw)## Rows: 6,062
## Columns: 52
## $ pid <chr> "bcsssd_f25_1", "bcsssd_f25_2", "bcsssd_f25_3", "bcsssd_f25…
## $ ls1 <dbl> 6, 4, 6, 5, 2, 3, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 4, 5, 3, 6,…
## $ ls2 <dbl> 6, 5, 6, 5, 4, 6, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 5, 5, 4, 6,…
## $ ls3 <dbl> 6, 5, 6, 6, 4, 5, 5, 5, 6, 5, 6, 6, 6, 6, 6, 6, 6, 6, 4, 6,…
## $ ls4 <dbl> 6, 5, 6, 5, 4, 5, 5, 5, 6, 6, 6, 4, 6, 6, 6, 6, 4, 5, 4, 6,…
## $ ls5 <dbl> 6, 5, 6, 5, 3, 5, 5, 4, 6, 6, 6, 5, 6, 6, 6, 6, 5, 6, 3, 5,…
## $ ls6 <dbl> 6, 6, 6, 6, 3, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 5, 6, 4, 6,…
## $ ls7 <dbl> 6, 6, 4, 4, 2, 4, 5, 3, 6, 5, 6, 6, 6, 6, 6, 6, 6, 5, 3, 5,…
## $ ls8 <dbl> 6, 6, NA, 5, 2, 4, 5, 3, 6, 5, 6, 6, 6, 6, 6, 6, 5, 5, 2, 6…
## $ ls9 <dbl> 6, 6, 6, 5, 3, 5, 5, 4, 6, 6, 6, 5, 6, 6, 6, 6, 5, 5, 4, 6,…
## $ ls10 <dbl> 6, 6, 6, 5, 3, 4, 5, 3, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 3, 6,…
## $ ls11 <dbl> 6, 6, 5, 6, 3, 5, 6, 3, 6, 6, 6, 6, 6, 6, 6, 6, 5, 5, 3, 5,…
## $ ls12 <dbl> 6, 6, 6, 5, 4, 6, 5, 4, 6, 6, 6, 6, 6, 6, 5, 6, 6, 6, 4, 6,…
## $ ls13 <dbl> 6, 6, 5, 5, 1, 3, 4, 5, 6, 6, 6, 4, 6, 6, 6, 6, 6, 4, 3, 6,…
## $ ls14 <dbl> 5, 6, 5, 4, 1, 4, 3, 2, 6, 5, 6, 5, 6, 6, 5, 6, 5, 2, 2, 6,…
## $ growth1 <dbl> 5, 6, 6, 5, 4, 6, 5, 4, 6, 6, 6, 5, 6, 6, 6, 6, 6, 6, 6, 5,…
## $ growth2 <dbl> 5, 6, 6, 5, 4, 5, 5, 4, 6, 6, 6, 5, 5, 6, 6, 6, 6, 6, 6, 6,…
## $ growth3 <dbl> 5, 6, 6, 5, 6, 6, 5, 4, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,…
## $ pwb1 <dbl> 6, 3, 5, 5, 5, 3, 5, 6, 5, 4, 5, 6, 5, 6, 6, 6, 4, 5, 6, 4,…
## $ pwb2 <dbl> 6, 6, 6, 5, 5, 4, 5, 6, 6, 5, 6, 6, 5, 6, 6, 6, 5, 5, 6, 4,…
## $ accept1 <dbl> 4, 5, 5, 5, 4, 6, 5, 5, 6, 5, 6, 6, 6, 6, 3, 6, 4, 5, 5, 5,…
## $ accept2 <dbl> 5, 5, 4, 5, 4, 5, 5, 5, 6, 5, 6, 4, 5, 5, 4, 6, 6, 5, 5, 5,…
## $ adapt1 <dbl> 5, 5, 4, 5, 5, 5, 5, 5, 6, 4, 5, 5, 5, 6, 5, 6, 6, 5, 5, 5,…
## $ adapt2 <dbl> 5, 5, 5, 5, 5, 6, 5, 5, 6, 6, 6, 4, 6, 6, 5, 6, 5, 5, 5, 5,…
## $ depl1 <dbl> 5, 6, 2, 3, 4, 4, 2, 4, 5, 1, 1, 3, 4, 4, 2, 2, 6, 6, 4, 3,…
## $ depl2 <dbl> 4, 6, 4, 4, 5, 3, 2, 3, 1, 3, 1, 4, 3, 2, 3, 1, 6, 5, 4, 4,…
## $ fs1 <dbl> 4, 3, 2, 2, 1, 2, 5, 1, 5, 5, 6, 3, 4, 4, 6, 5, 4, 1, 3, 3,…
## $ fs2 <dbl> 5, 3, 4, 4, 1, 4, 4, 4, 6, 5, 6, 3, 6, 6, 6, 6, 3, 1, 3, 5,…
## $ fs3 <dbl> 5, 2, 2, 2, 1, 1, 5, 1, 6, 2, 4, 3, 4, 3, 5, 5, 4, 1, 2, 5,…
## $ fs4 <dbl> 4, 3, 3, 4, 1, 6, 2, 4, 6, 5, 6, 5, 6, 5, 4, 6, 2, 1, 3, 5,…
## $ fs5 <dbl> 4, 5, 2, 4, 1, 4, 2, 5, 6, 4, 6, 3, 6, 6, 4, 6, 2, 1, 3, 5,…
## $ fs6 <dbl> 5, 5, 5, 5, 1, 6, 5, 4, 6, 6, 6, 4, 6, 6, 6, 6, 3, 1, 4, 5,…
## $ ret1 <chr> "5", "5", "5", "5", "5", "5", "3", "5", "5", "5", "5", "5",…
## $ ret2 <chr> NA, "5", "5", "5", "5", "3", "3", "5", "5", "5", "5", "5", …
## $ ret3 <chr> "4", "1", "2", "2", "1", "3", "2", "3", "5", "4", "4", "3",…
## $ job_sat <chr> "8", "8", "9", "8", "4", "5", "7", "9", "10", "10", "10", "…
## $ job_stress <chr> "4", "7", "4", "7", "8", "8", "7", "7", "4", "6", "2", "9",…
## $ npss <chr> "3", "7", "9", "7", "4", "6", "0", "7", "10", "8", "8", "6"…
## $ npsd <chr> "3", "7", "9", "7", "4", "6", "0", "7", "8", "8", "8", "7",…
## $ lmc1 <dbl> 6, 6, 6, 6, 6, 6, 5, 3, 6, 6, 6, 6, 5, 6, 6, 5, 6, 5, 6, 6,…
## $ lmc2 <dbl> 6, 6, 6, 5, 6, 5, 5, 4, 4, 5, 6, 4, 4, 4, 6, 5, 6, 4, 6, 5,…
## $ lmc3 <dbl> 5, 6, 6, 4, 4, 4, 5, 5, 3, 5, 6, 6, 5, 4, 6, 6, 6, 5, 6, 5,…
## $ lmc4 <dbl> 4, 6, 4, 5, 2, 5, 2, 3, 4, 4, 5, 3, 4, 5, 6, 5, 4, 4, 3, 5,…
## $ lmc5 <dbl> 5, 6, 6, 5, 5, 5, 5, 4, 4, 5, 6, 5, 6, 2, 6, 5, 4, 2, 4, 5,…
## $ lmc6 <dbl> 5, 6, 6, 6, 5, 5, 5, 5, 4, 6, 6, 5, 5, 2, 6, 5, 6, 5, 4, 6,…
## $ lmc7 <dbl> 4, 5, 6, 4, 2, 5, 2, 2, 4, 4, 5, 3, 4, 4, 5, 5, 5, 5, 4, 4,…
## $ lmc8 <dbl> 5, 5, 6, 5, 5, NA, 4, 5, 1, 4, 4, 4, 4, 5, 6, 6, 3, 4, 4, 5…
## $ lmc9 <dbl> 5, 4, 5, 5, 4, 6, 4, 4, 1, 4, 4, 3, 5, 5, 6, 5, 4, 2, 3, 4,…
## $ lmc10 <dbl> 5, 4, 4, 4, 2, NA, 4, 3, 1, 4, 4, 4, 4, 5, 6, 5, 4, 1, 3, 5…
## $ lmc11 <dbl> 6, 4, 6, 5, 5, NA, 5, 4, 1, 6, 6, 4, 5, 2, 6, 5, 5, 4, 3, 5…
## $ lmc12 <dbl> 5, 4, 6, 4, 2, NA, 5, 4, 1, 5, 4, 4, 6, 2, 6, 6, 5, 5, 3, 5…
## $ lmc13 <dbl> 6, 6, 6, 6, 2, NA, 5, 4, 1, 5, 4, 6, 5, 4, 6, 5, 6, 2, 3, 6…glimpse(demo_raw)## Rows: 6,062
## Columns: 14
## $ pid <chr> …
## $ district <chr> …
## $ school <chr> …
## $ role <chr> …
## $ gender <chr> …
## $ race <chr> …
## $ ethnicity_hispanic_spanish_latinx <chr> …
## $ years_experience <dbl> …
## $ years_experience_group <chr> …
## $ grade <chr> …
## $ special_education <chr> …
## $ hours_of_collaboration <chr> …
## $ union_member <chr> …
## $ how_much_interaction_do_you_have_with_the_local_education_association_union <chr> …# ---- Sanity checks before joining ----
# Confirm PID exists in both
stopifnot("pid" %in% names(survey_raw))
stopifnot("pid" %in% names(demo_raw))
# Check for duplicate PID rows (duplicates can inflate the dataset after join)
survey_dups <- survey_raw %>% count(pid) %>% filter(n > 1)
demo_dups <- demo_raw %>% count(pid) %>% filter(n > 1)
survey_dups %>% head()## # A tibble: 0 × 2
## # ℹ 2 variables: pid <chr>, n <int>demo_dups %>% head()## # A tibble: 0 × 2
## # ℹ 2 variables: pid <chr>, n <int># ---- Join ----
# Keep all survey respondents, attach demographics when available
df0 <- survey_raw %>%
left_join(demo_raw, by = "pid")
# ---- Post-join checks ----
cat("Rows survey_raw:", nrow(survey_raw), "\n")## Rows survey_raw: 6062cat("Rows df0:", nrow(df0), "\n")## Rows df0: 6062# If df0 has more rows than survey_raw, you likely had duplicates in demo_raw (or survey_raw)
if (nrow(df0) > nrow(survey_raw)) {
warning("Join increased row count. Check duplicate PIDs in survey_raw or demo_raw.")
}
glimpse(df0)## Rows: 6,062
## Columns: 65
## $ pid <chr> …
## $ ls1 <dbl> …
## $ ls2 <dbl> …
## $ ls3 <dbl> …
## $ ls4 <dbl> …
## $ ls5 <dbl> …
## $ ls6 <dbl> …
## $ ls7 <dbl> …
## $ ls8 <dbl> …
## $ ls9 <dbl> …
## $ ls10 <dbl> …
## $ ls11 <dbl> …
## $ ls12 <dbl> …
## $ ls13 <dbl> …
## $ ls14 <dbl> …
## $ growth1 <dbl> …
## $ growth2 <dbl> …
## $ growth3 <dbl> …
## $ pwb1 <dbl> …
## $ pwb2 <dbl> …
## $ accept1 <dbl> …
## $ accept2 <dbl> …
## $ adapt1 <dbl> …
## $ adapt2 <dbl> …
## $ depl1 <dbl> …
## $ depl2 <dbl> …
## $ fs1 <dbl> …
## $ fs2 <dbl> …
## $ fs3 <dbl> …
## $ fs4 <dbl> …
## $ fs5 <dbl> …
## $ fs6 <dbl> …
## $ ret1 <chr> …
## $ ret2 <chr> …
## $ ret3 <chr> …
## $ job_sat <chr> …
## $ job_stress <chr> …
## $ npss <chr> …
## $ npsd <chr> …
## $ lmc1 <dbl> …
## $ lmc2 <dbl> …
## $ lmc3 <dbl> …
## $ lmc4 <dbl> …
## $ lmc5 <dbl> …
## $ lmc6 <dbl> …
## $ lmc7 <dbl> …
## $ lmc8 <dbl> …
## $ lmc9 <dbl> …
## $ lmc10 <dbl> …
## $ lmc11 <dbl> …
## $ lmc12 <dbl> …
## $ lmc13 <dbl> …
## $ district <chr> …
## $ school <chr> …
## $ role <chr> …
## $ gender <chr> …
## $ race <chr> …
## $ ethnicity_hispanic_spanish_latinx <chr> …
## $ years_experience <dbl> …
## $ years_experience_group <chr> …
## $ grade <chr> …
## $ special_education <chr> …
## $ hours_of_collaboration <chr> …
## $ union_member <chr> …
## $ how_much_interaction_do_you_have_with_the_local_education_association_union <chr> …# Outcomes that must be numeric for scoring/modeling
outcome_numeric_vars <- c("ret1","ret2","ret3","job_sat","job_stress","npss","npsd")
# 1) Recode known non-numeric strings (expand this list if you find others)
df0 <- df0 %>%
mutate(
across(
all_of(outcome_numeric_vars),
~ na_if(., "N/A - I plan to retire")
)
)
# 2) Convert to numeric safely:
# - as.character() handles factors too
# - non-numeric values become NA (as they should)
df0 <- df0 %>%
mutate(
across(
all_of(outcome_numeric_vars),
~ suppressWarnings(as.numeric(as.character(.)))
)
)
# 3) Sanity checks: types + NA counts
df0 %>%
select(all_of(outcome_numeric_vars)) %>%
glimpse()## Rows: 6,062
## Columns: 7
## $ ret1 <dbl> 5, 5, 5, 5, 5, 5, 3, 5, 5, 5, 5, 5, 5, 5, NA, 5, 5, 3, NA, …
## $ ret2 <dbl> NA, 5, 5, 5, 5, 3, 3, 5, 5, 5, 5, 5, 5, 5, NA, 5, 3, 3, 5, …
## $ ret3 <dbl> 4, 1, 2, 2, 1, 3, 2, 3, 5, 4, 4, 3, 3, 4, 5, 5, 1, 1, 5, 3,…
## $ job_sat <dbl> 8, 8, 9, 8, 4, 5, 7, 9, 10, 10, 10, 8, 9, 8, 10, 9, 6, 4, 8…
## $ job_stress <dbl> 4, 7, 4, 7, 8, 8, 7, 7, 4, 6, 2, 9, 5, 8, 4, 4, 9, 10, 7, 5…
## $ npss <dbl> 3, 7, 9, 7, 4, 6, 0, 7, 10, 8, 8, 6, 7, 10, 10, 10, 4, 3, 7…
## $ npsd <dbl> 3, 7, 9, 7, 4, 6, 0, 7, 8, 8, 8, 7, 7, 6, 10, 9, 4, 3, 6, 8…df0 %>%
summarise(across(all_of(outcome_numeric_vars), ~ sum(is.na(.)))) %>%
pivot_longer(everything(), names_to = "var", values_to = "n_missing") %>%
arrange(desc(n_missing))## # A tibble: 7 × 2
## var n_missing
## <chr> <int>
## 1 ret1 1002
## 2 ret2 937
## 3 ret3 665
## 4 npss 613
## 5 npsd 602
## 6 job_sat 585
## 7 job_stress 582# Convert common grouping variables to factors (edit names if your columns differ)
df0 <- df0 %>%
mutate(
district = as.factor(district),
school = as.factor(school),
role = as.factor(role),
gender = as.factor(gender),
race = as.factor(race),
grade = as.factor(grade),
union_interaction = as.factor(how_much_interaction_do_you_have_with_the_local_education_association_union)
)
# Ordered version of union involvement (only works if labels match)
df0 <- df0 %>%
mutate(
union_interaction_ord = factor(
union_interaction,
levels = c("None", "A little", "Some", "A lot"),
ordered = TRUE
)
)
# Quick check of levels
df0 %>% count(union_interaction, sort = TRUE)## # A tibble: 5 × 2
## union_interaction n
## <fct> <int>
## 1 A little 2124
## 2 Some 1818
## 3 None 862
## 4 A lot 678
## 5 <NA> 580# -------------------------
# Teacher Well-Being (TWB)
# -------------------------
# Responsive Leadership & Supportive Culture (14 items)
items_twb_leader <- c(
"ls1","ls2","ls3","ls4","ls5","ls6","ls7",
"ls8","ls9","ls10","ls11","ls12","ls13","ls14"
)
# Growth Orientation (3 items)
items_twb_growth <- c("growth1","growth2","growth3")
# Personal Well-Being (2 items)
items_twb_wellbeing <- c("pwb1","pwb2")
# Acceptance (2 items)
items_twb_acceptance <- c("accept1","accept2")
# Adaptability (2 items)
items_twb_adaptability <- c("adapt1","adapt2")
# Depletion (2 items)
items_twb_depletion <- c("depl1","depl2")
# Foundational Supports (6 items)
items_twb_foundational <- c("fs1","fs2","fs3","fs4","fs5","fs6")
# All TWB items (31 total)
items_twb_all <- c(
items_twb_leader,
items_twb_growth,
items_twb_wellbeing,
items_twb_acceptance,
items_twb_adaptability,
items_twb_depletion,
items_twb_foundational
)
length(items_twb_all) # should be 31## [1] 31df <- df0 %>%
mutate(
# -------------------------
# TWB Subscales (means)
# -------------------------
twb_leader_score = rowMeans(
select(., all_of(items_twb_leader)),
na.rm = TRUE
),
growth_score = rowMeans(
select(., all_of(items_twb_growth)),
na.rm = TRUE
),
wellbeing_score = rowMeans(
select(., all_of(items_twb_wellbeing)),
na.rm = TRUE
),
acceptance_score = rowMeans(
select(., all_of(items_twb_acceptance)),
na.rm = TRUE
),
adaptability_score = rowMeans(
select(., all_of(items_twb_adaptability)),
na.rm = TRUE
),
depletion_score = rowMeans(
select(., all_of(items_twb_depletion)),
na.rm = TRUE
),
foundational_support_score = rowMeans(
select(., all_of(items_twb_foundational)),
na.rm = TRUE
),
# -------------------------
# TWB Composite (mean of all 31 items)
# -------------------------
twb_composite_score = rowMeans(
select(., all_of(items_twb_all)),
na.rm = TRUE
)
)
# Quick sanity check
df %>%
summarise(
n = n(),
twb_mean = mean(twb_composite_score, na.rm = TRUE),
leader_mean = mean(twb_leader_score, na.rm = TRUE),
depletion_mean = mean(depletion_score, na.rm = TRUE)
)## # A tibble: 1 × 4
## n twb_mean leader_mean depletion_mean
## <int> <dbl> <dbl> <dbl>
## 1 6062 4.43 4.52 3.73# -------------------------
# Labor-Management Collaboration (LMC)
# -------------------------
items_lmc_union <- c("lmc1","lmc2","lmc3")
items_lmc_school <- c("lmc4","lmc5","lmc6","lmc7")
items_lmc_district <- c("lmc8","lmc9","lmc10","lmc11","lmc12","lmc13")
items_lmc_all <- c(
items_lmc_union,
items_lmc_school,
items_lmc_district
)
df <- df %>%
mutate(
# LMC subscales
union_view_score = rowMeans(
select(., all_of(items_lmc_union)),
na.rm = TRUE
),
school_lmc_score = rowMeans(
select(., all_of(items_lmc_school)),
na.rm = TRUE
),
district_lmc_score = rowMeans(
select(., all_of(items_lmc_district)),
na.rm = TRUE
),
# LMC composite (mean of all 13 items)
lmc_composite_score = rowMeans(
select(., all_of(items_lmc_all)),
na.rm = TRUE
)
)
# Sanity check
df %>%
summarise(
lmc_mean = mean(lmc_composite_score, na.rm = TRUE),
school_mean = mean(school_lmc_score, na.rm = TRUE),
district_mean = mean(district_lmc_score, na.rm = TRUE)
)## # A tibble: 1 × 3
## lmc_mean school_mean district_mean
## <dbl> <dbl> <dbl>
## 1 4.52 4.46 4.42# Retention (mean of 3 items)
items_retention <- c("ret1","ret2","ret3")
df <- df %>%
mutate(
retention_score = rowMeans(
select(., all_of(items_retention)),
na.rm = TRUE
),
# NPS combined (mean of school + district)
nps_score = rowMeans(
select(., npss, npsd),
na.rm = TRUE
)
)
# Check ranges
df %>%
summarise(
retention_min = min(retention_score, na.rm = TRUE),
retention_max = max(retention_score, na.rm = TRUE),
nps_min = min(nps_score, na.rm = TRUE),
nps_max = max(nps_score, na.rm = TRUE)
)## # A tibble: 1 × 4
## retention_min retention_max nps_min nps_max
## <dbl> <dbl> <dbl> <dbl>
## 1 1 5 0 10# ------------------------------------------------------------
# Create item-level dataframes for reliability checks
# (mirrors Spring workflow exactly)
# ------------------------------------------------------------
twb_leader_items <- df %>%
select(all_of(items_twb_leader)) %>%
na.omit()
growth_items <- df %>%
select(all_of(items_twb_growth)) %>%
na.omit()
wellbeing_items <- df %>%
select(all_of(items_twb_wellbeing)) %>%
na.omit()
acceptance_items <- df %>%
select(all_of(items_twb_acceptance)) %>%
na.omit()
adaptability_items <- df %>%
select(all_of(items_twb_adaptability)) %>%
na.omit()
depletion_items <- df %>%
select(all_of(items_twb_depletion)) %>%
na.omit()
foundational_items <- df %>%
select(all_of(items_twb_foundational)) %>%
na.omit()
# LMC items
lmc_union_items <- df %>%
select(all_of(items_lmc_union)) %>%
na.omit()
lmc_school_items <- df %>%
select(all_of(items_lmc_school)) %>%
na.omit()
lmc_district_items <- df %>%
select(all_of(items_lmc_district)) %>%
na.omit()
lmc_all_items <- df %>%
select(all_of(items_lmc_all)) %>%
na.omit()
# ------------------------------------------------------------
# Cronbach's alpha
# ------------------------------------------------------------
psych::alpha(twb_leader_items)##
## Reliability analysis
## Call: psych::alpha(x = twb_leader_items)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.97 0.97 0.97 0.7 33 0.00058 4.5 1.2 0.71
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.97 0.97 0.97
## Duhachek 0.97 0.97 0.97
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## ls1 0.97 0.97 0.97 0.70 31 0.00063 0.0043 0.69
## ls2 0.97 0.97 0.97 0.70 30 0.00066 0.0038 0.69
## ls3 0.97 0.97 0.97 0.70 30 0.00065 0.0039 0.69
## ls4 0.97 0.97 0.97 0.70 30 0.00064 0.0039 0.69
## ls5 0.97 0.97 0.97 0.70 30 0.00065 0.0041 0.69
## ls6 0.97 0.97 0.97 0.70 30 0.00064 0.0042 0.69
## ls7 0.97 0.97 0.97 0.71 31 0.00062 0.0040 0.71
## ls8 0.97 0.97 0.97 0.70 30 0.00065 0.0045 0.69
## ls9 0.97 0.97 0.97 0.70 30 0.00064 0.0043 0.69
## ls10 0.97 0.97 0.97 0.70 31 0.00063 0.0046 0.69
## ls11 0.97 0.97 0.97 0.71 31 0.00062 0.0046 0.72
## ls12 0.97 0.97 0.97 0.72 33 0.00059 0.0032 0.72
## ls13 0.97 0.97 0.97 0.71 32 0.00061 0.0042 0.72
## ls14 0.97 0.97 0.97 0.71 32 0.00060 0.0040 0.72
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## ls1 5426 0.86 0.86 0.85 0.83 4.4 1.5
## ls2 5426 0.90 0.91 0.90 0.89 4.6 1.4
## ls3 5426 0.89 0.89 0.89 0.87 4.6 1.4
## ls4 5426 0.88 0.88 0.88 0.86 4.6 1.4
## ls5 5426 0.89 0.89 0.88 0.87 4.5 1.4
## ls6 5426 0.86 0.87 0.86 0.84 4.8 1.3
## ls7 5426 0.83 0.83 0.82 0.80 4.4 1.4
## ls8 5426 0.88 0.88 0.88 0.86 4.4 1.5
## ls9 5426 0.88 0.88 0.87 0.86 4.8 1.3
## ls10 5426 0.86 0.86 0.85 0.83 4.4 1.5
## ls11 5426 0.83 0.83 0.81 0.80 4.6 1.3
## ls12 5426 0.76 0.76 0.73 0.72 4.9 1.4
## ls13 5426 0.80 0.80 0.77 0.76 4.2 1.4
## ls14 5426 0.79 0.79 0.77 0.75 4.3 1.4
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 miss
## ls1 0.07 0.07 0.09 0.17 0.35 0.25 0
## ls2 0.05 0.06 0.07 0.17 0.33 0.31 0
## ls3 0.04 0.06 0.07 0.19 0.35 0.28 0
## ls4 0.05 0.05 0.08 0.20 0.36 0.26 0
## ls5 0.06 0.06 0.07 0.21 0.35 0.25 0
## ls6 0.04 0.05 0.06 0.17 0.36 0.33 0
## ls7 0.06 0.07 0.10 0.18 0.36 0.23 0
## ls8 0.07 0.07 0.08 0.20 0.32 0.26 0
## ls9 0.04 0.04 0.06 0.16 0.35 0.34 0
## ls10 0.07 0.07 0.09 0.19 0.30 0.29 0
## ls11 0.05 0.05 0.07 0.20 0.35 0.28 0
## ls12 0.05 0.05 0.04 0.13 0.34 0.40 0
## ls13 0.07 0.08 0.11 0.24 0.33 0.18 0
## ls14 0.06 0.07 0.11 0.23 0.34 0.19 0psych::alpha(growth_items)##
## Reliability analysis
## Call: psych::alpha(x = growth_items)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.81 0.81 0.75 0.59 4.4 0.0043 5.2 0.73 0.58
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.8 0.81 0.82
## Duhachek 0.8 0.81 0.82
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## growth1 0.73 0.74 0.58 0.58 2.8 0.0070 NA 0.58
## growth2 0.72 0.73 0.57 0.57 2.7 0.0073 NA 0.57
## growth3 0.76 0.77 0.62 0.62 3.3 0.0062 NA 0.62
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## growth1 5682 0.87 0.86 0.74 0.67 5.1 0.95
## growth2 5682 0.86 0.86 0.76 0.68 5.1 0.86
## growth3 5682 0.82 0.84 0.71 0.64 5.4 0.76
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 miss
## growth1 0.01 0.01 0.03 0.15 0.43 0.37 0
## growth2 0.01 0.01 0.02 0.14 0.46 0.37 0
## growth3 0.00 0.00 0.01 0.08 0.43 0.48 0psych::alpha(wellbeing_items)##
## Reliability analysis
## Call: psych::alpha(x = wellbeing_items)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.75 0.75 0.6 0.6 3 0.0065 5 0.94 0.6
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.74 0.75 0.77
## Duhachek 0.74 0.75 0.77
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## pwb1 0.62 0.6 0.36 0.6 1.5 NA 0 0.6
## pwb2 0.58 0.6 0.36 0.6 1.5 NA 0 0.6
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## pwb1 5706 0.90 0.9 0.7 0.6 5.0 1.1
## pwb2 5706 0.89 0.9 0.7 0.6 5.1 1.0
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 miss
## pwb1 0.01 0.03 0.05 0.17 0.38 0.37 0
## pwb2 0.01 0.02 0.04 0.12 0.39 0.42 0psych::alpha(acceptance_items)##
## Reliability analysis
## Call: psych::alpha(x = acceptance_items)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.73 0.73 0.57 0.57 2.7 0.0073 4.8 0.83 0.57
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.71 0.73 0.74
## Duhachek 0.71 0.73 0.74
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## accept1 0.57 0.57 0.33 0.57 1.3 NA 0 0.57
## accept2 0.57 0.57 0.33 0.57 1.3 NA 0 0.57
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## accept1 5647 0.89 0.89 0.67 0.57 4.8 0.94
## accept2 5647 0.89 0.89 0.67 0.57 4.7 0.93
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 miss
## accept1 0.01 0.02 0.05 0.23 0.47 0.23 0
## accept2 0.01 0.02 0.06 0.24 0.49 0.18 0psych::alpha(adaptability_items)##
## Reliability analysis
## Call: psych::alpha(x = adaptability_items)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.67 0.68 0.51 0.51 2.1 0.0086 5.1 0.65 0.51
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.66 0.67 0.69
## Duhachek 0.66 0.67 0.69
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## adapt1 0.55 0.51 0.26 0.51 1 NA 0 0.51
## adapt2 0.47 0.51 0.26 0.51 1 NA 0 0.51
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## adapt1 5635 0.88 0.87 0.62 0.51 5.1 0.78
## adapt2 5635 0.86 0.87 0.62 0.51 5.2 0.72
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 miss
## adapt1 0 0.01 0.02 0.13 0.54 0.29 0
## adapt2 0 0.00 0.01 0.12 0.54 0.33 0psych::alpha(depletion_items)##
## Reliability analysis
## Call: psych::alpha(x = depletion_items)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.74 0.74 0.59 0.59 2.8 0.007 3.7 1.2 0.59
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.72 0.74 0.75
## Duhachek 0.73 0.74 0.75
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## depl1 0.58 0.59 0.34 0.59 1.4 NA 0 0.59
## depl2 0.59 0.59 0.34 0.59 1.4 NA 0 0.59
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## depl1 5596 0.89 0.89 0.68 0.59 4.0 1.4
## depl2 5596 0.89 0.89 0.68 0.59 3.4 1.4
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 miss
## depl1 0.04 0.14 0.12 0.33 0.21 0.17 0
## depl2 0.07 0.25 0.16 0.29 0.16 0.07 0psych::alpha(foundational_items)##
## Reliability analysis
## Call: psych::alpha(x = foundational_items)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.84 0.84 0.83 0.47 5.2 0.0034 3.5 1.1 0.43
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.83 0.84 0.84
## Duhachek 0.83 0.84 0.84
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## fs1 0.84 0.84 0.82 0.51 5.2 0.0035 0.0092 0.51
## fs2 0.81 0.82 0.80 0.47 4.4 0.0040 0.0135 0.45
## fs3 0.82 0.82 0.80 0.48 4.6 0.0039 0.0146 0.48
## fs4 0.78 0.79 0.75 0.42 3.7 0.0046 0.0044 0.42
## fs5 0.79 0.79 0.76 0.43 3.8 0.0045 0.0050 0.43
## fs6 0.82 0.82 0.80 0.47 4.5 0.0040 0.0130 0.45
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## fs1 5456 0.65 0.64 0.52 0.48 3.0 1.5
## fs2 5456 0.73 0.73 0.65 0.60 3.5 1.4
## fs3 5456 0.72 0.71 0.62 0.57 3.1 1.5
## fs4 5456 0.83 0.84 0.83 0.74 3.4 1.4
## fs5 5456 0.81 0.82 0.80 0.71 3.6 1.4
## fs6 5456 0.72 0.73 0.65 0.59 4.2 1.4
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 miss
## fs1 0.24 0.20 0.16 0.19 0.18 0.03 0
## fs2 0.13 0.16 0.16 0.28 0.21 0.06 0
## fs3 0.19 0.20 0.17 0.20 0.18 0.05 0
## fs4 0.13 0.16 0.17 0.28 0.20 0.06 0
## fs5 0.11 0.15 0.17 0.30 0.23 0.05 0
## fs6 0.07 0.08 0.11 0.22 0.41 0.12 0# Full TWB scale (all 31 items)
psych::alpha(df %>% select(all_of(items_twb_all)) %>% na.omit())## Some items ( depl1 depl2 ) were negatively correlated with the first principal component and
## probably should be reversed.
## To do this, run the function again with the 'check.keys=TRUE' option##
## Reliability analysis
## Call: psych::alpha(x = df %>% select(all_of(items_twb_all)) %>% na.omit())
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.93 0.93 0.96 0.3 13 0.0012 4.4 0.74 0.25
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.93 0.93 0.94
## Duhachek 0.93 0.93 0.94
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## ls1 0.93 0.92 0.95 0.29 12 0.0013 0.068 0.25
## ls2 0.93 0.92 0.95 0.29 12 0.0013 0.067 0.25
## ls3 0.93 0.92 0.95 0.29 12 0.0013 0.067 0.25
## ls4 0.93 0.92 0.95 0.29 12 0.0013 0.068 0.25
## ls5 0.93 0.92 0.95 0.29 12 0.0013 0.067 0.25
## ls6 0.93 0.92 0.95 0.29 12 0.0013 0.068 0.25
## ls7 0.93 0.93 0.95 0.29 12 0.0013 0.069 0.25
## ls8 0.93 0.92 0.95 0.29 12 0.0013 0.068 0.25
## ls9 0.93 0.92 0.95 0.29 12 0.0013 0.068 0.25
## ls10 0.93 0.93 0.95 0.29 12 0.0013 0.068 0.25
## ls11 0.93 0.93 0.95 0.29 12 0.0013 0.069 0.25
## ls12 0.93 0.93 0.95 0.29 13 0.0012 0.070 0.25
## ls13 0.93 0.93 0.95 0.29 12 0.0013 0.069 0.25
## ls14 0.93 0.93 0.95 0.29 12 0.0012 0.070 0.25
## growth1 0.93 0.93 0.96 0.30 13 0.0012 0.073 0.24
## growth2 0.93 0.93 0.96 0.30 13 0.0012 0.073 0.25
## growth3 0.93 0.93 0.96 0.31 13 0.0012 0.073 0.26
## pwb1 0.93 0.93 0.96 0.31 14 0.0011 0.071 0.26
## pwb2 0.93 0.93 0.96 0.31 14 0.0011 0.071 0.26
## accept1 0.93 0.93 0.96 0.31 13 0.0012 0.072 0.26
## accept2 0.93 0.93 0.96 0.30 13 0.0012 0.073 0.26
## adapt1 0.93 0.93 0.96 0.31 13 0.0012 0.073 0.26
## adapt2 0.93 0.93 0.96 0.31 13 0.0012 0.073 0.26
## depl1 0.94 0.94 0.96 0.33 15 0.0011 0.058 0.26
## depl2 0.94 0.94 0.96 0.33 15 0.0011 0.060 0.26
## fs1 0.93 0.93 0.96 0.31 13 0.0011 0.073 0.26
## fs2 0.93 0.93 0.96 0.30 13 0.0012 0.073 0.25
## fs3 0.93 0.93 0.96 0.30 13 0.0012 0.073 0.25
## fs4 0.93 0.93 0.95 0.30 13 0.0012 0.072 0.25
## fs5 0.93 0.93 0.95 0.30 13 0.0012 0.072 0.25
## fs6 0.93 0.93 0.96 0.30 13 0.0012 0.073 0.25
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## ls1 5076 0.81 0.79 0.79 0.79 4.4 1.47
## ls2 5076 0.85 0.83 0.84 0.83 4.6 1.40
## ls3 5076 0.84 0.82 0.83 0.82 4.6 1.36
## ls4 5076 0.83 0.81 0.81 0.81 4.5 1.36
## ls5 5076 0.84 0.81 0.82 0.82 4.5 1.40
## ls6 5076 0.81 0.79 0.79 0.79 4.7 1.31
## ls7 5076 0.79 0.77 0.77 0.77 4.4 1.43
## ls8 5076 0.83 0.81 0.81 0.81 4.4 1.48
## ls9 5076 0.83 0.81 0.81 0.81 4.8 1.32
## ls10 5076 0.81 0.78 0.78 0.78 4.4 1.50
## ls11 5076 0.79 0.77 0.77 0.76 4.6 1.35
## ls12 5076 0.71 0.70 0.69 0.68 4.8 1.38
## ls13 5076 0.79 0.77 0.77 0.76 4.2 1.42
## ls14 5076 0.76 0.74 0.74 0.73 4.3 1.41
## growth1 5076 0.48 0.53 0.51 0.44 5.1 0.95
## growth2 5076 0.42 0.48 0.46 0.39 5.1 0.85
## growth3 5076 0.36 0.42 0.40 0.33 5.4 0.76
## pwb1 5076 0.23 0.27 0.24 0.18 5.0 1.07
## pwb2 5076 0.25 0.29 0.27 0.20 5.1 1.03
## accept1 5076 0.37 0.42 0.39 0.33 4.8 0.94
## accept2 5076 0.42 0.48 0.46 0.39 4.7 0.93
## adapt1 5076 0.36 0.42 0.39 0.33 5.1 0.78
## adapt2 5076 0.36 0.42 0.40 0.33 5.2 0.72
## depl1 5076 -0.18 -0.19 -0.23 -0.23 4.0 1.38
## depl2 5076 -0.12 -0.15 -0.18 -0.18 3.4 1.39
## fs1 5076 0.41 0.40 0.36 0.35 3.0 1.53
## fs2 5076 0.54 0.53 0.51 0.49 3.5 1.44
## fs3 5076 0.51 0.50 0.48 0.46 3.1 1.54
## fs4 5076 0.65 0.64 0.63 0.61 3.4 1.45
## fs5 5076 0.58 0.57 0.56 0.54 3.6 1.40
## fs6 5076 0.59 0.58 0.56 0.55 4.2 1.36
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 miss
## ls1 0.07 0.07 0.09 0.17 0.36 0.24 0
## ls2 0.05 0.06 0.07 0.18 0.33 0.31 0
## ls3 0.05 0.06 0.07 0.19 0.35 0.27 0
## ls4 0.05 0.05 0.08 0.20 0.36 0.26 0
## ls5 0.06 0.05 0.07 0.21 0.35 0.25 0
## ls6 0.04 0.05 0.06 0.17 0.36 0.33 0
## ls7 0.06 0.07 0.10 0.18 0.36 0.23 0
## ls8 0.07 0.07 0.08 0.20 0.32 0.26 0
## ls9 0.04 0.04 0.06 0.16 0.35 0.34 0
## ls10 0.07 0.07 0.09 0.19 0.30 0.29 0
## ls11 0.05 0.05 0.07 0.20 0.36 0.27 0
## ls12 0.05 0.05 0.04 0.12 0.34 0.40 0
## ls13 0.07 0.08 0.11 0.24 0.33 0.18 0
## ls14 0.07 0.07 0.11 0.23 0.34 0.19 0
## growth1 0.01 0.01 0.03 0.15 0.43 0.37 0
## growth2 0.01 0.01 0.02 0.14 0.46 0.37 0
## growth3 0.01 0.00 0.01 0.08 0.43 0.48 0
## pwb1 0.01 0.03 0.05 0.17 0.38 0.37 0
## pwb2 0.01 0.02 0.04 0.12 0.39 0.42 0
## accept1 0.01 0.02 0.06 0.23 0.47 0.22 0
## accept2 0.01 0.02 0.06 0.24 0.50 0.17 0
## adapt1 0.00 0.01 0.02 0.13 0.55 0.29 0
## adapt2 0.00 0.00 0.01 0.12 0.54 0.33 0
## depl1 0.04 0.14 0.12 0.33 0.21 0.17 0
## depl2 0.07 0.25 0.15 0.29 0.16 0.07 0
## fs1 0.24 0.20 0.16 0.19 0.18 0.03 0
## fs2 0.13 0.16 0.16 0.28 0.21 0.06 0
## fs3 0.20 0.20 0.18 0.20 0.18 0.05 0
## fs4 0.13 0.16 0.17 0.28 0.20 0.06 0
## fs5 0.11 0.15 0.17 0.29 0.23 0.05 0
## fs6 0.07 0.08 0.10 0.22 0.42 0.12 0# LMC alphas
psych::alpha(lmc_union_items)##
## Reliability analysis
## Call: psych::alpha(x = lmc_union_items)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.75 0.74 0.71 0.49 2.9 0.0058 4.8 1 0.4
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.74 0.75 0.76
## Duhachek 0.74 0.75 0.76
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## lmc1 0.85 0.85 0.74 0.74 5.8 0.004 NA 0.74
## lmc2 0.50 0.50 0.33 0.33 1.0 0.013 NA 0.33
## lmc3 0.57 0.57 0.40 0.40 1.3 0.012 NA 0.40
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## lmc1 5409 0.68 0.71 0.44 0.39 5.3 1.1
## lmc2 5409 0.89 0.88 0.84 0.72 4.6 1.3
## lmc3 5409 0.87 0.85 0.79 0.66 4.5 1.3
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 miss
## lmc1 0.03 0.02 0.02 0.07 0.30 0.58 0
## lmc2 0.04 0.05 0.06 0.25 0.34 0.25 0
## lmc3 0.04 0.04 0.07 0.25 0.34 0.25 0psych::alpha(lmc_school_items)##
## Reliability analysis
## Call: psych::alpha(x = lmc_school_items)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.84 0.84 0.81 0.57 5.4 0.0036 4.5 0.98 0.56
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.84 0.84 0.85
## Duhachek 0.84 0.84 0.85
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## lmc4 0.78 0.78 0.71 0.54 3.5 0.0053 0.0057 0.52
## lmc5 0.80 0.80 0.75 0.58 4.1 0.0048 0.0183 0.52
## lmc6 0.85 0.85 0.80 0.65 5.7 0.0036 0.0044 0.62
## lmc7 0.77 0.77 0.69 0.52 3.3 0.0056 0.0058 0.48
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## lmc4 5239 0.85 0.86 0.81 0.73 4.5 1.2
## lmc5 5239 0.82 0.82 0.73 0.67 4.5 1.2
## lmc6 5239 0.75 0.75 0.60 0.56 4.7 1.2
## lmc7 5239 0.87 0.87 0.83 0.76 4.2 1.2
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 miss
## lmc4 0.03 0.04 0.08 0.29 0.38 0.17 0
## lmc5 0.03 0.06 0.08 0.25 0.41 0.18 0
## lmc6 0.03 0.04 0.06 0.18 0.46 0.23 0
## lmc7 0.04 0.05 0.11 0.35 0.34 0.11 0psych::alpha(lmc_district_items)##
## Reliability analysis
## Call: psych::alpha(x = lmc_district_items)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.89 0.89 0.89 0.58 8.1 0.0024 4.4 0.9 0.54
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.88 0.89 0.89
## Duhachek 0.88 0.89 0.89
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## lmc8 0.86 0.87 0.86 0.56 6.4 0.0031 0.0108 0.53
## lmc9 0.86 0.86 0.85 0.56 6.2 0.0031 0.0082 0.53
## lmc10 0.86 0.86 0.84 0.55 6.1 0.0032 0.0064 0.53
## lmc11 0.88 0.89 0.88 0.61 7.9 0.0026 0.0130 0.63
## lmc12 0.86 0.87 0.86 0.56 6.4 0.0031 0.0159 0.48
## lmc13 0.89 0.89 0.88 0.61 7.9 0.0025 0.0131 0.63
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## lmc8 5172 0.83 0.83 0.80 0.74 4.4 1.1
## lmc9 5172 0.84 0.85 0.83 0.76 4.4 1.0
## lmc10 5172 0.86 0.86 0.85 0.78 4.3 1.1
## lmc11 5172 0.73 0.73 0.64 0.61 4.5 1.1
## lmc12 5172 0.83 0.83 0.79 0.75 4.4 1.2
## lmc13 5172 0.74 0.73 0.64 0.61 4.5 1.2
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 miss
## lmc8 0.03 0.04 0.09 0.32 0.41 0.13 0
## lmc9 0.02 0.03 0.08 0.34 0.41 0.12 0
## lmc10 0.03 0.04 0.10 0.38 0.36 0.09 0
## lmc11 0.03 0.04 0.07 0.25 0.46 0.15 0
## lmc12 0.03 0.04 0.09 0.32 0.39 0.13 0
## lmc13 0.03 0.05 0.08 0.24 0.41 0.19 0psych::alpha(lmc_all_items)##
## Reliability analysis
## Call: psych::alpha(x = lmc_all_items)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.91 0.91 0.93 0.45 11 0.0018 4.5 0.81 0.44
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.91 0.91 0.92
## Duhachek 0.91 0.91 0.92
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## lmc1 0.92 0.92 0.93 0.48 11.1 0.0018 0.012 0.47
## lmc2 0.91 0.91 0.92 0.45 9.9 0.0019 0.018 0.47
## lmc3 0.91 0.91 0.92 0.45 9.6 0.0020 0.018 0.44
## lmc4 0.90 0.90 0.92 0.44 9.5 0.0020 0.018 0.43
## lmc5 0.91 0.91 0.93 0.46 10.1 0.0019 0.017 0.46
## lmc6 0.91 0.91 0.92 0.45 9.9 0.0019 0.019 0.44
## lmc7 0.90 0.90 0.92 0.44 9.4 0.0020 0.017 0.43
## lmc8 0.90 0.90 0.92 0.44 9.4 0.0020 0.016 0.43
## lmc9 0.90 0.90 0.92 0.44 9.5 0.0020 0.016 0.44
## lmc10 0.90 0.90 0.92 0.44 9.4 0.0020 0.015 0.43
## lmc11 0.91 0.91 0.92 0.45 9.7 0.0019 0.019 0.43
## lmc12 0.90 0.90 0.92 0.44 9.3 0.0020 0.017 0.43
## lmc13 0.91 0.91 0.93 0.45 9.8 0.0019 0.020 0.44
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## lmc1 5101 0.46 0.46 0.38 0.37 5.3 1.1
## lmc2 5101 0.69 0.68 0.65 0.61 4.6 1.3
## lmc3 5101 0.73 0.72 0.70 0.66 4.6 1.3
## lmc4 5101 0.75 0.75 0.73 0.70 4.5 1.2
## lmc5 5101 0.64 0.63 0.60 0.56 4.5 1.2
## lmc6 5101 0.67 0.66 0.63 0.60 4.7 1.2
## lmc7 5101 0.76 0.76 0.74 0.70 4.2 1.2
## lmc8 5101 0.76 0.76 0.75 0.71 4.4 1.1
## lmc9 5101 0.74 0.75 0.74 0.69 4.4 1.0
## lmc10 5101 0.75 0.76 0.76 0.70 4.3 1.1
## lmc11 5101 0.70 0.71 0.68 0.64 4.5 1.1
## lmc12 5101 0.77 0.77 0.76 0.72 4.4 1.2
## lmc13 5101 0.69 0.69 0.65 0.62 4.5 1.2
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 miss
## lmc1 0.02 0.02 0.01 0.06 0.29 0.59 0
## lmc2 0.04 0.04 0.06 0.25 0.35 0.26 0
## lmc3 0.04 0.04 0.07 0.25 0.34 0.25 0
## lmc4 0.03 0.04 0.08 0.29 0.39 0.17 0
## lmc5 0.03 0.06 0.08 0.24 0.41 0.18 0
## lmc6 0.03 0.04 0.06 0.17 0.46 0.24 0
## lmc7 0.04 0.05 0.11 0.34 0.34 0.11 0
## lmc8 0.03 0.04 0.09 0.32 0.40 0.12 0
## lmc9 0.02 0.03 0.08 0.34 0.41 0.12 0
## lmc10 0.03 0.04 0.10 0.38 0.36 0.09 0
## lmc11 0.03 0.04 0.07 0.25 0.46 0.16 0
## lmc12 0.03 0.05 0.09 0.32 0.39 0.13 0
## lmc13 0.03 0.05 0.08 0.24 0.42 0.19 0# ============================================================
# CFA: TWB 7-factor model (Fall 2025)
# Purpose:
# 1) Validate the hypothesized 7-factor structure
# 2) Get item-level diagnostics (loadings, residuals) to shorten leadership
# ============================================================
# 7-factor TWB model specification (mirrors item sets)
mod_twb_7 <- '
leader =~ ls1 + ls2 + ls3 + ls4 + ls5 + ls6 + ls7 + ls8 + ls9 + ls10 + ls11 + ls12 + ls13 + ls14
growth =~ growth1 + growth2 + growth3
wellbeing =~ pwb1 + pwb2
acceptance =~ accept1 + accept2
adaptability =~ adapt1 + adapt2
depletion =~ depl1 + depl2
foundational =~ fs1 + fs2 + fs3 + fs4 + fs5 + fs6
'
fit_twb_7 <- lavaan::cfa(
model = mod_twb_7,
data = df,
estimator = "MLR", # robust SEs + fit (good default for Likert-ish data)
missing = "fiml", # uses all available data under MAR assumption
std.lv = TRUE # sets latent var variance = 1 (stable standardized loadings)
)
# Core output
summary(fit_twb_7, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6-19 ended normally after 113 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 114
##
## Used Total
## Number of observations 5987 6062
## Number of missing patterns 222
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 6004.780 4509.157
## Degrees of freedom 413 413
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.332
## Yuan-Bentler correction (Mplus variant)
##
## Model Test Baseline Model:
##
## Test statistic 123134.053 89115.378
## Degrees of freedom 465 465
## P-value 0.000 0.000
## Scaling correction factor 1.382
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.954 0.954
## Tucker-Lewis Index (TLI) 0.949 0.948
##
## Robust Comparative Fit Index (CFI) 0.955
## Robust Tucker-Lewis Index (TLI) 0.949
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -228179.389 -228179.389
## Scaling correction factor 1.455
## for the MLR correction
## Loglikelihood unrestricted model (H1) -225176.999 -225176.999
## Scaling correction factor 1.358
## for the MLR correction
##
## Akaike (AIC) 456586.779 456586.779
## Bayesian (BIC) 457350.276 457350.276
## Sample-size adjusted Bayesian (SABIC) 456988.016 456988.016
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.048 0.041
## 90 Percent confidence interval - lower 0.046 0.040
## 90 Percent confidence interval - upper 0.049 0.042
## P-value H_0: RMSEA <= 0.050 1.000 1.000
## P-value H_0: RMSEA >= 0.080 0.000 0.000
##
## Robust RMSEA 0.049
## 90 Percent confidence interval - lower 0.047
## 90 Percent confidence interval - upper 0.050
## P-value H_0: Robust RMSEA <= 0.050 0.967
## P-value H_0: Robust RMSEA >= 0.080 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.032 0.032
##
## Parameter Estimates:
##
## Standard errors Sandwich
## Information bread Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## leader =~
## ls1 1.240 0.016 78.269 0.000 1.240 0.847
## ls2 1.264 0.015 81.926 0.000 1.264 0.906
## ls3 1.205 0.015 78.106 0.000 1.205 0.890
## ls4 1.188 0.016 75.952 0.000 1.188 0.880
## ls5 1.233 0.016 79.550 0.000 1.233 0.887
## ls6 1.123 0.017 68.032 0.000 1.123 0.861
## ls7 1.165 0.016 72.212 0.000 1.165 0.818
## ls8 1.282 0.015 83.365 0.000 1.282 0.868
## ls9 1.140 0.017 67.329 0.000 1.140 0.871
## ls10 1.260 0.016 79.595 0.000 1.260 0.840
## ls11 1.076 0.017 62.958 0.000 1.076 0.798
## ls12 0.995 0.020 49.678 0.000 0.995 0.726
## ls13 1.086 0.016 67.649 0.000 1.086 0.770
## ls14 1.064 0.017 62.446 0.000 1.064 0.756
## growth =~
## growth1 0.752 0.016 46.462 0.000 0.752 0.792
## growth2 0.679 0.016 43.019 0.000 0.679 0.790
## growth3 0.553 0.017 33.465 0.000 0.553 0.729
## wellbeing =~
## pwb1 0.815 0.019 43.175 0.000 0.815 0.762
## pwb2 0.817 0.019 42.195 0.000 0.817 0.791
## acceptance =~
## accept1 0.671 0.015 44.702 0.000 0.671 0.713
## accept2 0.747 0.014 52.425 0.000 0.747 0.800
## adaptability =~
## adapt1 0.559 0.014 38.701 0.000 0.559 0.721
## adapt2 0.509 0.014 36.718 0.000 0.509 0.708
## depletion =~
## depl1 1.066 0.021 49.677 0.000 1.066 0.770
## depl2 1.058 0.021 50.537 0.000 1.058 0.761
## foundational =~
## fs1 0.774 0.019 39.870 0.000 0.774 0.503
## fs2 0.942 0.017 55.011 0.000 0.942 0.652
## fs3 0.931 0.018 51.805 0.000 0.931 0.605
## fs4 1.246 0.013 92.837 0.000 1.246 0.861
## fs5 1.149 0.015 78.909 0.000 1.149 0.819
## fs6 0.901 0.018 50.181 0.000 0.901 0.659
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## leader ~~
## growth 0.376 0.016 23.908 0.000 0.376 0.376
## wellbeing 0.149 0.017 8.604 0.000 0.149 0.149
## acceptance 0.342 0.017 20.357 0.000 0.342 0.342
## adaptability 0.332 0.017 19.317 0.000 0.332 0.332
## depletion -0.154 0.017 -9.104 0.000 -0.154 -0.154
## foundational 0.591 0.012 51.234 0.000 0.591 0.591
## growth ~~
## wellbeing 0.438 0.020 21.830 0.000 0.438 0.438
## acceptance 0.491 0.019 26.275 0.000 0.491 0.491
## adaptability 0.578 0.019 29.819 0.000 0.578 0.578
## depletion -0.219 0.019 -11.650 0.000 -0.219 -0.219
## foundational 0.374 0.015 24.216 0.000 0.374 0.374
## wellbeing ~~
## acceptance 0.432 0.020 21.751 0.000 0.432 0.432
## adaptability 0.446 0.022 20.710 0.000 0.446 0.446
## depletion -0.528 0.016 -33.132 0.000 -0.528 -0.528
## foundational 0.229 0.017 13.804 0.000 0.229 0.229
## acceptance ~~
## adaptability 0.925 0.014 63.831 0.000 0.925 0.925
## depletion -0.402 0.018 -22.032 0.000 -0.402 -0.402
## foundational 0.381 0.016 23.474 0.000 0.381 0.381
## adaptability ~~
## depletion -0.341 0.020 -16.979 0.000 -0.341 -0.341
## foundational 0.348 0.017 20.476 0.000 0.348 0.348
## depletion ~~
## foundational -0.339 0.017 -19.549 0.000 -0.339 -0.339
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .ls1 4.409 0.019 231.030 0.000 4.409 3.012
## .ls2 4.619 0.018 253.971 0.000 4.619 3.310
## .ls3 4.574 0.018 259.155 0.000 4.574 3.380
## .ls4 4.547 0.018 258.305 0.000 4.547 3.367
## .ls5 4.480 0.018 247.080 0.000 4.480 3.221
## .ls6 4.740 0.017 278.604 0.000 4.740 3.633
## .ls7 4.406 0.019 236.661 0.000 4.406 3.092
## .ls8 4.400 0.019 227.410 0.000 4.400 2.979
## .ls9 4.762 0.017 277.316 0.000 4.762 3.635
## .ls10 4.429 0.020 225.184 0.000 4.429 2.953
## .ls11 4.573 0.018 257.963 0.000 4.573 3.391
## .ls12 4.840 0.018 267.527 0.000 4.840 3.529
## .ls13 4.210 0.019 227.474 0.000 4.210 2.986
## .ls14 4.256 0.019 229.786 0.000 4.256 3.022
## .growth1 5.083 0.013 406.114 0.000 5.083 5.349
## .growth2 5.148 0.011 454.155 0.000 5.148 5.990
## .growth3 5.354 0.010 534.623 0.000 5.354 7.063
## .pwb1 4.983 0.014 352.924 0.000 4.983 4.661
## .pwb2 5.110 0.014 374.566 0.000 5.110 4.947
## .accept1 4.812 0.012 385.697 0.000 4.812 5.118
## .accept2 4.723 0.012 381.582 0.000 4.723 5.061
## .adapt1 5.089 0.010 494.769 0.000 5.089 6.566
## .adapt2 5.173 0.010 542.767 0.000 5.173 7.196
## .depl1 4.024 0.018 218.263 0.000 4.024 2.906
## .depl2 3.427 0.019 185.110 0.000 3.427 2.466
## .fs1 2.985 0.021 145.244 0.000 2.985 1.940
## .fs2 3.459 0.019 179.322 0.000 3.459 2.391
## .fs3 3.146 0.021 153.178 0.000 3.146 2.044
## .fs4 3.453 0.019 179.370 0.000 3.453 2.384
## .fs5 3.569 0.019 190.904 0.000 3.569 2.543
## .fs6 4.197 0.018 231.055 0.000 4.197 3.072
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .ls1 0.605 0.019 31.500 0.000 0.605 0.282
## .ls2 0.348 0.013 27.522 0.000 0.348 0.179
## .ls3 0.380 0.013 29.609 0.000 0.380 0.208
## .ls4 0.412 0.014 29.550 0.000 0.412 0.226
## .ls5 0.414 0.015 27.821 0.000 0.414 0.214
## .ls6 0.440 0.014 32.197 0.000 0.440 0.259
## .ls7 0.672 0.019 35.182 0.000 0.672 0.331
## .ls8 0.540 0.016 34.767 0.000 0.540 0.247
## .ls9 0.416 0.013 32.987 0.000 0.416 0.242
## .ls10 0.663 0.019 35.492 0.000 0.663 0.295
## .ls11 0.662 0.018 36.745 0.000 0.662 0.364
## .ls12 0.890 0.024 37.132 0.000 0.890 0.473
## .ls13 0.809 0.020 40.350 0.000 0.809 0.407
## .ls14 0.850 0.023 36.696 0.000 0.850 0.429
## .growth1 0.337 0.019 17.742 0.000 0.337 0.373
## .growth2 0.278 0.014 20.187 0.000 0.278 0.376
## .growth3 0.269 0.012 22.951 0.000 0.269 0.469
## .pwb1 0.478 0.024 20.072 0.000 0.478 0.419
## .pwb2 0.399 0.022 17.751 0.000 0.399 0.374
## .accept1 0.434 0.016 27.250 0.000 0.434 0.491
## .accept2 0.314 0.014 22.847 0.000 0.314 0.360
## .adapt1 0.288 0.012 23.299 0.000 0.288 0.480
## .adapt2 0.258 0.010 24.757 0.000 0.258 0.499
## .depl1 0.781 0.038 20.361 0.000 0.781 0.407
## .depl2 0.812 0.038 21.469 0.000 0.812 0.420
## .fs1 1.768 0.031 56.175 0.000 1.768 0.747
## .fs2 1.204 0.027 44.035 0.000 1.204 0.575
## .fs3 1.504 0.031 47.862 0.000 1.504 0.635
## .fs4 0.543 0.021 26.346 0.000 0.543 0.259
## .fs5 0.650 0.021 31.330 0.000 0.650 0.330
## .fs6 1.055 0.025 42.875 0.000 1.055 0.565
## leader 1.000 1.000 1.000
## growth 1.000 1.000 1.000
## wellbeing 1.000 1.000 1.000
## acceptance 1.000 1.000 1.000
## adaptability 1.000 1.000 1.000
## depletion 1.000 1.000 1.000
## foundational 1.000 1.000 1.000# Fit indices you’ll report (same set as Spring)
lavaan::fitMeasures(
fit_twb_7,
c("cfi","tli","rmsea","rmsea.ci.lower","rmsea.ci.upper","srmr","aic","bic")
)## cfi tli rmsea rmsea.ci.lower rmsea.ci.upper
## 0.954 0.949 0.048 0.046 0.049
## srmr aic bic
## 0.032 456586.779 457350.276CFI/TLI ~ .95 = excellent (especially for a 31-item multi-factor model)
RMSEA < .05 = excellent
SRMR ~ .03 = excellent
Lower-loading leadership items (still fine, but more expendable)
ls12 = .726 (lowest)
ls14 = .756
ls13 = .770
ls11 = .798
ls7 = .818
Higher residual variance = item is less explained by the factor.
Leadership residual variances (Std.all) include:
ls12 = .473 (highest residual; consistent with lowest loading)
ls14 = .429
ls13 = .407
ls11 = .364
ls7 = .331 These are the items most likely to be less central to the latent construct.
lead_tbl <- lavaan::standardizedSolution(fit_twb_7) %>%
dplyr::filter(op == "=~", lhs == "leader") %>%
dplyr::transmute(
item = rhs,
loading_std = est.std
) %>%
dplyr::arrange(dplyr::desc(loading_std))
lead_tbl## item loading_std
## 1 ls2 0.906
## 2 ls3 0.890
## 3 ls5 0.887
## 4 ls4 0.880
## 5 ls9 0.871
## 6 ls8 0.868
## 7 ls6 0.861
## 8 ls1 0.847
## 9 ls10 0.840
## 10 ls7 0.818
## 11 ls11 0.798
## 12 ls13 0.770
## 13 ls14 0.756
## 14 ls12 0.726resid_cor <- lavaan::resid(fit_twb_7, type = "cor")$cov
resid_lead <- resid_cor[items_twb_leader, items_twb_leader]
lead_pairs <- as.data.frame(as.table(resid_lead)) %>%
dplyr::rename(item1 = Var1, item2 = Var2, resid_cor = Freq) %>%
dplyr::filter(item1 < item2) %>%
dplyr::arrange(dplyr::desc(abs(resid_cor)))
lead_pairs %>% dplyr::slice(1:20)## [1] item1 item2 resid_cor
## <0 rows> (or 0-length row.names)Keeps the psychological core and preserves breadth.
Keep (8):
ls2 Takes my concerns seriously (strongest loading; core responsiveness)
ls3 Helps me problem-solve around student needs (instructional support)
ls4 Does what they say they will do (reliability/follow-through)
ls5 Integrates others’ input on meaningful issues (voice)
ls6 Responds to my requests directly (responsiveness/actionability; distinct from ls2)
ls8 Works to unify the staff (culture/cohesion)
ls9 Trust my admin to do the right thing for students (trust/values)
ls14 Staff share challenges without being judged (psychological safety)
Why not ls1/ls10/ls11/ls13/ls7/ls12 here?
ls1 overlaps conceptually with ls4 (consistency/follow-through)
ls10 overlaps with ls2/ls6 (helpfulness/responsiveness)
ls11 (boundaries) is valuable but can be “benefit/program dependent”; keep if it’s a priority
ls13 (info in time) is more operational than psychological (still useful, but different)
ls7 (fair discipline) is important but more student-discipline policy specific
ls12 (trust me to do my job) was the lowest loading; also narrower (autonomy)
items_twb_leader_short8 <- c("ls2","ls3","ls4","ls5","ls6","ls8","ls9","ls14")
items_twb_leader_short6 <- c("ls2","ls3","ls4","ls5","ls6","ls9")df <- df %>%
mutate(
twb_leader_score_short8 = rowMeans(select(., all_of(items_twb_leader_short8)), na.rm = TRUE),
twb_leader_score_short6 = rowMeans(select(., all_of(items_twb_leader_short6)), na.rm = TRUE)
)
cor(df$twb_leader_score, df$twb_leader_score_short8, use = "pairwise.complete.obs")## [1] 0.9882298cor(df$twb_leader_score, df$twb_leader_score_short6, use = "pairwise.complete.obs")## [1] 0.9748267mod_twb_7_short8 <- '
leader =~ ls2 + ls3 + ls4 + ls5 + ls6 + ls8 + ls9 + ls14
growth =~ growth1 + growth2 + growth3
wellbeing =~ pwb1 + pwb2
acceptance =~ accept1 + accept2
adaptability =~ adapt1 + adapt2
depletion =~ depl1 + depl2
foundational =~ fs1 + fs2 + fs3 + fs4 + fs5 + fs6
'
fit_twb_7_short8 <- lavaan::cfa(
model = mod_twb_7_short8, data = df,
estimator = "MLR", missing = "fiml", std.lv = TRUE
)
lavaan::fitMeasures(fit_twb_7_short8, c("cfi","tli","rmsea","srmr","aic","bic"))## cfi tli rmsea srmr aic bic
## 0.972 0.967 0.040 0.029 368899.907 369542.820- Direct comparison: Full vs. Short Leadership CFA 🔹 Full 14-item leadership model (from earlier)
CFI ≈ 0.954
TLI ≈ 0.949
RMSEA ≈ 0.048
SRMR ≈ 0.032
AIC / BIC ≈ 456,587 / 457,350
🔹 Short 8-item leadership model (new) CFI = 0.972 TLI = 0.967 RMSEA = 0.040 SRMR = 0.029 AIC = 368,899.9 BIC = 369,542.8
✅ What changed — and why it matters Metric What happened Interpretation CFI/TLI ↑ Increased substantially Less noise, clearer latent construct RMSEA ↓ Improved Better approximation of population covariance SRMR ↓ Improved Fewer residual misfits AIC/BIC ↓↓↓ Much lower Parsimony strongly favored
📌 Conclusion: The short-form leadership factor is statistically and substantively superior to the full version for this model.
- Why the short form works better (conceptually)
This is important for your theory section.
What you removed
Items that were:
More operational (e.g., information logistics)
More policy-specific (discipline fairness)
More administrative framing than psychological experience
What you kept
The short form isolates psychological leadership support, specifically:
Responsiveness
Trust
Voice
Follow-through
Psychological safety
Cultural cohesion
This aligns perfectly with your stated goal:
“a tighter psychological construct, not broader foundational conditions.”
You didn’t weaken the construct — you purified it.
To reduce respondent burden and sharpen construct validity, we evaluated a shortened version of the Responsive Leadership & Supportive Culture subscale. Guided by factor loadings, conceptual breadth, and psychological relevance, we retained eight items capturing responsiveness, trust, voice, follow-through, cohesion, and psychological safety. A confirmatory factor analysis indicated that the shortened subscale demonstrated improved model fit relative to the full 14-item version (CFI = .972, TLI = .967, RMSEA = .040, SRMR = .029), supporting its use as a parsimonious measure of psychologically supportive leadership.
cor(df$twb_leader_score, df$twb_leader_score_short8, use = "pairwise.complete.obs")## [1] 0.9882298A1. Overview of Measurement Development
The Teacher Well-Being (TWB) measure was designed as a multidimensional instrument assessing both psychological experiences of work and structural conditions supporting educators. The full instrument includes seven theoretically distinct subscales:
Responsive Leadership & Supportive Culture
Growth Orientation
Personal Well-Being
Acceptance
Adaptability
Depletion
Foundational Supports
All items were rated on a six-point Likert scale. Subscale scores were computed as mean scores of constituent items to preserve original scale metrics and avoid unequal weighting.
A2. Initial Reliability Assessment
Internal consistency was evaluated using Cronbach’s alpha for each subscale. Across Fall 2025 data, all subscales demonstrated acceptable to excellent reliability, consistent with prior administrations.
Leadership subscale (14 items): very high internal consistency
Growth Orientation (3 items): strong reliability
Personal Well-Being, Acceptance, Adaptability, Depletion (2 items each): adequate reliability given brevity
Foundational Supports (6 items): moderate reliability, consistent with its heterogeneous structural content
Item-level response distributions indicated no floor or ceiling effects that would threaten estimation.
A3. Full Measurement Model (Seven-Factor CFA)
A confirmatory factor analysis (CFA) specifying seven correlated latent factors was estimated using robust maximum likelihood (MLR) with full information maximum likelihood (FIML) to handle missing data.
Estimator: MLR Software: lavaan (v0.6-19) Sample size: N = 5,987 (Fall 2025)
Model fit (full leadership scale):
CFI = .954
TLI = .949
RMSEA = .048
SRMR = .032
All items loaded significantly on their intended factors (p < .001), supporting the theorized multidimensional structure.
A4. Rationale for Leadership Subscale Refinement
Despite strong psychometric performance, the Responsive Leadership & Supportive Culture subscale contained 14 items, making it disproportionately long relative to other TWB dimensions.
The refinement goal was to:
Reduce respondent burden
Improve parsimony
Retain a psychologically focused construct (responsiveness, trust, voice, safety)
Avoid conflating leadership experience with broader organizational logistics
Items were evaluated based on:
Standardized factor loadings
Conceptual redundancy
Coverage of core psychological facets
A5. Leadership Short-Form Construction
An 8-item short form was selected to preserve conceptual breadth while eliminating redundancy and operational specificity.
Retained facets:
Responsiveness to concerns and requests
Trust and ethical leadership
Voice and inclusion
Follow-through and reliability
Psychological safety
Cultural cohesion
The short form excludes items that primarily reflected administrative logistics (e.g., information flow timing) or policy-specific practices (e.g., disciplinary procedures), which were conceptually distinct from the targeted psychological leadership experience.
A6. Short-Form Measurement Model Results
The seven-factor CFA was re-estimated with the 8-item leadership short form, leaving all other subscales unchanged.
Model fit (short-form leadership):
CFI = .972
TLI = .967
RMSEA = .040
SRMR = .029
AIC = 368,899.91
BIC = 369,542.82
Interpretation:
Model fit improved across all indices relative to the full leadership model, indicating that the short form provides a cleaner and more parsimonious representation of the underlying leadership construct.
A7. Construct Validity and Parsimony
The improved fit indices and reduced information criteria suggest that the leadership short form:
Minimizes measurement noise
Reduces redundancy
Enhances construct coherence
Importantly, shortening the leadership scale did not alter the broader factor structure of the TWB measure, nor did it compromise interpretability of other subscales.
A8. Treatment of Foundational Supports
The Foundational Supports subscale was intentionally retained in its original 6-item form. Unlike psychological well-being constructs, foundational supports capture structural and policy-level conditions (e.g., compensation, staffing, professional development, equity of opportunity) that are inherently heterogeneous.
Maintaining full coverage was prioritized over parsimony to avoid construct underrepresentation and to preserve policy relevance.
A9. Summary and Recommendations
The Teacher Well-Being measure demonstrates strong psychometric properties across Fall 2025 data.
A shortened 8-item leadership subscale improves overall model fit while retaining conceptual breadth.
The refined measure balances psychological precision, respondent burden, and practical interpretability.
The short-form leadership subscale is recommended for future survey administrations and longitudinal analyses.
Supplementary Table S1 Leadership Subscale Refinement: Item Retention and Exclusion Item ID Item Content (Abbreviated) Retained in Short Form Primary Rationale ls1 Expectations of staff are consistent ❌ Removed Conceptual redundancy with follow-through (ls4) ls2 Takes my concerns seriously ✅ Retained Core psychological responsiveness; very high loading ls3 Helps problem-solve around student needs ✅ Retained Instructional and relational leadership ls4 Does what they say they will do ✅ Retained Reliability / trustworthiness ls5 Integrates others’ input on meaningful issues ✅ Retained Voice and inclusion ls6 Responds to my requests directly ✅ Retained Action-oriented responsiveness (distinct from ls2) ls7 Fair discipline processes ❌ Removed Policy-specific; less psychological in focus ls8 Works to unify the staff ✅ Retained Cultural cohesion ls9 Trust administrator to do the right thing ✅ Retained Moral trust / values ls10 Asks “How can I help?” ❌ Removed Redundant with responsiveness items (ls2, ls6) ls11 Encourages healthy boundaries ❌ Removed Valuable but narrower; context-dependent ls12 Trusts me to do my job ❌ Removed Lower loading; autonomy-specific ls13 Info provided in time to plan ❌ Removed Operational rather than psychological ls14 Staff share challenges without judgment ✅ Retained Psychological safety
Final short-form leadership subscale: 8 items (ls2, ls3, ls4, ls5, ls6, ls8, ls9, ls14)
Supplementary Table S2 Teacher Well-Being Subscales and Item Counts Subscale Items Construct Focus Responsive Leadership (Short Form) 8 Psychological leadership support Growth Orientation 3 Learning and development mindset Personal Well-Being 2 Life satisfaction outside work Acceptance 2 Emotional regulation Adaptability 2 Flexibility and composure Depletion 2 Cognitive and physical exhaustion Foundational Supports 6 Structural and policy conditions Supplementary Table S3 CFA Model Fit Comparison: Full vs. Short Leadership Model CFI TLI RMSEA SRMR AIC BIC Full Leadership (14 items) .954 .949 .048 .032 456,586.78 457,350.28 Short Leadership (8 items) .972 .967 .040 .029 368,899.91 369,542.82
Note. All models estimated using robust maximum likelihood (MLR) with FIML for missing data.
Supplementary Figure S1 (Description)
Seven-Factor CFA Model with Short Leadership Subscale
A path diagram depicting the seven correlated latent factors of Teacher Well-Being, including the refined 8-item Responsive Leadership factor. All standardized factor loadings exceed .70 for leadership items, and all inter-factor correlations are theoretically consistent (e.g., leadership positively associated with growth, well-being, and acceptance, and negatively associated with depletion).
(Figure can be generated using semPlot::semPaths() if desired.)
Results Section (1-Paragraph Summary)
A confirmatory factor analysis tested a seven-factor model of Teacher Well-Being. To reduce respondent burden and improve construct parsimony, the Responsive Leadership & Supportive Culture subscale was refined from 14 to 8 items based on factor loadings and conceptual relevance. The shortened model demonstrated excellent fit (CFI = .972, TLI = .967, RMSEA = .040, SRMR = .029) and improved upon the full leadership model across all fit indices. The retained leadership items captured core psychological dimensions of responsiveness, trust, voice, cohesion, and psychological safety, while excluding more operational or policy-specific indicators. All other subscales were retained in full. These results support the refined leadership subscale as a psychometrically strong and parsimonious measure within the broader Teacher Well-Being framework.
#create 14,8,6 item scale for the leadership factor
# Leadership item sets
items_leader_14 <- paste0("ls", 1:14)
items_leader_8 <- c("ls2","ls3","ls4","ls5","ls6","ls8","ls9","ls14")
items_leader_6 <- c("ls2","ls3","ls4","ls5","ls6","ls9")
# Compute mean scores (not sums) for each version
df <- df %>%
mutate(
twb_leader_14 = rowMeans(select(., all_of(items_leader_14)), na.rm = TRUE),
twb_leader_8 = rowMeans(select(., all_of(items_leader_8)), na.rm = TRUE),
twb_leader_6 = rowMeans(select(., all_of(items_leader_6)), na.rm = TRUE)
)
# Correlations with full 14-item score (key validity check)
cor(df$twb_leader_14, df$twb_leader_8, use = "pairwise.complete.obs")## [1] 0.9882298cor(df$twb_leader_14, df$twb_leader_6, use = "pairwise.complete.obs")## [1] 0.9748267# Leadership item sets (explicit)
items_leader_14 <- c("ls1","ls2","ls3","ls4","ls5","ls6","ls7","ls8","ls9","ls10","ls11","ls12","ls13","ls14")
items_leader_8 <- c("ls2","ls3","ls4","ls5","ls6","ls8","ls9","ls14")
items_leader_6 <- c("ls2","ls3","ls4","ls5","ls6","ls9")
# -----------------------------
# Model 1: Leader = 14 items
# -----------------------------
mod_twb_7_leader14 <- '
leader =~ ls1 + ls2 + ls3 + ls4 + ls5 + ls6 + ls7 + ls8 + ls9 + ls10 + ls11 + ls12 + ls13 + ls14
growth =~ growth1 + growth2 + growth3
wellbeing =~ pwb1 + pwb2
acceptance =~ accept1 + accept2
adaptability =~ adapt1 + adapt2
depletion =~ depl1 + depl2
foundational =~ fs1 + fs2 + fs3 + fs4 + fs5 + fs6
'
fit_leader14 <- lavaan::cfa(
model = mod_twb_7_leader14,
data = df,
estimator = "MLR",
missing = "fiml",
std.lv = TRUE
)
# -----------------------------
# Model 2: Leader = 8 items
# -----------------------------
mod_twb_7_leader8 <- '
leader =~ ls2 + ls3 + ls4 + ls5 + ls6 + ls8 + ls9 + ls14
growth =~ growth1 + growth2 + growth3
wellbeing =~ pwb1 + pwb2
acceptance =~ accept1 + accept2
adaptability =~ adapt1 + adapt2
depletion =~ depl1 + depl2
foundational =~ fs1 + fs2 + fs3 + fs4 + fs5 + fs6
'
fit_leader8 <- lavaan::cfa(
model = mod_twb_7_leader8,
data = df,
estimator = "MLR",
missing = "fiml",
std.lv = TRUE
)
# -----------------------------
# Model 3: Leader = 6 items
# -----------------------------
mod_twb_7_leader6 <- '
leader =~ ls2 + ls3 + ls4 + ls5 + ls6 + ls9
growth =~ growth1 + growth2 + growth3
wellbeing =~ pwb1 + pwb2
acceptance =~ accept1 + accept2
adaptability =~ adapt1 + adapt2
depletion =~ depl1 + depl2
foundational =~ fs1 + fs2 + fs3 + fs4 + fs5 + fs6
'
fit_leader6 <- lavaan::cfa(
model = mod_twb_7_leader6,
data = df,
estimator = "MLR",
missing = "fiml",
std.lv = TRUE
)
# Quick check: models exist
fit_leader14## lavaan 0.6-19 ended normally after 113 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 114
##
## Used Total
## Number of observations 5987 6062
## Number of missing patterns 222
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 6004.780 4509.157
## Degrees of freedom 413 413
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.332
## Yuan-Bentler correction (Mplus variant)fit_leader8## lavaan 0.6-19 ended normally after 97 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 96
##
## Used Total
## Number of observations 5985 6062
## Number of missing patterns 173
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 2670.424 2089.008
## Degrees of freedom 254 254
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.278
## Yuan-Bentler correction (Mplus variant)fit_leader6## lavaan 0.6-19 ended normally after 93 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 90
##
## Used Total
## Number of observations 5985 6062
## Number of missing patterns 145
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 1616.891 1281.297
## Degrees of freedom 209 209
## P-value (Chi-square) 0.000 0.000
## Scaling correction factor 1.262
## Yuan-Bentler correction (Mplus variant)# Fit measures to extract (same for all models)
fit_set <- c("cfi","tli","rmsea","rmsea.ci.lower","rmsea.ci.upper","srmr","aic","bic")
fm14 <- lavaan::fitMeasures(fit_leader14, fit_set)
fm8 <- lavaan::fitMeasures(fit_leader8, fit_set)
fm6 <- lavaan::fitMeasures(fit_leader6, fit_set)
# Convert to base matrix (guaranteed stable)
fit_mat <- rbind(
Leader_14 = fm14,
Leader_8 = fm8,
Leader_6 = fm6
)
fit_df <- as.data.frame(round(fit_mat, 3))
# Print the table
fit_df## cfi tli rmsea rmsea.ci.lower rmsea.ci.upper srmr aic
## Leader_14 0.954 0.949 0.048 0.046 0.049 0.032 456586.8
## Leader_8 0.972 0.967 0.040 0.039 0.041 0.029 368899.9
## Leader_6 0.981 0.977 0.034 0.032 0.035 0.027 339279.7
## bic
## Leader_14 457350.3
## Leader_8 369542.8
## Leader_6 339882.5As expected with a large sample, chi-square tests were statistically significant for all models. Accordingly, model evaluation focused on comparative and approximate fit indices. Reducing the Responsive Leadership & Supportive Culture subscale from 14 to 8 items resulted in improved fit across all indices (CFI, TLI, RMSEA, SRMR) and substantially lower information criteria, indicating a more parsimonious and coherent representation of psychologically supportive leadership. A further reduction to 6 items yielded marginal additional gains in fit but at the cost of reduced conceptual breadth. The 8-item leadership subscale was therefore selected as the optimal balance of psychometric performance and construct coverage. Does shortening the leadership factor weaken, preserve, or improve the ability of overall Teacher Well-Being to predict key outcomes?
# Full TWB leadership (14)
items_leader_14 <- paste0("ls", 1:14)
# Short TWB leadership (8)
items_leader_8 <- c("ls2","ls3","ls4","ls5","ls6","ls8","ls9","ls14")
# Non-leadership TWB items (identical across versions)
items_twb_nonleader <- c(
items_twb_growth,
items_twb_wellbeing,
items_twb_acceptance,
items_twb_adaptability,
items_twb_depletion,
items_twb_foundational
)df <- df %>%
mutate(
# TWB composite with 14-item leadership
twb_total_14 = rowMeans(
select(., all_of(c(items_leader_14, items_twb_nonleader))),
na.rm = TRUE
),
# TWB composite with 8-item leadership
twb_total_8 = rowMeans(
select(., all_of(c(items_leader_8, items_twb_nonleader))),
na.rm = TRUE
)
)
# Sanity check: they should be extremely highly correlated
cor(df$twb_total_14, df$twb_total_8, use = "pairwise.complete.obs")## [1] 0.9861736m_js_14 <- lm(job_sat ~ twb_total_14, data = df)
m_js_8 <- lm(job_sat ~ twb_total_8, data = df)
summary(m_js_14)##
## Call:
## lm(formula = job_sat ~ twb_total_14, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.5612 -0.8189 0.1700 0.9981 8.3214
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.31083 0.12749 -2.438 0.0148 *
## twb_total_14 1.66686 0.02847 58.550 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.558 on 5475 degrees of freedom
## (585 observations deleted due to missingness)
## Multiple R-squared: 0.385, Adjusted R-squared: 0.3849
## F-statistic: 3428 on 1 and 5475 DF, p-value: < 2.2e-16summary(m_js_8)##
## Call:
## lm(formula = job_sat ~ twb_total_8, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.6455 -0.8079 0.1719 0.9700 8.9692
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.32973 0.14064 -9.455 <2e-16 ***
## twb_total_8 1.90368 0.03159 60.254 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.54 on 5475 degrees of freedom
## (585 observations deleted due to missingness)
## Multiple R-squared: 0.3987, Adjusted R-squared: 0.3986
## F-statistic: 3631 on 1 and 5475 DF, p-value: < 2.2e-16broom::glance(m_js_14)[c("r.squared","adj.r.squared")]## # A tibble: 1 × 2
## r.squared adj.r.squared
## <dbl> <dbl>
## 1 0.385 0.385broom::glance(m_js_8)[c("r.squared","adj.r.squared")]## # A tibble: 1 × 2
## r.squared adj.r.squared
## <dbl> <dbl>
## 1 0.399 0.399m_stress_14 <- lm(job_stress ~ twb_total_14, data = df)
m_stress_8 <- lm(job_stress ~ twb_total_8, data = df)
summary(m_stress_14)##
## Call:
## lm(formula = job_stress ~ twb_total_14, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.587 -1.214 0.317 1.538 3.959
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.1586 0.1715 59.24 <2e-16 ***
## twb_total_14 -0.6863 0.0383 -17.92 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.099 on 5478 degrees of freedom
## (582 observations deleted due to missingness)
## Multiple R-squared: 0.05537, Adjusted R-squared: 0.0552
## F-statistic: 321.1 on 1 and 5478 DF, p-value: < 2.2e-16summary(m_stress_8)##
## Call:
## lm(formula = job_stress ~ twb_total_8, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.6697 -1.2238 0.3105 1.5518 4.1278
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.58701 0.19113 55.39 <2e-16 ***
## twb_total_8 -0.78580 0.04294 -18.30 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.097 on 5478 degrees of freedom
## (582 observations deleted due to missingness)
## Multiple R-squared: 0.05761, Adjusted R-squared: 0.05744
## F-statistic: 334.9 on 1 and 5478 DF, p-value: < 2.2e-16broom::glance(m_stress_14)[c("r.squared","adj.r.squared")]## # A tibble: 1 × 2
## r.squared adj.r.squared
## <dbl> <dbl>
## 1 0.0554 0.0552broom::glance(m_stress_8)[c("r.squared","adj.r.squared")]## # A tibble: 1 × 2
## r.squared adj.r.squared
## <dbl> <dbl>
## 1 0.0576 0.0574m_ret_14 <- lm(retention_score ~ twb_total_14, data = df)
m_ret_8 <- lm(retention_score ~ twb_total_8, data = df)
summary(m_ret_14)##
## Call:
## lm(formula = retention_score ~ twb_total_14, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4109 -0.4594 0.1856 0.6073 2.6384
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.37353 0.07104 19.34 <2e-16 ***
## twb_total_14 0.56383 0.01588 35.51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8651 on 5433 degrees of freedom
## (627 observations deleted due to missingness)
## Multiple R-squared: 0.1884, Adjusted R-squared: 0.1882
## F-statistic: 1261 on 1 and 5433 DF, p-value: < 2.2e-16summary(m_ret_8)##
## Call:
## lm(formula = retention_score ~ twb_total_8, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4678 -0.4678 0.1814 0.5940 2.9101
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.96684 0.07858 12.30 <2e-16 ***
## twb_total_8 0.65807 0.01767 37.25 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8571 on 5433 degrees of freedom
## (627 observations deleted due to missingness)
## Multiple R-squared: 0.2034, Adjusted R-squared: 0.2033
## F-statistic: 1387 on 1 and 5433 DF, p-value: < 2.2e-16broom::glance(m_ret_14)[c("r.squared","adj.r.squared")]## # A tibble: 1 × 2
## r.squared adj.r.squared
## <dbl> <dbl>
## 1 0.188 0.188broom::glance(m_ret_8)[c("r.squared","adj.r.squared")]## # A tibble: 1 × 2
## r.squared adj.r.squared
## <dbl> <dbl>
## 1 0.203 0.203m_nps_14 <- lm(nps_score ~ twb_total_14, data = df)
m_nps_8 <- lm(nps_score ~ twb_total_8, data = df)
summary(m_nps_14)##
## Call:
## lm(formula = nps_score ~ twb_total_14, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.6502 -1.1766 0.2956 1.3420 10.6624
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.35435 0.16725 -20.06 <2e-16 ***
## twb_total_14 2.25541 0.03736 60.37 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.047 on 5467 degrees of freedom
## (593 observations deleted due to missingness)
## Multiple R-squared: 0.4, Adjusted R-squared: 0.3999
## F-statistic: 3645 on 1 and 5467 DF, p-value: < 2.2e-16summary(m_nps_8)##
## Call:
## lm(formula = nps_score ~ twb_total_8, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.7440 -1.1514 0.2766 1.3383 11.4617
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.62544 0.18569 -24.91 <2e-16 ***
## twb_total_8 2.55142 0.04172 61.15 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.036 on 5467 degrees of freedom
## (593 observations deleted due to missingness)
## Multiple R-squared: 0.4062, Adjusted R-squared: 0.4061
## F-statistic: 3739 on 1 and 5467 DF, p-value: < 2.2e-16broom::glance(m_nps_14)[c("r.squared","adj.r.squared")]## # A tibble: 1 × 2
## r.squared adj.r.squared
## <dbl> <dbl>
## 1 0.400 0.400broom::glance(m_nps_8)[c("r.squared","adj.r.squared")]## # A tibble: 1 × 2
## r.squared adj.r.squared
## <dbl> <dbl>
## 1 0.406 0.406pred_tbl <- tibble(
outcome = c("Job Satisfaction","Job Stress","Retention","NPS"),
r2_14 = c(
summary(m_js_14)$r.squared,
summary(m_stress_14)$r.squared,
summary(m_ret_14)$r.squared,
summary(m_nps_14)$r.squared
),
r2_8 = c(
summary(m_js_8)$r.squared,
summary(m_stress_8)$r.squared,
summary(m_ret_8)$r.squared,
summary(m_nps_8)$r.squared
),
beta_14 = c(
coef(m_js_14)[2],
coef(m_stress_14)[2],
coef(m_ret_14)[2],
coef(m_nps_14)[2]
),
beta_8 = c(
coef(m_js_8)[2],
coef(m_stress_8)[2],
coef(m_ret_8)[2],
coef(m_nps_8)[2]
)
) %>%
mutate(
delta_r2 = r2_8 - r2_14
)
print(pred_tbl)## # A tibble: 4 × 6
## outcome r2_14 r2_8 beta_14 beta_8 delta_r2
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Job Satisfaction 0.385 0.399 1.67 1.90 0.0137
## 2 Job Stress 0.0554 0.0576 -0.686 -0.786 0.00224
## 3 Retention 0.188 0.203 0.564 0.658 0.0150
## 4 NPS 0.400 0.406 2.26 2.55 0.00616knitr::kable(
pred_tbl,
digits = 3,
caption = "Predictive Validity of TWB Composite Scores: 14-item vs 8-item Leadership"
)| outcome | r2_14 | r2_8 | beta_14 | beta_8 | delta_r2 |
|---|---|---|---|---|---|
| Job Satisfaction | 0.385 | 0.399 | 1.667 | 1.904 | 0.014 |
| Job Stress | 0.055 | 0.058 | -0.686 | -0.786 | 0.002 |
| Retention | 0.188 | 0.203 | 0.564 | 0.658 | 0.015 |
| NPS | 0.400 | 0.406 | 2.255 | 2.551 | 0.006 |
To evaluate predictive validity, we compared two versions of the Teacher Well-Being (TWB) composite score: one incorporating the original 14-item Responsive Leadership subscale and one incorporating a refined 8-item version. Across all four criterion outcomes—job satisfaction, job stress, retention intentions, and net promoter scores—the shortened leadership version demonstrated equal or superior predictive validity. The 8-item composite explained more variance in job satisfaction (R² = .399 vs. .385), retention intentions (R² = .203 vs. .188), and net promoter scores (R² = .406 vs. .400), with a small improvement for job stress (R² = .058 vs. .055). Regression coefficients were consistently larger in magnitude for the refined scale. These findings indicate that shortening the leadership subscale did not compromise predictive validity and instead yielded a more parsimonious and psychologically precise measure of teacher well-being.
Use the 8-item leadership scale as the default going forward.
Main text: 8-item version
Supplement: 14-item version documented
Policy briefs & dashboards: 8-item version
Longitudinal work: 8-item version (less fatigue, more stability)
# -----------------------------------
# Outcome variables
# -----------------------------------
outcomes <- c(
"job_sat", # Job satisfaction
"job_stress", # Job stress
"retention_score",
"nps_score"
)
df %>%
summarise(across(all_of(outcomes), ~ mean(.x, na.rm = TRUE)))## # A tibble: 1 × 4
## job_sat job_stress retention_score nps_score
## <dbl> <dbl> <dbl> <dbl>
## 1 7.05 7.13 3.86 6.60# -----------------------------------
# TWB predictors
# -----------------------------------
twb_predictors <- c(
"twb_leader_score_short8",
"growth_score",
"wellbeing_score",
"acceptance_score",
"adaptability_score",
"depletion_score",
"foundational_support_score",
"twb_composite_score"
)
# -----------------------------------
# LMC predictors
# -----------------------------------
lmc_predictors <- c(
"union_view_score",
"school_lmc_score",
"district_lmc_score",
"lmc_composite_score"
)# -----------------------------------
# Correlation matrix: TWB + LMC + outcomes
# -----------------------------------
corr_vars <- c(
twb_predictors,
lmc_predictors,
outcomes
)
cor_mat <- cor(
df[, corr_vars],
use = "pairwise.complete.obs"
)
round(cor_mat, 3)## twb_leader_score_short8 growth_score wellbeing_score
## twb_leader_score_short8 1.000 0.328 0.118
## growth_score 0.328 1.000 0.343
## wellbeing_score 0.118 0.343 1.000
## acceptance_score 0.280 0.370 0.320
## adaptability_score 0.262 0.428 0.317
## depletion_score -0.120 -0.171 -0.394
## foundational_support_score 0.536 0.316 0.201
## twb_composite_score 0.927 0.493 0.270
## union_view_score 0.194 0.156 0.120
## school_lmc_score 0.533 0.238 0.123
## district_lmc_score 0.387 0.251 0.140
## lmc_composite_score 0.440 0.253 0.151
## job_sat 0.524 0.349 0.236
## job_stress -0.170 -0.109 -0.208
## retention_score 0.344 0.275 0.132
## nps_score 0.558 0.329 0.135
## acceptance_score adaptability_score depletion_score
## twb_leader_score_short8 0.280 0.262 -0.120
## growth_score 0.370 0.428 -0.171
## wellbeing_score 0.320 0.317 -0.394
## acceptance_score 1.000 0.649 -0.296
## adaptability_score 0.649 1.000 -0.241
## depletion_score -0.296 -0.241 1.000
## foundational_support_score 0.310 0.271 -0.289
## twb_composite_score 0.444 0.422 -0.155
## union_view_score 0.089 0.108 -0.041
## school_lmc_score 0.204 0.185 -0.092
## district_lmc_score 0.189 0.189 -0.113
## lmc_composite_score 0.188 0.188 -0.099
## job_sat 0.313 0.281 -0.333
## job_stress -0.247 -0.195 0.485
## retention_score 0.202 0.163 -0.187
## nps_score 0.241 0.236 -0.229
## foundational_support_score twb_composite_score
## twb_leader_score_short8 0.536 0.927
## growth_score 0.316 0.493
## wellbeing_score 0.201 0.270
## acceptance_score 0.310 0.444
## adaptability_score 0.271 0.422
## depletion_score -0.289 -0.155
## foundational_support_score 1.000 0.736
## twb_composite_score 0.736 1.000
## union_view_score 0.218 0.237
## school_lmc_score 0.422 0.557
## district_lmc_score 0.490 0.472
## lmc_composite_score 0.453 0.500
## job_sat 0.618 0.621
## job_stress -0.337 -0.235
## retention_score 0.482 0.434
## nps_score 0.608 0.632
## union_view_score school_lmc_score district_lmc_score
## twb_leader_score_short8 0.194 0.533 0.387
## growth_score 0.156 0.238 0.251
## wellbeing_score 0.120 0.123 0.140
## acceptance_score 0.089 0.204 0.189
## adaptability_score 0.108 0.185 0.189
## depletion_score -0.041 -0.092 -0.113
## foundational_support_score 0.218 0.422 0.490
## twb_composite_score 0.237 0.557 0.472
## union_view_score 1.000 0.529 0.579
## school_lmc_score 0.529 1.000 0.645
## district_lmc_score 0.579 0.645 1.000
## lmc_composite_score 0.786 0.850 0.910
## job_sat 0.217 0.397 0.392
## job_stress -0.023 -0.101 -0.135
## retention_score 0.205 0.284 0.324
## nps_score 0.261 0.420 0.440
## lmc_composite_score job_sat job_stress
## twb_leader_score_short8 0.440 0.524 -0.170
## growth_score 0.253 0.349 -0.109
## wellbeing_score 0.151 0.236 -0.208
## acceptance_score 0.188 0.313 -0.247
## adaptability_score 0.188 0.281 -0.195
## depletion_score -0.099 -0.333 0.485
## foundational_support_score 0.453 0.618 -0.337
## twb_composite_score 0.500 0.621 -0.235
## union_view_score 0.786 0.217 -0.023
## school_lmc_score 0.850 0.397 -0.101
## district_lmc_score 0.910 0.392 -0.135
## lmc_composite_score 1.000 0.399 -0.106
## job_sat 0.399 1.000 -0.394
## job_stress -0.106 -0.394 1.000
## retention_score 0.325 0.525 -0.220
## nps_score 0.443 0.722 -0.291
## retention_score nps_score
## twb_leader_score_short8 0.344 0.558
## growth_score 0.275 0.329
## wellbeing_score 0.132 0.135
## acceptance_score 0.202 0.241
## adaptability_score 0.163 0.236
## depletion_score -0.187 -0.229
## foundational_support_score 0.482 0.608
## twb_composite_score 0.434 0.632
## union_view_score 0.205 0.261
## school_lmc_score 0.284 0.420
## district_lmc_score 0.324 0.440
## lmc_composite_score 0.325 0.443
## job_sat 0.525 0.722
## job_stress -0.220 -0.291
## retention_score 1.000 0.536
## nps_score 0.536 1.000Zero-order correlations supported the theorized nomological network. The TWB composite was strongly related to responsive leadership (r = .927) and foundational supports (r = .736), indicating that educator well-being in this sample was anchored in both psychological leadership climate and structural working conditions. LMC showed moderate associations with TWB (r = .500) and particularly with leadership (r = .440) and foundational supports (r = .453), consistent with LMC functioning as a structural-contextual indicator of system-level collaboration. In relation to outcomes, TWB and foundational supports were strongly correlated with job satisfaction (rs = .621 and .618) and net promoter scores (rs = .632 and .608), while retention intentions were most closely associated with foundational supports (r = .482). Job stress demonstrated a distinct pattern, correlating most strongly with depletion (r = .485) and more modestly with leadership and LMC. Together, these patterns suggest that leadership and foundational conditions are central correlates of global job attitudes and advocacy outcomes, whereas stress is more tightly linked to exhaustion and personal bandwidth.
m_js_sub <- lm(
job_sat ~
twb_leader_score_short8 +
growth_score +
wellbeing_score +
acceptance_score +
adaptability_score +
depletion_score +
foundational_support_score +
lmc_composite_score,
data = df
)
summary(m_js_sub)##
## Call:
## lm(formula = job_sat ~ twb_leader_score_short8 + growth_score +
## wellbeing_score + acceptance_score + adaptability_score +
## depletion_score + foundational_support_score + lmc_composite_score,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.3100 -0.7722 0.1210 0.8961 8.2091
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1868010 0.2298249 5.164 2.51e-07 ***
## twb_leader_score_short8 0.3795162 0.0200224 18.955 < 2e-16 ***
## growth_score 0.2477704 0.0315322 7.858 4.69e-15 ***
## wellbeing_score 0.0353135 0.0240144 1.471 0.14148
## acceptance_score 0.0845680 0.0319543 2.647 0.00816 **
## adaptability_score -0.0008296 0.0413625 -0.020 0.98400
## depletion_score -0.2533059 0.0179619 -14.102 < 2e-16 ***
## foundational_support_score 0.6801340 0.0232615 29.239 < 2e-16 ***
## lmc_composite_score 0.1908785 0.0274834 6.945 4.23e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.429 on 5382 degrees of freedom
## (671 observations deleted due to missingness)
## Multiple R-squared: 0.4818, Adjusted R-squared: 0.4811
## F-statistic: 625.6 on 8 and 5382 DF, p-value: < 2.2e-16In a multivariate regression predicting job satisfaction, the full model explained 48% of the variance (adjusted R² = .48). Foundational supports emerged as the strongest predictor (β = .68, p < .001), followed by responsive leadership (β = .38, p < .001). Depletion was a substantial negative predictor (β = −.25, p < .001), while growth orientation (β = .25, p < .001) and acceptance (β = .09, p = .008) showed smaller positive associations. Importantly, labor–management collaboration remained a significant independent predictor of job satisfaction (β = .19, p < .001), even after controlling for leadership, foundational supports, and individual well-being factors. Personal well-being outside of work and adaptability were not significant in the full model, suggesting that job satisfaction is primarily driven by organizational conditions, leadership climate, and energy depletion rather than general dispositional traits.
m_js_comp <- lm(
job_sat ~ twb_composite_score + lmc_composite_score,
data = df
)
summary(m_js_comp)##
## Call:
## lm(formula = job_sat ~ twb_composite_score + lmc_composite_score,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.7688 -0.8334 0.1546 0.9844 8.8008
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.88277 0.14005 -6.303 3.15e-10 ***
## twb_composite_score 1.51114 0.03282 46.047 < 2e-16 ***
## lmc_composite_score 0.27838 0.02919 9.537 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.545 on 5433 degrees of freedom
## (626 observations deleted due to missingness)
## Multiple R-squared: 0.3952, Adjusted R-squared: 0.395
## F-statistic: 1775 on 2 and 5433 DF, p-value: < 2.2e-16In a parsimonious regression model, overall Teacher Well-Being and Labor-Management Collaboration jointly explained 39.5% of the variance in job satisfaction. Teacher Well-Being was the strongest predictor (β = 1.51, p < .001), indicating that educators with higher overall well-being reported substantially greater job satisfaction. Labor-Management Collaboration also showed a significant independent association with job satisfaction (β = 0.28, p < .001), even after accounting for overall well-being. These findings suggest that collaborative labor–management structures contribute uniquely to job satisfaction above and beyond individual and climate-based well-being.
m_stress_sub <- lm(
job_stress ~
twb_leader_score_short8 +
growth_score +
wellbeing_score +
acceptance_score +
adaptability_score +
depletion_score +
foundational_support_score +
lmc_composite_score,
data = df
)
summary(m_stress_sub)##
## Call:
## lm(formula = job_stress ~ twb_leader_score_short8 + growth_score +
## wellbeing_score + acceptance_score + adaptability_score +
## depletion_score + foundational_support_score + lmc_composite_score,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.7169 -1.0678 0.2155 1.2459 5.9745
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.956543 0.292918 20.335 < 2e-16 ***
## twb_leader_score_short8 -0.029154 0.025536 -1.142 0.25363
## growth_score 0.170729 0.040221 4.245 2.22e-05 ***
## wellbeing_score 0.005715 0.030628 0.187 0.85198
## acceptance_score -0.178675 0.040748 -4.385 1.18e-05 ***
## adaptability_score -0.074817 0.052736 -1.419 0.15604
## depletion_score 0.715365 0.022909 31.226 < 2e-16 ***
## foundational_support_score -0.444782 0.029669 -14.991 < 2e-16 ***
## lmc_composite_score 0.109912 0.035040 3.137 0.00172 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.823 on 5384 degrees of freedom
## (669 observations deleted due to missingness)
## Multiple R-squared: 0.2862, Adjusted R-squared: 0.2851
## F-statistic: 269.8 on 8 and 5384 DF, p-value: < 2.2e-16m_stress_comp <- lm(
job_stress ~ twb_composite_score + lmc_composite_score,
data = df
)
summary(m_stress_comp)##
## Call:
## lm(formula = job_stress ~ twb_composite_score + lmc_composite_score,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.522 -1.207 0.310 1.534 3.930
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.06385 0.19006 52.952 <2e-16 ***
## twb_composite_score -0.70639 0.04451 -15.870 <2e-16 ***
## lmc_composite_score 0.04075 0.03964 1.028 0.304
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.1 on 5435 degrees of freedom
## (624 observations deleted due to missingness)
## Multiple R-squared: 0.05502, Adjusted R-squared: 0.05467
## F-statistic: 158.2 on 2 and 5435 DF, p-value: < 2.2e-16In a multivariate regression predicting job stress, the model explained 28.5% of the variance. Depletion emerged as the dominant predictor (β = .72, p < .001), indicating that exhaustion and cognitive overload were the primary drivers of perceived stress. Foundational supports were also strongly associated with lower stress (β = −.44, p < .001), suggesting that structural conditions buffer stress beyond individual experience. Acceptance showed a modest protective association (β = −.18, p < .001), while growth orientation was positively associated with stress (β = .17, p < .001), consistent with strain associated with high engagement. In contrast, responsive leadership and labor–management collaboration were not significant predictors once depletion and foundational conditions were included, indicating that their associations with stress operate indirectly through these more proximal mechanisms.
In a parsimonious model including only the TWB composite and labor–management collaboration, the explained variance in job stress was substantially lower (adjusted R² = .055). While overall well-being was negatively associated with stress, labor–management collaboration was not a significant predictor. This contrast highlights the importance of modeling specific well-being mechanisms—particularly depletion—when examining stress outcomes.
m_ret_sub <- lm(
retention_score ~
twb_leader_score_short8 +
growth_score +
wellbeing_score +
acceptance_score +
adaptability_score +
depletion_score +
foundational_support_score +
lmc_composite_score,
data = df
)
summary(m_ret_sub)##
## Call:
## lm(formula = retention_score ~ twb_leader_score_short8 + growth_score +
## wellbeing_score + acceptance_score + adaptability_score +
## depletion_score + foundational_support_score + lmc_composite_score,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6306 -0.4318 0.1489 0.5621 2.3920
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.64608 0.13276 12.399 < 2e-16 ***
## twb_leader_score_short8 0.05938 0.01156 5.136 2.90e-07 ***
## growth_score 0.15434 0.01819 8.484 < 2e-16 ***
## wellbeing_score -0.01941 0.01385 -1.401 0.16117
## acceptance_score 0.03523 0.01839 1.915 0.05549 .
## adaptability_score -0.06951 0.02378 -2.923 0.00348 **
## depletion_score -0.04263 0.01037 -4.112 3.98e-05 ***
## foundational_support_score 0.31182 0.01344 23.201 < 2e-16 ***
## lmc_composite_score 0.11378 0.01596 7.127 1.16e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.821 on 5329 degrees of freedom
## (724 observations deleted due to missingness)
## Multiple R-squared: 0.2658, Adjusted R-squared: 0.2647
## F-statistic: 241.2 on 8 and 5329 DF, p-value: < 2.2e-16m_ret_comp <- lm(
retention_score ~ twb_composite_score + lmc_composite_score,
data = df
)
summary(m_ret_comp)##
## Call:
## lm(formula = retention_score ~ twb_composite_score + lmc_composite_score,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5627 -0.4519 0.1858 0.5992 2.6741
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.04333 0.07799 13.38 <2e-16 ***
## twb_composite_score 0.47155 0.01834 25.72 <2e-16 ***
## lmc_composite_score 0.16318 0.01635 9.98 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8562 on 5376 degrees of freedom
## (683 observations deleted due to missingness)
## Multiple R-squared: 0.2035, Adjusted R-squared: 0.2032
## F-statistic: 686.7 on 2 and 5376 DF, p-value: < 2.2e-16In a multivariate regression predicting retention intentions, the model explained 26.5% of the variance. Foundational supports emerged as the strongest predictor (β = .31, p < .001), indicating that structural conditions such as compensation, staffing, and professional support are central to educators’ decisions to remain in their roles. Labor–management collaboration was also a significant positive predictor (β = .11, p < .001), demonstrating an independent contribution of system-level collaboration beyond school leadership and individual well-being. Responsive leadership showed a smaller but significant association (β = .06, p < .001). Growth orientation was positively associated with retention (β = .15, p < .001), whereas depletion was negatively associated (β = −.04, p < .001). Together, these findings suggest that retention reflects a combination of structural feasibility, institutional trust, and motivational commitment.
In a more parsimonious model including only the TWB composite and labor–management collaboration, explained variance in retention was lower (adjusted R² = .203). Although both predictors were significant, this model obscured the distinct psychological and structural mechanisms identified in the decomposed model, underscoring the value of modeling specific well-being dimensions when examining retention outcomes. Stress = energy depletion
Satisfaction = relational + structural appraisal
Retention = feasibility + commitment + system trust
m_nps_sub <- lm(
nps_score ~
twb_leader_score_short8 +
growth_score +
wellbeing_score +
acceptance_score +
adaptability_score +
depletion_score +
foundational_support_score +
lmc_composite_score,
data = df
)
summary(m_nps_sub)##
## Call:
## lm(formula = nps_score ~ twb_leader_score_short8 + growth_score +
## wellbeing_score + acceptance_score + adaptability_score +
## depletion_score + foundational_support_score + lmc_composite_score,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.5539 -1.0942 0.2119 1.2465 10.0452
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.36075 0.30758 -4.424 9.88e-06 ***
## twb_leader_score_short8 0.59387 0.02681 22.154 < 2e-16 ***
## growth_score 0.35713 0.04224 8.454 < 2e-16 ***
## wellbeing_score -0.15659 0.03215 -4.870 1.15e-06 ***
## acceptance_score -0.06972 0.04272 -1.632 0.103
## adaptability_score 0.05020 0.05530 0.908 0.364
## depletion_score -0.19226 0.02404 -7.996 1.56e-15 ***
## foundational_support_score 0.87543 0.03113 28.122 < 2e-16 ***
## lmc_composite_score 0.43429 0.03680 11.801 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.915 on 5385 degrees of freedom
## (668 observations deleted due to missingness)
## Multiple R-squared: 0.4746, Adjusted R-squared: 0.4738
## F-statistic: 608.1 on 8 and 5385 DF, p-value: < 2.2e-16m_nps_comp <- lm(
nps_score ~ twb_composite_score + lmc_composite_score,
data = df
)
summary(m_nps_comp)##
## Call:
## lm(formula = nps_score ~ twb_composite_score + lmc_composite_score,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.0534 -1.1708 0.2511 1.3196 11.5964
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.47202 0.18194 -24.58 <2e-16 ***
## twb_composite_score 1.96337 0.04259 46.10 <2e-16 ***
## lmc_composite_score 0.53225 0.03795 14.02 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.01 on 5436 degrees of freedom
## (623 observations deleted due to missingness)
## Multiple R-squared: 0.4222, Adjusted R-squared: 0.422
## F-statistic: 1986 on 2 and 5436 DF, p-value: < 2.2e-16model_compare <- tibble(
outcome = c("Job Satisfaction","Job Stress","Retention","NPS"),
r2_sub = c(
summary(m_js_sub)$r.squared,
summary(m_stress_sub)$r.squared,
summary(m_ret_sub)$r.squared,
summary(m_nps_sub)$r.squared
),
r2_comp = c(
summary(m_js_comp)$r.squared,
summary(m_stress_comp)$r.squared,
summary(m_ret_comp)$r.squared,
summary(m_nps_comp)$r.squared
)
) %>%
mutate(delta_r2 = r2_sub - r2_comp)
model_compare %>%
mutate(across(where(is.numeric), ~ round(.x, 3)))## # A tibble: 4 × 4
## outcome r2_sub r2_comp delta_r2
## <chr> <dbl> <dbl> <dbl>
## 1 Job Satisfaction 0.482 0.395 0.087
## 2 Job Stress 0.286 0.055 0.231
## 3 Retention 0.266 0.203 0.062
## 4 NPS 0.475 0.422 0.052| Construct | Satisfaction | Stress | Retention | NPS |
|---|---|---|---|---|
| Leadership | High | Weak | Modest | Very high |
| Foundations | High | High | Highest | Highest |
| LMC | Moderate | Small | Moderate | High |
| Growth | Moderate | + (stress) | High | High |
| Depletion | High (−) | Dominant | Small | Moderate |
A multivariate regression predicting Net Promoter Score (NPS) explained 47.4% of the variance. Foundational supports emerged as the strongest predictor (β = .88, p < .001), followed by responsive leadership (β = .59, p < .001) and labor–management collaboration (β = .43, p < .001). Growth orientation was positively associated with NPS (β = .36, p < .001), whereas depletion was negatively associated (β = −.19, p < .001). Together, these results indicate that educators’ willingness to recommend their institution reflects a combination of strong structural conditions, psychologically supportive leadership, institutional trust, and sustained motivational engagement. Comparing decomposed and composite models revealed that modeling specific dimensions of teacher well-being substantially improved predictive validity across outcomes. Decomposed models explained an additional 8.7% of variance in job satisfaction, 23.1% in job stress, 6.2% in retention, and 5.2% in NPS relative to composite-only models. These findings demonstrate that global well-being scores obscure distinct psychological and structural pathways linking work conditions to educator outcomes. Teacher well-being is not a single continuum from “bad” to “good.” Different aspects of well-being support different forms of functioning — emotional, behavioral, and institutional.
┌──────────────────────────────┐
│ FOUNDATIONAL SUPPORTS │
│ (Pay, Staffing, PD, Equity) │
└──────────────┬───────────────┘
│
│ (Strong for all outcomes)
▼┌──────────────────────────────┐ ┌──────────────────────────────┐ │ RESPONSIVE LEADERSHIP │────▶│ JOB SATISFACTION │ │ (Trust, Voice, Safety, │ │ (Affective Evaluation) │ │ Follow-Through) │ └──────────────────────────────┘ └──────────────┬───────────────┘ │ │ (Strong for advocacy & meaning) ▼ ┌──────────────────────────────┐ ┌──────────────────────────────┐ │ GROWTH ORIENTATION │────▶│ RETENTION INTENTIONS │ │ (Learning, Improvement) │ │ (Commitment / Staying) │ └──────────────┬───────────────┘ └──────────────────────────────┘ │ │ ▼ ┌──────────────────────────────┐ ┌──────────────────────────────┐ │ DEPLETION │────▶│ JOB STRESS │ │ (Cognitive & Emotional │ │ (Strain / Load) │ │ Exhaustion) │ └──────────────────────────────┘ └──────────────────────────────┘
┌──────────────────────────────┐ ┌──────────────────────────────┐ │ LABOR–MANAGEMENT │────▶│ NET PROMOTER SCORE (NPS) │ │ COLLABORATION (LMC) │ │ (Advocacy & Trust) │ │ (Institutional Trust) │ └──────────────────────────────┘ └──────────────────────────────┘
───────────────────────────────────────────────────────────────────── Key Insight: • Satisfaction reflects experience • Stress reflects strain • Retention reflects feasibility + commitment • NPS reflects identity, trust, and advocacy Figure 1. A multidimensional model linking teacher well-being dimensions to distinct organizational outcomes. Structural supports and leadership broadly predict satisfaction, retention, and advocacy, while depletion uniquely drives job stress. Labor–management collaboration functions as an institutional trust signal, particularly for advocacy outcomes.
A Multidimensional Model of Teacher Well-Being and Organizational Outcomes
Teacher well-being is often conceptualized as a single continuum ranging from distress to flourishing. However, such unidimensional approaches obscure the fact that different psychological and structural components of work support distinct forms of functioning. Building on stress, motivation, and organizational trust literatures, we propose a multidimensional model in which specific well-being dimensions differentially predict satisfaction, stress, retention, and institutional advocacy.
Foundational Supports as Structural Preconditions
Foundational supports—including compensation, staffing adequacy, professional development, and equity of opportunity—represent structural conditions that enable sustainable engagement. These supports consistently predicted all outcomes, particularly job satisfaction, retention, and Net Promoter Score (NPS). This pattern suggests that foundational supports function as preconditions for positive experiences rather than discretionary enhancers; when these conditions are weak, even strong leadership or individual resilience may be insufficient to sustain engagement or advocacy.
Responsive Leadership as a Psychological Experience
Responsive leadership—characterized by trust, follow-through, voice, cohesion, and psychological safety—emerged as a central predictor of satisfaction and institutional advocacy. Importantly, leadership showed weaker associations with job stress once depletion was accounted for, indicating that leadership primarily shapes how educators interpret and relate to their work rather than the objective strain they experience. This distinction clarifies why leadership is often strongly linked to morale and commitment but inconsistently linked to stress reduction.
Growth Orientation and Motivational Commitment
Growth orientation uniquely predicted retention and advocacy, even when controlling for leadership and structural supports. This suggests that growth reflects motivational commitment and meaning-making rather than comfort or ease. Growth-oriented educators appear more willing to remain in and publicly endorse their institutions, even when work is demanding—highlighting growth as a psychological resource that supports persistence rather than stress avoidance.
Depletion as the Primary Driver of Stress
Depletion—capturing cognitive and emotional exhaustion—was the dominant predictor of job stress and negatively associated with all positive outcomes. Unlike leadership or growth, depletion reflects accumulated load rather than relational or motivational processes. This finding reinforces theoretical distinctions between stress as strain and well-being as evaluation, demonstrating that reducing stress requires addressing workload and recovery rather than solely improving leadership climate.
Labor–Management Collaboration as Institutional Trust
Labor–management collaboration (LMC) contributed independently to satisfaction, retention, and especially NPS, even after accounting for leadership and well-being. This pattern positions LMC as an institutional-level trust signal, reflecting educators’ beliefs that systems for shared governance function effectively. Advocacy outcomes such as NPS appear particularly sensitive to these broader institutional dynamics.
Implications for Theory and Practice
Together, these findings demonstrate that teacher well-being is not a single latent state but a coordinated system of psychological, relational, and structural components. Different outcomes draw on different components of this system. Stress reflects depletion, satisfaction reflects experience, retention reflects feasibility and commitment, and advocacy reflects trust and identity. Models that collapse these dimensions into a single composite risk obscuring both mechanisms and intervention targets.
df %>%
summarise(across(
c(job_sat, job_stress, retention_score, nps_score,
lmc_composite_score, school_lmc_score, district_lmc_score,
twb_composite_score, twb_leader_score_short8, foundational_support_score),
list(mean = ~mean(.x, na.rm = TRUE),
sd = ~sd(.x, na.rm = TRUE))
)) %>%
pivot_longer(everything(),
names_to = c("variable","stat"),
names_sep = "_(?=[^_]+$)") %>%
pivot_wider(names_from = stat, values_from = value)## # A tibble: 10 × 3
## variable mean sd
## <chr> <dbl> <dbl>
## 1 job_sat 7.05 1.99
## 2 job_stress 7.13 2.16
## 3 retention_score 3.86 0.960
## 4 nps_score 6.60 2.64
## 5 lmc_composite_score 4.52 0.831
## 6 school_lmc_score 4.46 0.986
## 7 district_lmc_score 4.42 0.906
## 8 twb_composite_score 4.43 0.762
## 9 twb_leader_score_short8 4.55 1.21
## 10 foundational_support_score 3.47 1.09library(tidyverse)
vars <- c(
"job_sat", "job_stress", "retention_score", "nps_score",
"lmc_composite_score", "school_lmc_score", "district_lmc_score",
"twb_composite_score", "twb_leader_score_short8", "foundational_support_score"
)
var_labels <- c(
job_sat = "Job Satisfaction",
job_stress = "Job Stress",
retention_score = "Retention Intentions",
nps_score = "Net Promoter (School+District)",
lmc_composite_score = "LMC Composite",
school_lmc_score = "School-Level LMC",
district_lmc_score = "District-Level LMC",
twb_composite_score = "Teacher Well-Being (Composite)",
twb_leader_score_short8 = "TWB Leadership (8-item)",
foundational_support_score = "TWB Foundational Supports"
)
cohort_summary <- df %>%
summarise(across(all_of(vars),
list(
n = ~sum(!is.na(.x)),
mean = ~mean(.x, na.rm = TRUE),
sd = ~sd(.x, na.rm = TRUE)
),
.names = "{.col}__{.fn}"
)) %>%
pivot_longer(everything(),
names_to = c("variable","stat"),
names_sep = "__",
values_to = "value") %>%
pivot_wider(names_from = stat, values_from = value) %>%
mutate(
se = sd / sqrt(n),
ci_low = mean - 1.96 * se,
ci_high = mean + 1.96 * se,
variable_label = recode(variable, !!!var_labels)
) %>%
arrange(mean)
cohort_summary## # A tibble: 10 × 8
## variable n mean sd se ci_low ci_high variable_label
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 foundational_support… 5597 3.47 1.09 0.0145 3.44 3.50 TWB Foundatio…
## 2 retention_score 5435 3.86 0.960 0.0130 3.84 3.89 Retention Int…
## 3 district_lmc_score 5290 4.42 0.906 0.0125 4.40 4.44 District-Leve…
## 4 twb_composite_score 5987 4.43 0.762 0.00985 4.41 4.45 Teacher Well-…
## 5 school_lmc_score 5353 4.46 0.986 0.0135 4.43 4.49 School-Level …
## 6 lmc_composite_score 5454 4.52 0.831 0.0113 4.49 4.54 LMC Composite
## 7 twb_leader_score_sho… 5881 4.55 1.21 0.0158 4.52 4.58 TWB Leadershi…
## 8 nps_score 5469 6.60 2.64 0.0357 6.53 6.67 Net Promoter …
## 9 job_sat 5477 7.05 1.99 0.0268 7.00 7.10 Job Satisfact…
## 10 job_stress 5480 7.13 2.16 0.0292 7.07 7.19 Job Stressggplot(cohort_summary, aes(x = mean, y = reorder(variable_label, mean))) +
geom_errorbarh(aes(xmin = ci_low, xmax = ci_high), height = 0.15) +
geom_point(size = 2.2) +
geom_text(aes(label = paste0("n=", n)), hjust = -0.05, size = 3) +
labs(
title = "Cohort Averages Across Key Outcomes",
subtitle = "Points are means; error bars are 95% CIs. Labels show non-missing n.",
x = "Mean (original scale)",
y = NULL
) +
theme_minimal(base_size = 12) +
theme(
panel.grid.minor = element_blank(),
panel.grid.major.y = element_blank(),
axis.title.x = element_text(face = "bold"),
plot.title = element_text(face = "bold")
)
ggplot(df, aes(lmc_composite_score, job_sat)) +
geom_point(alpha = 0.2) +
geom_smooth(method = "lm", se = TRUE) +
labs(
x = "Labor-Management Collaboration (Composite)",
y = "Job Satisfaction",
title = "Higher Collaboration Is Associated with Higher Job Satisfaction"
) +
theme_minimal()
df_long <- df %>%
select(job_sat,
union_view_score,
school_lmc_score,
district_lmc_score) %>%
pivot_longer(-job_sat,
names_to = "lmc_level",
values_to = "lmc_score")
ggplot(df_long, aes(lmc_score, job_sat)) +
geom_point(alpha = 0.2) +
geom_smooth(method = "lm", se = FALSE) +
facet_wrap(~ lmc_level, scales = "free_x") +
labs(
x = "LMC Score",
y = "Job Satisfaction",
title = "Job Satisfaction by Level of Collaboration"
) +
theme_minimal()
df %>%
group_by(district) %>%
summarise(
lmc = mean(lmc_composite_score, na.rm = TRUE),
twb = mean(twb_composite_score, na.rm = TRUE),
job_sat = mean(job_sat, na.rm = TRUE),
stress = mean(job_stress, na.rm = TRUE),
retention = mean(retention_score, na.rm = TRUE)
) %>%
pivot_longer(-district) %>%
ggplot(aes(district, value, fill = district)) +
geom_col(show.legend = FALSE) +
facet_wrap(~ name, scales = "free_y") +
coord_flip() +
labs(
title = "District Differences in Collaboration and Outcomes",
x = "District",
y = "Mean Score"
) +
theme_minimal()
ggplot(df, aes(role, lmc_composite_score, fill = role)) +
geom_boxplot(alpha = 0.7, outlier.alpha = 0.2) +
coord_flip() +
labs(
title = "Labor-Management Collaboration by Role",
y = "LMC Composite",
x = "Role"
) +
theme_minimal()
df %>%
group_by(race) %>%
summarise(
n = n(),
lmc = mean(lmc_composite_score, na.rm = TRUE),
twb = mean(twb_composite_score, na.rm = TRUE),
ret = mean(retention_score, na.rm = TRUE),
sat = mean(job_sat, na.rm = TRUE),
stress = mean(job_stress, na.rm = TRUE),
nps = mean(nps_score, na.rm = TRUE)
) %>%
filter(n > 30)## # A tibble: 8 × 8
## race n lmc twb ret sat stress nps
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 American Indian/Native American or… 37 4.45 4.34 3.84 7.43 6.54 7.23
## 2 Asian 167 4.52 4.49 4.03 7.46 6.78 6.96
## 3 Black or African American 417 4.65 4.49 3.95 7.40 6.47 6.94
## 4 Other 344 4.46 4.45 3.92 7.44 6.79 7.01
## 5 Prefer not to say 759 4.20 4.10 3.58 6.37 7.49 5.47
## 6 White or Caucasian 3332 4.58 4.46 3.90 7.08 7.21 6.75
## 7 White or Caucasian,Other 32 4.71 4.71 4.08 7.78 6.62 7.62
## 8 <NA> 811 4.36 4.55 3.83 7.17 7.04 6.52df %>%
group_by(district) %>%
summarise(
n = n(),
lmc = mean(lmc_composite_score, na.rm = TRUE),
twb = mean(twb_composite_score, na.rm = TRUE),
ret = mean(retention_score, na.rm = TRUE),
sat = mean(job_sat, na.rm = TRUE),
stress = mean(job_stress, na.rm = TRUE),
nps = mean(nps_score, na.rm = TRUE)
) %>%
filter(n > 30)## # A tibble: 10 × 8
## district n lmc twb ret sat stress nps
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ABC 983 4.32 4.35 3.74 6.94 7.17 6.28
## 2 Boise 724 4.65 4.62 3.96 7.28 7.19 7.49
## 3 Burlington 404 4.42 4.46 4.01 7.35 6.35 6.39
## 4 Cherry Hill 539 4.59 4.36 3.75 6.80 7.20 5.96
## 5 Cotati-Rohnert Park 523 4.72 4.78 3.92 7.55 7.10 7.27
## 6 Dakota 970 4.52 4.31 3.91 6.76 7.59 6.82
## 7 Gateway 117 4.60 4.69 4.17 7.56 6.30 7.76
## 8 Rahway 208 4.57 4.20 3.41 6.15 7.15 4.52
## 9 Salinas 343 4.46 4.19 3.57 6.70 7.46 6.03
## 10 Trenton 1251 4.48 4.43 3.96 7.21 6.90 6.59df %>%
group_by(gender) %>%
summarise(
n = n(),
lmc = mean(lmc_composite_score, na.rm = TRUE),
twb = mean(twb_composite_score, na.rm = TRUE),
ret = mean(retention_score, na.rm = TRUE),
sat = mean(job_sat, na.rm = TRUE),
stress = mean(job_stress, na.rm = TRUE),
nps = mean(nps_score, na.rm = TRUE)
) %>%
filter(n > 30)## # A tibble: 4 × 8
## gender n lmc twb ret sat stress nps
## <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Female 3975 4.57 4.45 3.89 7.09 7.19 6.73
## 2 Male 875 4.49 4.46 3.91 7.35 6.59 6.79
## 3 Prefer not to say 422 4.13 4.01 3.47 6.11 7.62 5.16
## 4 <NA> 767 4.37 4.55 3.83 7.05 7.13 6.36# Visualization polish + model summaries
library(ggplot2)
library(scales)
library(broom)# Output folder for figures
fig_dir <- file.path(data_dir, "figures_fall_2025")
if (!dir.exists(fig_dir)) dir.create(fig_dir, recursive = TRUE)
# A clean, consistent theme
theme_pub <- function(base_size = 12) {
theme_minimal(base_size = base_size) +
theme(
panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "grey30"),
axis.title = element_text(face = "bold"),
legend.title = element_text(face = "bold")
)
}
# Consistent ggsave defaults (print-friendly)
save_pub <- function(plot_obj, filename, width = 7.5, height = 5.0) {
ggsave(
filename = file.path(fig_dir, filename),
plot = plot_obj,
width = width, height = height, dpi = 300, bg = "white"
)
}# Outcomes (numeric)
outcomes <- c("job_sat","job_stress","retention_score","nps_score")
# LMC predictors at different levels
lmc_levels <- c("union_view_score","school_lmc_score","district_lmc_score","lmc_composite_score")
# TWB predictors (final structure)
twb_factors <- c(
"twb_leader_score_short8",
"growth_score",
"wellbeing_score",
"acceptance_score",
"adaptability_score",
"depletion_score",
"foundational_support_score",
"twb_composite_score"
)
# Labels for partner-facing plots
outcome_labels <- c(
job_sat = "Job Satisfaction",
job_stress = "Job Stress",
retention_score = "Retention Intentions",
nps_score = "Net Promoter (School+District)"
)
lmc_labels <- c(
union_view_score = "Union View",
school_lmc_score = "School-Level LMC",
district_lmc_score = "District-Level LMC",
lmc_composite_score = "LMC Composite"
)# Build a long dataframe for plotting
df_plot <- df %>%
select(all_of(outcomes), all_of(lmc_levels)) %>%
pivot_longer(cols = all_of(outcomes), names_to = "outcome", values_to = "y") %>%
pivot_longer(cols = all_of(lmc_levels), names_to = "lmc_level", values_to = "lmc") %>%
mutate(
outcome_lab = outcome_labels[outcome],
lmc_lab = lmc_labels[lmc_level]
)
# Scatter + linear smooth, faceted
p_lmc_outcomes <- ggplot(df_plot, aes(x = lmc, y = y)) +
geom_point(alpha = 0.12, size = 1) +
geom_smooth(method = "lm", se = TRUE) +
facet_grid(outcome_lab ~ lmc_lab, scales = "free_y") +
labs(
title = "Labor-Management Collaboration (LMC) and Outcomes",
subtitle = "Each panel shows a linear association with robust sample coverage (point alpha reduces overplotting)",
x = "LMC Score",
y = NULL
) +
theme_pub(11)
p_lmc_outcomes
save_pub(p_lmc_outcomes, "Fig1_LMC_by_level_vs_Outcomes.png", width = 12, height = 9)desc_vars <- c(outcomes, lmc_levels, twb_factors)
desc_tbl <- df %>%
summarise(across(all_of(desc_vars),
list(
n = ~sum(!is.na(.x)),
mean = ~mean(.x, na.rm = TRUE),
sd = ~sd(.x, na.rm = TRUE)
))) %>%
pivot_longer(everything(),
names_to = c("var","stat"),
names_sep = "_(?=[^_]+$)") %>%
pivot_wider(names_from = stat, values_from = value) %>%
mutate(
se = sd / sqrt(n),
ci_low = mean - 1.96 * se,
ci_high = mean + 1.96 * se
)
# Labeling
pretty_names <- c(
job_sat = "Job Satisfaction",
job_stress = "Job Stress",
retention_score = "Retention Intentions",
nps_score = "NPS (School+District)",
union_view_score = "Union View",
school_lmc_score = "School LMC",
district_lmc_score = "District LMC",
lmc_composite_score = "LMC Composite",
twb_leader_score_short8 = "TWB Leadership (8-item)",
growth_score = "TWB Growth",
wellbeing_score = "TWB Well-being",
acceptance_score = "TWB Acceptance",
adaptability_score = "TWB Adaptability",
depletion_score = "TWB Depletion",
foundational_support_score = "TWB Foundational Supports",
twb_composite_score = "TWB Composite"
)
desc_tbl <- desc_tbl %>%
mutate(label = pretty_names[var]) %>%
arrange(desc(mean))
p_desc <- ggplot(desc_tbl, aes(x = mean, y = reorder(label, mean))) +
geom_point(size = 2) +
geom_errorbarh(aes(xmin = ci_low, xmax = ci_high), height = 0.2) +
labs(
title = "Descriptive Snapshot (Means with 95% CI)",
x = "Mean (with 95% CI)",
y = NULL
) +
theme_pub(11)
p_desc
save_pub(p_desc, "Fig2_Descriptives_Means_CI.png", width = 9, height = 7)vars_cor <- c(lmc_levels, twb_factors, outcomes)
cor_mat <- cor(df %>% select(all_of(vars_cor)), use = "pairwise.complete.obs")
cor_long <- as.data.frame(as.table(cor_mat)) %>%
rename(var1 = Var1, var2 = Var2, r = Freq) %>%
mutate(
var1_lab = pretty_names[var1],
var2_lab = pretty_names[var2]
)
p_cor <- ggplot(cor_long, aes(x = var1_lab, y = var2_lab, fill = r)) +
geom_tile() +
scale_fill_gradient2(limits = c(-1,1), breaks = c(-1,-0.5,0,0.5,1)) +
labs(
title = "Correlation Map: LMC, TWB, and Outcomes",
x = NULL, y = NULL, fill = "r"
) +
theme_pub(10) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
p_cor
save_pub(p_cor, "Fig3_Correlation_Heatmap.png", width = 11, height = 9)library(tidyverse)
# ------------------------------------------------------------
# 1) Variables included in the correlation set
# (assumes you already defined: lmc_levels, twb_factors, outcomes, pretty_names)
# ------------------------------------------------------------
vars_cor <- c(lmc_levels, twb_factors, outcomes)
# keep only variables that actually exist (prevents silent failures)
vars_cor <- vars_cor[vars_cor %in% names(df)]
# ------------------------------------------------------------
# 2) Correlation matrix (r)
# ------------------------------------------------------------
cor_mat <- cor(df %>% select(all_of(vars_cor)), use = "pairwise.complete.obs")
# Long form of r (for plotting)
cor_long <- as.data.frame(as.table(cor_mat)) %>%
rename(var1 = Var1, var2 = Var2, r = Freq) %>%
mutate(
var1 = as.character(var1),
var2 = as.character(var2),
var1_lab = pretty_names[var1],
var2_lab = pretty_names[var2]
)
# ------------------------------------------------------------
# 3) Compute p-values for each pair using cor.test
# (pairwise complete)
# ------------------------------------------------------------
pair_grid <- expand.grid(var1 = vars_cor, var2 = vars_cor, stringsAsFactors = FALSE) %>%
filter(var1 != var2)
pair_grid <- pair_grid %>%
mutate(
r = map2_dbl(var1, var2, ~ cor(df[[.x]], df[[.y]], use = "pairwise.complete.obs")),
p = map2_dbl(var1, var2, ~ cor.test(df[[.x]], df[[.y]],
use = "pairwise.complete.obs")$p.value)
) %>%
mutate(
p_adj = p.adjust(p, method = "BH"),
sig = case_when(
p_adj < .001 ~ "***",
p_adj < .01 ~ "**",
p_adj < .05 ~ "*",
TRUE ~ ""
),
var1_lab = pretty_names[var1],
var2_lab = pretty_names[var2]
)
# ------------------------------------------------------------
# 4) Table output (unique pairs only; no duplicates)
# ------------------------------------------------------------
corr_tbl <- pair_grid %>%
mutate(pair_id = map2_chr(var1, var2, ~ paste(sort(c(.x, .y)), collapse = " | "))) %>%
group_by(pair_id) %>%
slice(1) %>% # keep one direction only
ungroup() %>%
select(var1_lab, var2_lab, r, p, p_adj, sig) %>%
arrange(desc(abs(r)))
knitr::kable(
corr_tbl,
digits = 3,
caption = "Correlations (Pearson r) with BH-adjusted p-values"
)| var1_lab | var2_lab | r | p | p_adj | sig |
|---|---|---|---|---|---|
| TWB Composite | TWB Leadership (8-item) | 0.927 | 0.000 | 0.000 | *** |
| LMC Composite | District LMC | 0.910 | 0.000 | 0.000 | *** |
| LMC Composite | School LMC | 0.850 | 0.000 | 0.000 | *** |
| LMC Composite | Union View | 0.786 | 0.000 | 0.000 | *** |
| TWB Composite | TWB Foundational Supports | 0.736 | 0.000 | 0.000 | *** |
| NPS (School+District) | Job Satisfaction | 0.722 | 0.000 | 0.000 | *** |
| TWB Adaptability | TWB Acceptance | 0.649 | 0.000 | 0.000 | *** |
| District LMC | School LMC | 0.645 | 0.000 | 0.000 | *** |
| NPS (School+District) | TWB Composite | 0.632 | 0.000 | 0.000 | *** |
| Job Satisfaction | TWB Composite | 0.621 | 0.000 | 0.000 | *** |
| Job Satisfaction | TWB Foundational Supports | 0.618 | 0.000 | 0.000 | *** |
| NPS (School+District) | TWB Foundational Supports | 0.608 | 0.000 | 0.000 | *** |
| District LMC | Union View | 0.579 | 0.000 | 0.000 | *** |
| NPS (School+District) | TWB Leadership (8-item) | 0.558 | 0.000 | 0.000 | *** |
| TWB Composite | School LMC | 0.557 | 0.000 | 0.000 | *** |
| NPS (School+District) | Retention Intentions | 0.536 | 0.000 | 0.000 | *** |
| TWB Foundational Supports | TWB Leadership (8-item) | 0.536 | 0.000 | 0.000 | *** |
| TWB Leadership (8-item) | School LMC | 0.533 | 0.000 | 0.000 | *** |
| School LMC | Union View | 0.529 | 0.000 | 0.000 | *** |
| Retention Intentions | Job Satisfaction | 0.525 | 0.000 | 0.000 | *** |
| Job Satisfaction | TWB Leadership (8-item) | 0.524 | 0.000 | 0.000 | *** |
| TWB Composite | LMC Composite | 0.500 | 0.000 | 0.000 | *** |
| TWB Composite | TWB Growth | 0.493 | 0.000 | 0.000 | *** |
| TWB Foundational Supports | District LMC | 0.490 | 0.000 | 0.000 | *** |
| Job Stress | TWB Depletion | 0.485 | 0.000 | 0.000 | *** |
| Retention Intentions | TWB Foundational Supports | 0.482 | 0.000 | 0.000 | *** |
| TWB Composite | District LMC | 0.472 | 0.000 | 0.000 | *** |
| TWB Foundational Supports | LMC Composite | 0.453 | 0.000 | 0.000 | *** |
| TWB Composite | TWB Acceptance | 0.444 | 0.000 | 0.000 | *** |
| NPS (School+District) | LMC Composite | 0.443 | 0.000 | 0.000 | *** |
| NPS (School+District) | District LMC | 0.440 | 0.000 | 0.000 | *** |
| TWB Leadership (8-item) | LMC Composite | 0.440 | 0.000 | 0.000 | *** |
| Retention Intentions | TWB Composite | 0.434 | 0.000 | 0.000 | *** |
| TWB Adaptability | TWB Growth | 0.428 | 0.000 | 0.000 | *** |
| TWB Composite | TWB Adaptability | 0.422 | 0.000 | 0.000 | *** |
| TWB Foundational Supports | School LMC | 0.422 | 0.000 | 0.000 | *** |
| NPS (School+District) | School LMC | 0.420 | 0.000 | 0.000 | *** |
| Job Satisfaction | LMC Composite | 0.399 | 0.000 | 0.000 | *** |
| Job Satisfaction | School LMC | 0.397 | 0.000 | 0.000 | *** |
| TWB Depletion | TWB Well-being | -0.394 | 0.000 | 0.000 | *** |
| Job Stress | Job Satisfaction | -0.394 | 0.000 | 0.000 | *** |
| Job Satisfaction | District LMC | 0.392 | 0.000 | 0.000 | *** |
| TWB Leadership (8-item) | District LMC | 0.387 | 0.000 | 0.000 | *** |
| TWB Acceptance | TWB Growth | 0.370 | 0.000 | 0.000 | *** |
| Job Satisfaction | TWB Growth | 0.349 | 0.000 | 0.000 | *** |
| Retention Intentions | TWB Leadership (8-item) | 0.344 | 0.000 | 0.000 | *** |
| TWB Well-being | TWB Growth | 0.343 | 0.000 | 0.000 | *** |
| Job Stress | TWB Foundational Supports | -0.337 | 0.000 | 0.000 | *** |
| Job Satisfaction | TWB Depletion | -0.333 | 0.000 | 0.000 | *** |
| NPS (School+District) | TWB Growth | 0.329 | 0.000 | 0.000 | *** |
| TWB Growth | TWB Leadership (8-item) | 0.328 | 0.000 | 0.000 | *** |
| Retention Intentions | LMC Composite | 0.325 | 0.000 | 0.000 | *** |
| Retention Intentions | District LMC | 0.324 | 0.000 | 0.000 | *** |
| TWB Acceptance | TWB Well-being | 0.320 | 0.000 | 0.000 | *** |
| TWB Adaptability | TWB Well-being | 0.317 | 0.000 | 0.000 | *** |
| TWB Foundational Supports | TWB Growth | 0.316 | 0.000 | 0.000 | *** |
| Job Satisfaction | TWB Acceptance | 0.313 | 0.000 | 0.000 | *** |
| TWB Foundational Supports | TWB Acceptance | 0.310 | 0.000 | 0.000 | *** |
| TWB Depletion | TWB Acceptance | -0.296 | 0.000 | 0.000 | *** |
| NPS (School+District) | Job Stress | -0.291 | 0.000 | 0.000 | *** |
| TWB Foundational Supports | TWB Depletion | -0.289 | 0.000 | 0.000 | *** |
| Retention Intentions | School LMC | 0.284 | 0.000 | 0.000 | *** |
| Job Satisfaction | TWB Adaptability | 0.281 | 0.000 | 0.000 | *** |
| TWB Acceptance | TWB Leadership (8-item) | 0.280 | 0.000 | 0.000 | *** |
| Retention Intentions | TWB Growth | 0.275 | 0.000 | 0.000 | *** |
| TWB Foundational Supports | TWB Adaptability | 0.271 | 0.000 | 0.000 | *** |
| TWB Composite | TWB Well-being | 0.270 | 0.000 | 0.000 | *** |
| TWB Adaptability | TWB Leadership (8-item) | 0.262 | 0.000 | 0.000 | *** |
| NPS (School+District) | Union View | 0.261 | 0.000 | 0.000 | *** |
| TWB Growth | LMC Composite | 0.253 | 0.000 | 0.000 | *** |
| TWB Growth | District LMC | 0.251 | 0.000 | 0.000 | *** |
| Job Stress | TWB Acceptance | -0.247 | 0.000 | 0.000 | *** |
| TWB Depletion | TWB Adaptability | -0.241 | 0.000 | 0.000 | *** |
| NPS (School+District) | TWB Acceptance | 0.241 | 0.000 | 0.000 | *** |
| TWB Growth | School LMC | 0.238 | 0.000 | 0.000 | *** |
| TWB Composite | Union View | 0.237 | 0.000 | 0.000 | *** |
| Job Satisfaction | TWB Well-being | 0.236 | 0.000 | 0.000 | *** |
| NPS (School+District) | TWB Adaptability | 0.236 | 0.000 | 0.000 | *** |
| Job Stress | TWB Composite | -0.235 | 0.000 | 0.000 | *** |
| NPS (School+District) | TWB Depletion | -0.229 | 0.000 | 0.000 | *** |
| Retention Intentions | Job Stress | -0.220 | 0.000 | 0.000 | *** |
| TWB Foundational Supports | Union View | 0.218 | 0.000 | 0.000 | *** |
| Job Satisfaction | Union View | 0.217 | 0.000 | 0.000 | *** |
| Job Stress | TWB Well-being | -0.208 | 0.000 | 0.000 | *** |
| Retention Intentions | Union View | 0.205 | 0.000 | 0.000 | *** |
| TWB Acceptance | School LMC | 0.204 | 0.000 | 0.000 | *** |
| Retention Intentions | TWB Acceptance | 0.202 | 0.000 | 0.000 | *** |
| TWB Foundational Supports | TWB Well-being | 0.201 | 0.000 | 0.000 | *** |
| Job Stress | TWB Adaptability | -0.195 | 0.000 | 0.000 | *** |
| TWB Leadership (8-item) | Union View | 0.194 | 0.000 | 0.000 | *** |
| TWB Acceptance | District LMC | 0.189 | 0.000 | 0.000 | *** |
| TWB Adaptability | District LMC | 0.189 | 0.000 | 0.000 | *** |
| TWB Adaptability | LMC Composite | 0.188 | 0.000 | 0.000 | *** |
| TWB Acceptance | LMC Composite | 0.188 | 0.000 | 0.000 | *** |
| Retention Intentions | TWB Depletion | -0.187 | 0.000 | 0.000 | *** |
| TWB Adaptability | School LMC | 0.185 | 0.000 | 0.000 | *** |
| TWB Depletion | TWB Growth | -0.171 | 0.000 | 0.000 | *** |
| Job Stress | TWB Leadership (8-item) | -0.170 | 0.000 | 0.000 | *** |
| Retention Intentions | TWB Adaptability | 0.163 | 0.000 | 0.000 | *** |
| TWB Growth | Union View | 0.156 | 0.000 | 0.000 | *** |
| TWB Composite | TWB Depletion | -0.155 | 0.000 | 0.000 | *** |
| TWB Well-being | LMC Composite | 0.151 | 0.000 | 0.000 | *** |
| TWB Well-being | District LMC | 0.140 | 0.000 | 0.000 | *** |
| NPS (School+District) | TWB Well-being | 0.135 | 0.000 | 0.000 | *** |
| Job Stress | District LMC | -0.135 | 0.000 | 0.000 | *** |
| Retention Intentions | TWB Well-being | 0.132 | 0.000 | 0.000 | *** |
| TWB Well-being | School LMC | 0.123 | 0.000 | 0.000 | *** |
| TWB Depletion | TWB Leadership (8-item) | -0.120 | 0.000 | 0.000 | *** |
| TWB Well-being | Union View | 0.120 | 0.000 | 0.000 | *** |
| TWB Well-being | TWB Leadership (8-item) | 0.118 | 0.000 | 0.000 | *** |
| TWB Depletion | District LMC | -0.113 | 0.000 | 0.000 | *** |
| Job Stress | TWB Growth | -0.109 | 0.000 | 0.000 | *** |
| TWB Adaptability | Union View | 0.108 | 0.000 | 0.000 | *** |
| Job Stress | LMC Composite | -0.106 | 0.000 | 0.000 | *** |
| Job Stress | School LMC | -0.101 | 0.000 | 0.000 | *** |
| TWB Depletion | LMC Composite | -0.099 | 0.000 | 0.000 | *** |
| TWB Depletion | School LMC | -0.092 | 0.000 | 0.000 | *** |
| TWB Acceptance | Union View | 0.089 | 0.000 | 0.000 | *** |
| TWB Depletion | Union View | -0.041 | 0.003 | 0.003 | ** |
| Job Stress | Union View | -0.023 | 0.086 | 0.086 |
# ------------------------------------------------------------
# 5) Heatmap (option A: plain tiles, like your original)
# ------------------------------------------------------------
p_cor <- ggplot(cor_long, aes(x = var1_lab, y = var2_lab, fill = r)) +
geom_tile(color = "white") +
scale_fill_gradient2(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
labs(
title = "Correlation Map: LMC, TWB, and Outcomes",
subtitle = "Pearson correlations (r). See table below for BH-adjusted significance.",
x = NULL, y = NULL, fill = "r"
) +
theme_pub(10) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
p_cor
# ------------------------------------------------------------
# 6) OPTIONAL: Heatmap with numbers + significance stars in tiles
# (comment out if too busy)
# ------------------------------------------------------------
cor_plot_tbl <- pair_grid %>%
mutate(label = paste0(round(r, 2), sig))
p_cor_labeled <- ggplot(cor_plot_tbl, aes(x = var1_lab, y = var2_lab, fill = r)) +
geom_tile(color = "white") +
geom_text(aes(label = label), size = 3) +
scale_fill_gradient2(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1)) +
labs(
title = "Correlation Map (Labeled): LMC, TWB, and Outcomes",
subtitle = "Labels show r with BH-adjusted significance: * p<.05, ** p<.01, *** p<.001",
x = NULL, y = NULL, fill = "r"
) +
theme_pub(10) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
p_cor_labeled
# Save (if you have save_pub defined)
# save_pub(p_cor, "Fig3_Correlation_Heatmap.png", width = 11, height = 9)
# save_pub(p_cor_labeled, "Fig3_Correlation_Heatmap_Labeled.png", width = 11, height = 9)library(tidyverse)
# ============================================================
# Correlation heatmap + significance table
# - Drops rows with any missing across vars_cor (complete cases)
# - Orders variables into an interpretable block structure
# - Produces (1) readable heatmap (upper triangle only)
# - Produces (2) sorted table of correlations with BH-adjusted p
# ============================================================
# ---------------------------
# 1) Define variable blocks (edit order here)
# ---------------------------
lmc_block <- c("union_view_score", "school_lmc_score", "district_lmc_score", "lmc_composite_score")
twb_block <- c(
"twb_leader_score_short8",
"growth_score",
"acceptance_score",
"adaptability_score",
"wellbeing_score",
"depletion_score",
"foundational_support_score",
"twb_composite_score"
)
outcome_block <- c("job_sat", "job_stress", "retention_score", "nps_score")
# Combine and keep only those that exist
vars_cor <- c(lmc_block, twb_block, outcome_block)
vars_cor <- vars_cor[vars_cor %in% names(df)]
# ---------------------------
# 2) Complete-case dataframe (removes NAs from *everything* used)
# ---------------------------
df_cc <- df %>%
select(all_of(vars_cor)) %>%
drop_na()
# quick check
cat("Complete-case N used for correlations:", nrow(df_cc), "\n")## Complete-case N used for correlations: 5140# ---------------------------
# 3) Correlations + p-values
# ---------------------------
# matrix of r
cor_mat <- cor(df_cc, use = "everything")
# helper to compute p for a single pair
get_p <- function(x, y) cor.test(x, y)$p.value
# all unique pairs (upper triangle only)
pairs <- expand.grid(var1 = vars_cor, var2 = vars_cor, stringsAsFactors = FALSE) %>%
mutate(i = match(var1, vars_cor),
j = match(var2, vars_cor)) %>%
filter(i < j)
pairs <- pairs %>%
mutate(
r = map2_dbl(var1, var2, ~ cor(df_cc[[.x]], df_cc[[.y]], use = "everything")),
p = map2_dbl(var1, var2, ~ get_p(df_cc[[.x]], df_cc[[.y]]))
) %>%
mutate(
p_adj = p.adjust(p, method = "BH"),
sig = case_when(
p_adj < .001 ~ "***",
p_adj < .01 ~ "**",
p_adj < .05 ~ "*",
TRUE ~ ""
),
var1_lab = if ("pretty_names" %in% ls()) pretty_names[var1] else var1,
var2_lab = if ("pretty_names" %in% ls()) pretty_names[var2] else var2
)
# ---------------------------
# 4) Table: strongest correlations first
# ---------------------------
corr_tbl <- pairs %>%
select(var1_lab, var2_lab, r, p, p_adj, sig) %>%
arrange(desc(abs(r)))
knitr::kable(
corr_tbl,
digits = 3,
caption = "Correlations (Pearson r) using complete cases; BH-adjusted p-values"
)| var1_lab | var2_lab | r | p | p_adj | sig |
|---|---|---|---|---|---|
| twb_leader_score_short8 | twb_composite_score | 0.929 | 0.000 | 0.000 | *** |
| district_lmc_score | lmc_composite_score | 0.910 | 0.000 | 0.000 | *** |
| school_lmc_score | lmc_composite_score | 0.849 | 0.000 | 0.000 | *** |
| union_view_score | lmc_composite_score | 0.775 | 0.000 | 0.000 | *** |
| foundational_support_score | twb_composite_score | 0.731 | 0.000 | 0.000 | *** |
| job_sat | nps_score | 0.725 | 0.000 | 0.000 | *** |
| school_lmc_score | district_lmc_score | 0.645 | 0.000 | 0.000 | *** |
| acceptance_score | adaptability_score | 0.643 | 0.000 | 0.000 | *** |
| twb_composite_score | nps_score | 0.633 | 0.000 | 0.000 | *** |
| twb_composite_score | job_sat | 0.621 | 0.000 | 0.000 | *** |
| foundational_support_score | job_sat | 0.619 | 0.000 | 0.000 | *** |
| foundational_support_score | nps_score | 0.610 | 0.000 | 0.000 | *** |
| union_view_score | district_lmc_score | 0.579 | 0.000 | 0.000 | *** |
| school_lmc_score | twb_composite_score | 0.564 | 0.000 | 0.000 | *** |
| twb_leader_score_short8 | nps_score | 0.557 | 0.000 | 0.000 | *** |
| retention_score | nps_score | 0.537 | 0.000 | 0.000 | *** |
| school_lmc_score | twb_leader_score_short8 | 0.536 | 0.000 | 0.000 | *** |
| twb_leader_score_short8 | foundational_support_score | 0.528 | 0.000 | 0.000 | *** |
| union_view_score | school_lmc_score | 0.526 | 0.000 | 0.000 | *** |
| job_sat | retention_score | 0.526 | 0.000 | 0.000 | *** |
| twb_leader_score_short8 | job_sat | 0.524 | 0.000 | 0.000 | *** |
| lmc_composite_score | twb_composite_score | 0.523 | 0.000 | 0.000 | *** |
| district_lmc_score | foundational_support_score | 0.496 | 0.000 | 0.000 | *** |
| foundational_support_score | retention_score | 0.485 | 0.000 | 0.000 | *** |
| growth_score | twb_composite_score | 0.483 | 0.000 | 0.000 | *** |
| depletion_score | job_stress | 0.480 | 0.000 | 0.000 | *** |
| lmc_composite_score | foundational_support_score | 0.477 | 0.000 | 0.000 | *** |
| district_lmc_score | twb_composite_score | 0.477 | 0.000 | 0.000 | *** |
| lmc_composite_score | nps_score | 0.458 | 0.000 | 0.000 | *** |
| lmc_composite_score | twb_leader_score_short8 | 0.456 | 0.000 | 0.000 | *** |
| district_lmc_score | nps_score | 0.443 | 0.000 | 0.000 | *** |
| acceptance_score | twb_composite_score | 0.435 | 0.000 | 0.000 | *** |
| twb_composite_score | retention_score | 0.432 | 0.000 | 0.000 | *** |
| school_lmc_score | foundational_support_score | 0.430 | 0.000 | 0.000 | *** |
| growth_score | adaptability_score | 0.423 | 0.000 | 0.000 | *** |
| school_lmc_score | nps_score | 0.421 | 0.000 | 0.000 | *** |
| lmc_composite_score | job_sat | 0.414 | 0.000 | 0.000 | *** |
| adaptability_score | twb_composite_score | 0.408 | 0.000 | 0.000 | *** |
| school_lmc_score | job_sat | 0.401 | 0.000 | 0.000 | *** |
| wellbeing_score | depletion_score | -0.399 | 0.000 | 0.000 | *** |
| district_lmc_score | job_sat | 0.395 | 0.000 | 0.000 | *** |
| job_sat | job_stress | -0.394 | 0.000 | 0.000 | *** |
| district_lmc_score | twb_leader_score_short8 | 0.392 | 0.000 | 0.000 | *** |
| growth_score | acceptance_score | 0.370 | 0.000 | 0.000 | *** |
| growth_score | job_sat | 0.347 | 0.000 | 0.000 | *** |
| twb_leader_score_short8 | retention_score | 0.341 | 0.000 | 0.000 | *** |
| foundational_support_score | job_stress | -0.339 | 0.000 | 0.000 | *** |
| growth_score | wellbeing_score | 0.338 | 0.000 | 0.000 | *** |
| lmc_composite_score | retention_score | 0.330 | 0.000 | 0.000 | *** |
| depletion_score | job_sat | -0.330 | 0.000 | 0.000 | *** |
| district_lmc_score | retention_score | 0.326 | 0.000 | 0.000 | *** |
| growth_score | nps_score | 0.324 | 0.000 | 0.000 | *** |
| twb_leader_score_short8 | growth_score | 0.316 | 0.000 | 0.000 | *** |
| adaptability_score | wellbeing_score | 0.315 | 0.000 | 0.000 | *** |
| acceptance_score | wellbeing_score | 0.314 | 0.000 | 0.000 | *** |
| acceptance_score | job_sat | 0.313 | 0.000 | 0.000 | *** |
| growth_score | foundational_support_score | 0.311 | 0.000 | 0.000 | *** |
| acceptance_score | foundational_support_score | 0.301 | 0.000 | 0.000 | *** |
| acceptance_score | depletion_score | -0.297 | 0.000 | 0.000 | *** |
| job_stress | nps_score | -0.295 | 0.000 | 0.000 | *** |
| depletion_score | foundational_support_score | -0.288 | 0.000 | 0.000 | *** |
| school_lmc_score | retention_score | 0.285 | 0.000 | 0.000 | *** |
| twb_leader_score_short8 | acceptance_score | 0.274 | 0.000 | 0.000 | *** |
| adaptability_score | job_sat | 0.271 | 0.000 | 0.000 | *** |
| growth_score | retention_score | 0.269 | 0.000 | 0.000 | *** |
| union_view_score | nps_score | 0.269 | 0.000 | 0.000 | *** |
| lmc_composite_score | growth_score | 0.262 | 0.000 | 0.000 | *** |
| adaptability_score | foundational_support_score | 0.259 | 0.000 | 0.000 | *** |
| wellbeing_score | twb_composite_score | 0.259 | 0.000 | 0.000 | *** |
| district_lmc_score | growth_score | 0.252 | 0.000 | 0.000 | *** |
| twb_leader_score_short8 | adaptability_score | 0.251 | 0.000 | 0.000 | *** |
| union_view_score | twb_composite_score | 0.246 | 0.000 | 0.000 | *** |
| adaptability_score | depletion_score | -0.243 | 0.000 | 0.000 | *** |
| acceptance_score | job_stress | -0.242 | 0.000 | 0.000 | *** |
| school_lmc_score | growth_score | 0.238 | 0.000 | 0.000 | *** |
| acceptance_score | nps_score | 0.237 | 0.000 | 0.000 | *** |
| twb_composite_score | job_stress | -0.234 | 0.000 | 0.000 | *** |
| wellbeing_score | job_sat | 0.233 | 0.000 | 0.000 | *** |
| depletion_score | nps_score | -0.230 | 0.000 | 0.000 | *** |
| union_view_score | foundational_support_score | 0.230 | 0.000 | 0.000 | *** |
| adaptability_score | nps_score | 0.226 | 0.000 | 0.000 | *** |
| union_view_score | job_sat | 0.225 | 0.000 | 0.000 | *** |
| job_stress | retention_score | -0.219 | 0.000 | 0.000 | *** |
| school_lmc_score | acceptance_score | 0.209 | 0.000 | 0.000 | *** |
| wellbeing_score | job_stress | -0.207 | 0.000 | 0.000 | *** |
| union_view_score | retention_score | 0.205 | 0.000 | 0.000 | *** |
| lmc_composite_score | acceptance_score | 0.201 | 0.000 | 0.000 | *** |
| union_view_score | twb_leader_score_short8 | 0.200 | 0.000 | 0.000 | *** |
| acceptance_score | retention_score | 0.200 | 0.000 | 0.000 | *** |
| lmc_composite_score | adaptability_score | 0.196 | 0.000 | 0.000 | *** |
| wellbeing_score | foundational_support_score | 0.196 | 0.000 | 0.000 | *** |
| adaptability_score | job_stress | -0.190 | 0.000 | 0.000 | *** |
| district_lmc_score | acceptance_score | 0.188 | 0.000 | 0.000 | *** |
| district_lmc_score | adaptability_score | 0.187 | 0.000 | 0.000 | *** |
| school_lmc_score | adaptability_score | 0.185 | 0.000 | 0.000 | *** |
| depletion_score | retention_score | -0.185 | 0.000 | 0.000 | *** |
| twb_leader_score_short8 | job_stress | -0.168 | 0.000 | 0.000 | *** |
| growth_score | depletion_score | -0.167 | 0.000 | 0.000 | *** |
| union_view_score | growth_score | 0.160 | 0.000 | 0.000 | *** |
| depletion_score | twb_composite_score | -0.156 | 0.000 | 0.000 | *** |
| adaptability_score | retention_score | 0.156 | 0.000 | 0.000 | *** |
| lmc_composite_score | wellbeing_score | 0.153 | 0.000 | 0.000 | *** |
| district_lmc_score | wellbeing_score | 0.140 | 0.000 | 0.000 | *** |
| district_lmc_score | job_stress | -0.138 | 0.000 | 0.000 | *** |
| wellbeing_score | nps_score | 0.133 | 0.000 | 0.000 | *** |
| wellbeing_score | retention_score | 0.127 | 0.000 | 0.000 | *** |
| school_lmc_score | wellbeing_score | 0.126 | 0.000 | 0.000 | *** |
| lmc_composite_score | job_stress | -0.122 | 0.000 | 0.000 | *** |
| twb_leader_score_short8 | depletion_score | -0.118 | 0.000 | 0.000 | *** |
| union_view_score | wellbeing_score | 0.118 | 0.000 | 0.000 | *** |
| district_lmc_score | depletion_score | -0.113 | 0.000 | 0.000 | *** |
| union_view_score | adaptability_score | 0.113 | 0.000 | 0.000 | *** |
| school_lmc_score | job_stress | -0.111 | 0.000 | 0.000 | *** |
| twb_leader_score_short8 | wellbeing_score | 0.109 | 0.000 | 0.000 | *** |
| growth_score | job_stress | -0.109 | 0.000 | 0.000 | *** |
| lmc_composite_score | depletion_score | -0.107 | 0.000 | 0.000 | *** |
| union_view_score | acceptance_score | 0.098 | 0.000 | 0.000 | *** |
| school_lmc_score | depletion_score | -0.097 | 0.000 | 0.000 | *** |
| union_view_score | depletion_score | -0.043 | 0.002 | 0.002 | ** |
| union_view_score | job_stress | -0.032 | 0.022 | 0.022 | * |
# ---------------------------
# 5) Heatmap data (upper triangle only, ordered)
# ---------------------------
# Build full long matrix for plotting; blank lower triangle for readability
cor_long <- as.data.frame(as.table(cor_mat)) %>%
rename(var1 = Var1, var2 = Var2, r = Freq) %>%
mutate(
var1 = as.character(var1),
var2 = as.character(var2),
i = match(var1, vars_cor),
j = match(var2, vars_cor),
# keep upper triangle incl diagonal; drop lower triangle
r_plot = ifelse(i <= j, r, NA_real_),
var1 = factor(var1, levels = vars_cor),
var2 = factor(var2, levels = vars_cor),
var1_lab = if ("pretty_names" %in% ls()) pretty_names[as.character(var1)] else as.character(var1),
var2_lab = if ("pretty_names" %in% ls()) pretty_names[as.character(var2)] else as.character(var2)
)
# Add stars to matching pairs (upper triangle only)
stars_df <- pairs %>%
mutate(
var1 = factor(var1, levels = vars_cor),
var2 = factor(var2, levels = vars_cor),
label = paste0(round(r, 2), sig),
var1_lab = if ("pretty_names" %in% ls()) pretty_names[as.character(var1)] else as.character(var1),
var2_lab = if ("pretty_names" %in% ls()) pretty_names[as.character(var2)] else as.character(var2)
)
# ---------------------------
# 6) Heatmap (readable, labeled)
# ---------------------------
p_cor <- ggplot(cor_long, aes(x = var1_lab, y = var2_lab, fill = r_plot)) +
geom_tile(color = "white") +
geom_text(
data = stars_df,
aes(x = var1_lab, y = var2_lab, label = label),
inherit.aes = FALSE,
size = 3
) +
scale_fill_gradient2(limits = c(-1, 1), breaks = c(-1, -0.5, 0, 0.5, 1), na.value = "white") +
labs(
title = "Correlation Map: LMC, TWB, and Outcomes (Complete Cases)",
subtitle = "Upper triangle shown. Labels are r with BH-adjusted significance (* p<.05, ** p<.01, *** p<.001).",
x = NULL, y = NULL, fill = "r"
) +
theme_pub(10) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
panel.grid = element_blank()
)
p_cor
# Optional save (if you have save_pub defined)
# save_pub(p_cor, "Fig_Correlation_Heatmap_Ordered_CC.png", width = 12, height = 9)# Key variables to test by district
vars_district <- c(outcomes, lmc_levels, twb_factors)
district_tests <- lapply(vars_district, function(v) {
f <- as.formula(paste0(v, " ~ district"))
m <- lm(f, data = df)
a <- anova(m)
tibble(
variable = v,
F = a$`F value`[1],
df1 = a$Df[1],
df2 = a$Df[2],
p = a$`Pr(>F)`[1],
r2 = summary(m)$r.squared
)
})
district_tests_tbl <- bind_rows(district_tests) %>%
mutate(
variable_label = pretty_names[variable],
p_adj = p.adjust(p, method = "BH")
) %>%
arrange(p_adj)
district_tests_tbl %>% mutate(across(c(F,p,p_adj,r2), ~round(.x, 4)))## # A tibble: 16 × 8
## variable F df1 df2 p r2 variable_label p_adj
## <chr> <dbl> <int> <int> <dbl> <dbl> <chr> <dbl>
## 1 foundational_support_sco… 46.9 9 5587 0 0.0702 TWB Foundatio… 0
## 2 nps_score 35.5 9 5459 0 0.0553 NPS (School+D… 0
## 3 twb_composite_score 30.4 9 5977 0 0.0437 TWB Composite 0
## 4 twb_leader_score_short8 28.0 9 5871 0 0.0411 TWB Leadershi… 0
## 5 district_lmc_score 17.8 9 5280 0 0.0295 District LMC 0
## 6 school_lmc_score 16.3 9 5343 0 0.0267 School LMC 0
## 7 job_sat 15.9 9 5467 0 0.0254 Job Satisfact… 0
## 8 depletion_score 15.8 9 5610 0 0.0247 TWB Depletion 0
## 9 union_view_score 15.6 9 5437 0 0.0252 Union View 0
## 10 retention_score 14.8 9 5425 0 0.024 Retention Int… 0
## 11 job_stress 14.2 9 5470 0 0.0228 Job Stress 0
## 12 lmc_composite_score 11.9 9 5444 0 0.0192 LMC Composite 0
## 13 acceptance_score 11.9 9 5656 0 0.0185 TWB Acceptance 0
## 14 growth_score 9.53 9 5735 0 0.0147 TWB Growth 0
## 15 adaptability_score 6.70 9 5656 0 0.0105 TWB Adaptabil… 0
## 16 wellbeing_score 3.71 9 5727 1e-4 0.0058 TWB Well-being 1e-4library(tidyverse)
# ============================================================
# District "who is better/worse?" dashboard (per outcome)
# 1) Overall district effect (ANOVA + BH p_adj)
# 2) District means + 95% CI + sample size
# 3) Pairwise district comparisons vs GRAND MEAN (one-sample t on centered scores)
# - This avoids O(K^2) pairwise district tests
# - Flags districts significantly above/below overall mean for each outcome
# - BH-adjusts p within each outcome
#
# Notes:
# - Uses complete-case within each outcome separately (keeps more data).
# - Treats district as fixed effect (descriptive inference).
# ============================================================
# ---- choose outcomes/vars you want to audit by district
vars_district <- c(
"job_sat","job_stress","retention_score","nps_score",
"lmc_composite_score","school_lmc_score","district_lmc_score","union_view_score",
"twb_composite_score","twb_leader_score_short8","foundational_support_score"
)
vars_district <- vars_district[vars_district %in% names(df)]
min_n_dist <- 30 # set your preferred threshold
# Helper: district summaries with CI
district_summary_ci <- function(dat, var, min_n = 30) {
dat %>%
select(district, all_of(var)) %>%
drop_na() %>%
group_by(district) %>%
summarise(
n = n(),
mean = mean(.data[[var]], na.rm = TRUE),
sd = sd(.data[[var]], na.rm = TRUE),
se = sd / sqrt(n),
ci_low = mean - 1.96 * se,
ci_high = mean + 1.96 * se,
.groups = "drop"
) %>%
filter(n >= min_n) %>%
mutate(variable = var)
}
# Helper: overall district effect test (ANOVA from lm)
district_anova <- function(dat, var) {
d <- dat %>% select(district, all_of(var)) %>% drop_na()
if (nrow(d) < 10) return(tibble(variable = var, F = NA, df1 = NA, df2 = NA, p = NA, r2 = NA, n = nrow(d)))
m <- lm(reformulate("district", response = var), data = d)
a <- anova(m)
tibble(
variable = var,
F = a$`F value`[1],
df1 = a$Df[1],
df2 = a$Df[2],
p = a$`Pr(>F)`[1],
r2 = summary(m)$r.squared,
n = nrow(d)
)
}
# Helper: district vs grand mean test (center outcome; one-sample t within district)
district_vs_grandmean <- function(dat, var, min_n = 30) {
d <- dat %>% select(district, all_of(var)) %>% drop_na()
grand_mean <- mean(d[[var]], na.rm = TRUE)
d %>%
mutate(centered = .data[[var]] - grand_mean) %>%
group_by(district) %>%
summarise(
n = n(),
mean = mean(.data[[var]], na.rm = TRUE),
grand_mean = grand_mean,
diff = mean - grand_mean,
sd = sd(centered, na.rm = TRUE),
se = sd / sqrt(n),
t = diff / se,
p = 2 * pt(abs(t), df = n - 1, lower.tail = FALSE),
.groups = "drop"
) %>%
filter(n >= min_n) %>%
mutate(
variable = var,
p_adj = p.adjust(p, method = "BH"),
flag = case_when(
p_adj < .05 & diff > 0 ~ "Above grand mean",
p_adj < .05 & diff < 0 ~ "Below grand mean",
TRUE ~ "Not different"
)
)
}
# ---------------------------
# 1) Overall ANOVA tests by district
# ---------------------------
district_tests_tbl <- map_dfr(vars_district, ~ district_anova(df, .x)) %>%
mutate(
variable_label = if ("pretty_names" %in% ls()) pretty_names[variable] else variable,
p_adj = p.adjust(p, method = "BH")
) %>%
arrange(p_adj)
knitr::kable(
district_tests_tbl %>% mutate(across(c(F, p, p_adj, r2), ~ round(.x, 4))),
caption = paste0("District effects (ANOVA via lm). BH-adjusted across variables. min district n = ", min_n_dist)
)| variable | F | df1 | df2 | p | r2 | n | variable_label | p_adj |
|---|---|---|---|---|---|---|---|---|
| foundational_support_score | 46.8506 | 9 | 5587 | 0 | 0.0702 | 5597 | foundational_support_score | 0 |
| nps_score | 35.5180 | 9 | 5459 | 0 | 0.0553 | 5469 | nps_score | 0 |
| twb_composite_score | 30.3564 | 9 | 5977 | 0 | 0.0437 | 5987 | twb_composite_score | 0 |
| twb_leader_score_short8 | 27.9757 | 9 | 5871 | 0 | 0.0411 | 5881 | twb_leader_score_short8 | 0 |
| district_lmc_score | 17.8183 | 9 | 5280 | 0 | 0.0295 | 5290 | district_lmc_score | 0 |
| school_lmc_score | 16.2813 | 9 | 5343 | 0 | 0.0267 | 5353 | school_lmc_score | 0 |
| job_sat | 15.8530 | 9 | 5467 | 0 | 0.0254 | 5477 | job_sat | 0 |
| union_view_score | 15.6380 | 9 | 5437 | 0 | 0.0252 | 5447 | union_view_score | 0 |
| retention_score | 14.8304 | 9 | 5425 | 0 | 0.0240 | 5435 | retention_score | 0 |
| job_stress | 14.2121 | 9 | 5470 | 0 | 0.0228 | 5480 | job_stress | 0 |
| lmc_composite_score | 11.8674 | 9 | 5444 | 0 | 0.0192 | 5454 | lmc_composite_score | 0 |
# ---------------------------
# 2) District means + CI for each variable (for plotting / ranking)
# ---------------------------
district_means_tbl <- map_dfr(vars_district, ~ district_summary_ci(df, .x, min_n = min_n_dist)) %>%
mutate(
variable_label = if ("pretty_names" %in% ls()) pretty_names[variable] else variable
)
# Example: show top/bottom districts for each variable
district_rank_tbl <- district_means_tbl %>%
group_by(variable_label) %>%
arrange(desc(mean), .by_group = TRUE) %>%
mutate(rank = row_number()) %>%
ungroup()
# Display: top 5 and bottom 5 per outcome
top_bottom_tbl <- district_rank_tbl %>%
group_by(variable_label) %>%
filter(rank <= 5 | rank > (max(rank) - 5)) %>%
arrange(variable_label, rank)
knitr::kable(
top_bottom_tbl %>% select(variable_label, district, n, mean, ci_low, ci_high, rank) %>%
mutate(across(c(mean, ci_low, ci_high), ~ round(.x, 3))),
caption = "Top 5 and Bottom 5 districts by mean (with 95% CI), filtered to districts with adequate n"
)| variable_label | district | n | mean | ci_low | ci_high | rank |
|---|---|---|---|---|---|---|
| district_lmc_score | Cotati-Rohnert Park | 477 | 4.704 | 4.636 | 4.772 | 1 |
| district_lmc_score | Boise | 656 | 4.658 | 4.591 | 4.725 | 2 |
| district_lmc_score | Gateway | 89 | 4.606 | 4.420 | 4.793 | 3 |
| district_lmc_score | Rahway | 178 | 4.463 | 4.341 | 4.585 | 4 |
| district_lmc_score | Dakota | 890 | 4.437 | 4.387 | 4.487 | 5 |
| district_lmc_score | Cherry Hill | 483 | 4.397 | 4.314 | 4.480 | 6 |
| district_lmc_score | Trenton | 1116 | 4.338 | 4.280 | 4.396 | 7 |
| district_lmc_score | Burlington | 358 | 4.325 | 4.229 | 4.421 | 8 |
| district_lmc_score | ABC | 749 | 4.248 | 4.180 | 4.316 | 9 |
| district_lmc_score | Salinas | 294 | 4.190 | 4.085 | 4.294 | 10 |
| foundational_support_score | Gateway | 115 | 4.400 | 4.247 | 4.552 | 1 |
| foundational_support_score | Cotati-Rohnert Park | 501 | 3.884 | 3.799 | 3.969 | 2 |
| foundational_support_score | Burlington | 375 | 3.752 | 3.640 | 3.865 | 3 |
| foundational_support_score | Boise | 682 | 3.581 | 3.505 | 3.656 | 4 |
| foundational_support_score | Trenton | 1162 | 3.573 | 3.507 | 3.639 | 5 |
| foundational_support_score | ABC | 848 | 3.377 | 3.302 | 3.452 | 6 |
| foundational_support_score | Dakota | 910 | 3.328 | 3.267 | 3.389 | 7 |
| foundational_support_score | Cherry Hill | 503 | 3.196 | 3.103 | 3.290 | 8 |
| foundational_support_score | Rahway | 184 | 2.973 | 2.828 | 3.118 | 9 |
| foundational_support_score | Salinas | 317 | 2.854 | 2.742 | 2.967 | 10 |
| job_sat | Gateway | 93 | 7.559 | 7.173 | 7.946 | 1 |
| job_sat | Cotati-Rohnert Park | 499 | 7.553 | 7.409 | 7.697 | 2 |
| job_sat | Burlington | 371 | 7.345 | 7.149 | 7.541 | 3 |
| job_sat | Boise | 677 | 7.276 | 7.138 | 7.414 | 4 |
| job_sat | Trenton | 1152 | 7.214 | 7.092 | 7.337 | 5 |
| job_sat | ABC | 792 | 6.942 | 6.797 | 7.087 | 6 |
| job_sat | Cherry Hill | 499 | 6.800 | 6.623 | 6.976 | 7 |
| job_sat | Dakota | 905 | 6.760 | 6.640 | 6.880 | 8 |
| job_sat | Salinas | 308 | 6.701 | 6.471 | 6.931 | 9 |
| job_sat | Rahway | 181 | 6.149 | 5.859 | 6.439 | 10 |
| job_stress | Dakota | 906 | 7.591 | 7.475 | 7.706 | 1 |
| job_stress | Salinas | 308 | 7.461 | 7.221 | 7.701 | 2 |
| job_stress | Cherry Hill | 500 | 7.204 | 7.031 | 7.377 | 3 |
| job_stress | Boise | 677 | 7.195 | 7.046 | 7.344 | 4 |
| job_stress | ABC | 792 | 7.168 | 7.009 | 7.327 | 5 |
| job_stress | Rahway | 182 | 7.154 | 6.824 | 7.484 | 6 |
| job_stress | Cotati-Rohnert Park | 499 | 7.096 | 6.917 | 7.275 | 7 |
| job_stress | Trenton | 1151 | 6.904 | 6.767 | 7.042 | 8 |
| job_stress | Burlington | 372 | 6.347 | 6.111 | 6.583 | 9 |
| job_stress | Gateway | 93 | 6.301 | 5.831 | 6.771 | 10 |
| lmc_composite_score | Cotati-Rohnert Park | 493 | 4.719 | 4.650 | 4.787 | 1 |
| lmc_composite_score | Boise | 675 | 4.647 | 4.584 | 4.710 | 2 |
| lmc_composite_score | Gateway | 92 | 4.598 | 4.400 | 4.796 | 3 |
| lmc_composite_score | Cherry Hill | 496 | 4.591 | 4.521 | 4.662 | 4 |
| lmc_composite_score | Rahway | 182 | 4.567 | 4.456 | 4.679 | 5 |
| lmc_composite_score | Dakota | 904 | 4.516 | 4.473 | 4.558 | 6 |
| lmc_composite_score | Trenton | 1150 | 4.485 | 4.433 | 4.537 | 7 |
| lmc_composite_score | Salinas | 308 | 4.465 | 4.380 | 4.550 | 8 |
| lmc_composite_score | Burlington | 372 | 4.422 | 4.337 | 4.508 | 9 |
| lmc_composite_score | ABC | 782 | 4.315 | 4.250 | 4.380 | 10 |
| nps_score | Gateway | 93 | 7.763 | 7.294 | 8.233 | 1 |
| nps_score | Boise | 676 | 7.487 | 7.313 | 7.660 | 2 |
| nps_score | Cotati-Rohnert Park | 497 | 7.268 | 7.079 | 7.457 | 3 |
| nps_score | Dakota | 904 | 6.821 | 6.664 | 6.978 | 4 |
| nps_score | Trenton | 1151 | 6.589 | 6.435 | 6.743 | 5 |
| nps_score | Burlington | 368 | 6.391 | 6.105 | 6.677 | 6 |
| nps_score | ABC | 791 | 6.278 | 6.084 | 6.472 | 7 |
| nps_score | Salinas | 308 | 6.026 | 5.747 | 6.305 | 8 |
| nps_score | Cherry Hill | 499 | 5.959 | 5.716 | 6.201 | 9 |
| nps_score | Rahway | 182 | 4.525 | 4.088 | 4.961 | 10 |
| retention_score | Gateway | 93 | 4.167 | 4.010 | 4.323 | 1 |
| retention_score | Burlington | 359 | 4.006 | 3.910 | 4.101 | 2 |
| retention_score | Trenton | 1129 | 3.961 | 3.904 | 4.019 | 3 |
| retention_score | Boise | 677 | 3.957 | 3.885 | 4.029 | 4 |
| retention_score | Cotati-Rohnert Park | 498 | 3.923 | 3.851 | 3.994 | 5 |
| retention_score | Dakota | 889 | 3.907 | 3.845 | 3.970 | 6 |
| retention_score | Cherry Hill | 500 | 3.745 | 3.656 | 3.834 | 7 |
| retention_score | ABC | 794 | 3.742 | 3.675 | 3.809 | 8 |
| retention_score | Salinas | 314 | 3.567 | 3.471 | 3.664 | 9 |
| retention_score | Rahway | 182 | 3.406 | 3.253 | 3.559 | 10 |
| school_lmc_score | Gateway | 90 | 4.763 | 4.569 | 4.956 | 1 |
| school_lmc_score | Cherry Hill | 490 | 4.751 | 4.677 | 4.825 | 2 |
| school_lmc_score | Cotati-Rohnert Park | 483 | 4.748 | 4.671 | 4.825 | 3 |
| school_lmc_score | Boise | 665 | 4.510 | 4.437 | 4.584 | 4 |
| school_lmc_score | Trenton | 1122 | 4.432 | 4.369 | 4.495 | 5 |
| school_lmc_score | Salinas | 300 | 4.413 | 4.308 | 4.518 | 6 |
| school_lmc_score | Rahway | 180 | 4.402 | 4.262 | 4.542 | 7 |
| school_lmc_score | Burlington | 367 | 4.397 | 4.305 | 4.489 | 8 |
| school_lmc_score | Dakota | 897 | 4.300 | 4.242 | 4.358 | 9 |
| school_lmc_score | ABC | 759 | 4.295 | 4.216 | 4.373 | 10 |
| twb_composite_score | Cotati-Rohnert Park | 523 | 4.782 | 4.731 | 4.834 | 1 |
| twb_composite_score | Gateway | 117 | 4.687 | 4.574 | 4.799 | 2 |
| twb_composite_score | Boise | 723 | 4.621 | 4.571 | 4.670 | 3 |
| twb_composite_score | Burlington | 404 | 4.463 | 4.385 | 4.541 | 4 |
| twb_composite_score | Trenton | 1250 | 4.431 | 4.384 | 4.478 | 5 |
| twb_composite_score | Cherry Hill | 539 | 4.360 | 4.299 | 4.421 | 6 |
| twb_composite_score | ABC | 910 | 4.352 | 4.300 | 4.404 | 7 |
| twb_composite_score | Dakota | 970 | 4.310 | 4.267 | 4.353 | 8 |
| twb_composite_score | Rahway | 208 | 4.201 | 4.095 | 4.307 | 9 |
| twb_composite_score | Salinas | 343 | 4.190 | 4.115 | 4.265 | 10 |
| twb_leader_score_short8 | Cotati-Rohnert Park | 523 | 5.128 | 5.052 | 5.204 | 1 |
| twb_leader_score_short8 | Boise | 723 | 4.860 | 4.777 | 4.944 | 2 |
| twb_leader_score_short8 | Gateway | 117 | 4.741 | 4.585 | 4.897 | 3 |
| twb_leader_score_short8 | Cherry Hill | 536 | 4.549 | 4.454 | 4.645 | 4 |
| twb_leader_score_short8 | Trenton | 1246 | 4.511 | 4.440 | 4.582 | 5 |
| twb_leader_score_short8 | Burlington | 403 | 4.492 | 4.373 | 4.611 | 6 |
| twb_leader_score_short8 | Dakota | 970 | 4.387 | 4.315 | 4.459 | 7 |
| twb_leader_score_short8 | ABC | 812 | 4.337 | 4.245 | 4.429 | 8 |
| twb_leader_score_short8 | Rahway | 208 | 4.290 | 4.130 | 4.450 | 9 |
| twb_leader_score_short8 | Salinas | 343 | 4.271 | 4.142 | 4.400 | 10 |
| union_view_score | Salinas | 308 | 5.073 | 4.965 | 5.180 | 1 |
| union_view_score | Rahway | 182 | 5.007 | 4.884 | 5.130 | 2 |
| union_view_score | Dakota | 904 | 4.953 | 4.900 | 5.006 | 3 |
| union_view_score | Trenton | 1146 | 4.842 | 4.783 | 4.901 | 4 |
| union_view_score | Boise | 675 | 4.837 | 4.758 | 4.916 | 5 |
| union_view_score | Cotati-Rohnert Park | 491 | 4.811 | 4.718 | 4.904 | 6 |
| union_view_score | Cherry Hill | 496 | 4.797 | 4.717 | 4.878 | 7 |
| union_view_score | Burlington | 372 | 4.683 | 4.581 | 4.784 | 8 |
| union_view_score | Gateway | 91 | 4.535 | 4.282 | 4.787 | 9 |
| union_view_score | ABC | 782 | 4.492 | 4.411 | 4.573 | 10 |
# ---------------------------
# 3) Significance flags: which districts differ from GRAND MEAN (BH within variable)
# ---------------------------
district_flags_tbl <- map_dfr(vars_district, ~ district_vs_grandmean(df, .x, min_n = min_n_dist)) %>%
mutate(
variable_label = if ("pretty_names" %in% ls()) pretty_names[variable] else variable
)
# Show only significant districts (above/below) for easy reading
district_sig_tbl <- district_flags_tbl %>%
filter(flag != "Not different") %>%
arrange(variable_label, p_adj, desc(abs(diff)))
knitr::kable(
district_sig_tbl %>%
select(variable_label, district, n, mean, grand_mean, diff, p, p_adj, flag) %>%
mutate(across(c(mean, grand_mean, diff), ~ round(.x, 3)),
across(c(p, p_adj), ~ signif(.x, 3))),
caption = "Districts significantly above/below the overall grand mean (BH-adjusted within variable)"
)| variable_label | district | n | mean | grand_mean | diff | p | p_adj | flag |
|---|---|---|---|---|---|---|---|---|
| district_lmc_score | Cotati-Rohnert Park | 477 | 4.704 | 4.420 | 0.284 | 0.00e+00 | 0.00e+00 | Above grand mean |
| district_lmc_score | Boise | 656 | 4.658 | 4.420 | 0.238 | 0.00e+00 | 0.00e+00 | Above grand mean |
| district_lmc_score | ABC | 749 | 4.248 | 4.420 | -0.171 | 1.00e-06 | 3.40e-06 | Below grand mean |
| district_lmc_score | Salinas | 294 | 4.190 | 4.420 | -0.230 | 2.21e-05 | 5.51e-05 | Below grand mean |
| district_lmc_score | Trenton | 1116 | 4.338 | 4.420 | -0.082 | 5.95e-03 | 1.19e-02 | Below grand mean |
| foundational_support_score | Salinas | 317 | 2.854 | 3.467 | -0.613 | 0.00e+00 | 0.00e+00 | Below grand mean |
| foundational_support_score | Gateway | 115 | 4.400 | 3.467 | 0.933 | 0.00e+00 | 0.00e+00 | Above grand mean |
| foundational_support_score | Cotati-Rohnert Park | 501 | 3.884 | 3.467 | 0.417 | 0.00e+00 | 0.00e+00 | Above grand mean |
| foundational_support_score | Rahway | 184 | 2.973 | 3.467 | -0.494 | 0.00e+00 | 0.00e+00 | Below grand mean |
| foundational_support_score | Cherry Hill | 503 | 3.196 | 3.467 | -0.271 | 0.00e+00 | 0.00e+00 | Below grand mean |
| foundational_support_score | Burlington | 375 | 3.752 | 3.467 | 0.285 | 9.00e-07 | 1.50e-06 | Above grand mean |
| foundational_support_score | Dakota | 910 | 3.328 | 3.467 | -0.139 | 9.00e-06 | 1.28e-05 | Below grand mean |
| foundational_support_score | Trenton | 1162 | 3.573 | 3.467 | 0.106 | 1.66e-03 | 2.08e-03 | Above grand mean |
| foundational_support_score | Boise | 682 | 3.581 | 3.467 | 0.114 | 3.28e-03 | 3.64e-03 | Above grand mean |
| foundational_support_score | ABC | 848 | 3.377 | 3.467 | -0.090 | 1.83e-02 | 1.83e-02 | Below grand mean |
| job_sat | Cotati-Rohnert Park | 499 | 7.553 | 7.051 | 0.502 | 0.00e+00 | 0.00e+00 | Above grand mean |
| job_sat | Rahway | 181 | 6.149 | 7.051 | -0.902 | 0.00e+00 | 0.00e+00 | Below grand mean |
| job_sat | Dakota | 905 | 6.760 | 7.051 | -0.291 | 2.20e-06 | 7.40e-06 | Below grand mean |
| job_sat | Boise | 677 | 7.276 | 7.051 | 0.225 | 1.46e-03 | 3.65e-03 | Above grand mean |
| job_sat | Salinas | 308 | 6.701 | 7.051 | -0.350 | 3.07e-03 | 5.97e-03 | Below grand mean |
| job_sat | Burlington | 371 | 7.345 | 7.051 | 0.294 | 3.58e-03 | 5.97e-03 | Above grand mean |
| job_sat | Cherry Hill | 499 | 6.800 | 7.051 | -0.252 | 5.49e-03 | 7.84e-03 | Below grand mean |
| job_sat | Trenton | 1152 | 7.214 | 7.051 | 0.163 | 9.24e-03 | 1.15e-02 | Above grand mean |
| job_sat | Gateway | 93 | 7.559 | 7.051 | 0.508 | 1.16e-02 | 1.29e-02 | Above grand mean |
| job_stress | Dakota | 906 | 7.591 | 7.128 | 0.462 | 0.00e+00 | 0.00e+00 | Above grand mean |
| job_stress | Burlington | 372 | 6.347 | 7.128 | -0.781 | 0.00e+00 | 0.00e+00 | Below grand mean |
| job_stress | Gateway | 93 | 6.301 | 7.128 | -0.827 | 8.47e-04 | 2.82e-03 | Below grand mean |
| job_stress | Trenton | 1151 | 6.904 | 7.128 | -0.224 | 1.46e-03 | 3.64e-03 | Below grand mean |
| job_stress | Salinas | 308 | 7.461 | 7.128 | 0.333 | 6.90e-03 | 1.38e-02 | Above grand mean |
| lmc_composite_score | ABC | 782 | 4.315 | 4.516 | -0.201 | 0.00e+00 | 0.00e+00 | Below grand mean |
| lmc_composite_score | Cotati-Rohnert Park | 493 | 4.719 | 4.516 | 0.203 | 0.00e+00 | 1.00e-07 | Above grand mean |
| lmc_composite_score | Boise | 675 | 4.647 | 4.516 | 0.131 | 4.73e-05 | 1.58e-04 | Above grand mean |
| nps_score | Boise | 676 | 7.487 | 6.604 | 0.883 | 0.00e+00 | 0.00e+00 | Above grand mean |
| nps_score | Rahway | 182 | 4.525 | 6.604 | -2.079 | 0.00e+00 | 0.00e+00 | Below grand mean |
| nps_score | Cotati-Rohnert Park | 497 | 7.268 | 6.604 | 0.664 | 0.00e+00 | 0.00e+00 | Above grand mean |
| nps_score | Cherry Hill | 499 | 5.959 | 6.604 | -0.645 | 3.00e-07 | 7.00e-07 | Below grand mean |
| nps_score | Gateway | 93 | 7.763 | 6.604 | 1.160 | 5.20e-06 | 1.03e-05 | Above grand mean |
| nps_score | Salinas | 308 | 6.026 | 6.604 | -0.578 | 6.12e-05 | 1.02e-04 | Below grand mean |
| nps_score | ABC | 791 | 6.278 | 6.604 | -0.326 | 1.03e-03 | 1.47e-03 | Below grand mean |
| nps_score | Dakota | 904 | 6.821 | 6.604 | 0.217 | 6.90e-03 | 8.62e-03 | Above grand mean |
| retention_score | Salinas | 314 | 3.567 | 3.862 | -0.294 | 0.00e+00 | 1.00e-07 | Below grand mean |
| retention_score | Rahway | 182 | 3.406 | 3.862 | -0.456 | 0.00e+00 | 1.00e-07 | Below grand mean |
| retention_score | Gateway | 93 | 4.167 | 3.862 | 0.305 | 2.45e-04 | 8.18e-04 | Above grand mean |
| retention_score | ABC | 794 | 3.742 | 3.862 | -0.119 | 4.84e-04 | 1.21e-03 | Below grand mean |
| retention_score | Trenton | 1129 | 3.961 | 3.862 | 0.100 | 6.90e-04 | 1.38e-03 | Above grand mean |
| retention_score | Burlington | 359 | 4.006 | 3.862 | 0.144 | 3.23e-03 | 5.38e-03 | Above grand mean |
| retention_score | Cherry Hill | 500 | 3.745 | 3.862 | -0.116 | 1.08e-02 | 1.35e-02 | Below grand mean |
| retention_score | Boise | 677 | 3.957 | 3.862 | 0.096 | 9.73e-03 | 1.35e-02 | Above grand mean |
| school_lmc_score | Cherry Hill | 490 | 4.751 | 4.459 | 0.292 | 0.00e+00 | 0.00e+00 | Above grand mean |
| school_lmc_score | Cotati-Rohnert Park | 483 | 4.748 | 4.459 | 0.289 | 0.00e+00 | 0.00e+00 | Above grand mean |
| school_lmc_score | Dakota | 897 | 4.300 | 4.459 | -0.159 | 1.00e-07 | 4.00e-07 | Below grand mean |
| school_lmc_score | ABC | 759 | 4.295 | 4.459 | -0.164 | 4.35e-05 | 1.09e-04 | Below grand mean |
| school_lmc_score | Gateway | 90 | 4.763 | 4.459 | 0.304 | 2.75e-03 | 5.50e-03 | Above grand mean |
| twb_composite_score | Cotati-Rohnert Park | 523 | 4.782 | 4.432 | 0.350 | 0.00e+00 | 0.00e+00 | Above grand mean |
| twb_composite_score | Boise | 723 | 4.621 | 4.432 | 0.189 | 0.00e+00 | 0.00e+00 | Above grand mean |
| twb_composite_score | Salinas | 343 | 4.190 | 4.432 | -0.242 | 0.00e+00 | 0.00e+00 | Below grand mean |
| twb_composite_score | Dakota | 970 | 4.310 | 4.432 | -0.122 | 0.00e+00 | 1.00e-07 | Below grand mean |
| twb_composite_score | Gateway | 117 | 4.687 | 4.432 | 0.255 | 2.08e-05 | 4.16e-05 | Above grand mean |
| twb_composite_score | Rahway | 208 | 4.201 | 4.432 | -0.231 | 3.11e-05 | 5.19e-05 | Below grand mean |
| twb_composite_score | ABC | 910 | 4.352 | 4.432 | -0.080 | 2.72e-03 | 3.88e-03 | Below grand mean |
| twb_composite_score | Cherry Hill | 539 | 4.360 | 4.432 | -0.072 | 2.07e-02 | 2.59e-02 | Below grand mean |
| twb_leader_score_short8 | Cotati-Rohnert Park | 523 | 5.128 | 4.549 | 0.579 | 0.00e+00 | 0.00e+00 | Above grand mean |
| twb_leader_score_short8 | Boise | 723 | 4.860 | 4.549 | 0.311 | 0.00e+00 | 0.00e+00 | Above grand mean |
| twb_leader_score_short8 | ABC | 812 | 4.337 | 4.549 | -0.212 | 7.40e-06 | 2.47e-05 | Below grand mean |
| twb_leader_score_short8 | Dakota | 970 | 4.387 | 4.549 | -0.162 | 1.17e-05 | 2.92e-05 | Below grand mean |
| twb_leader_score_short8 | Salinas | 343 | 4.271 | 4.549 | -0.278 | 3.08e-05 | 6.15e-05 | Below grand mean |
| twb_leader_score_short8 | Rahway | 208 | 4.290 | 4.549 | -0.259 | 1.71e-03 | 2.85e-03 | Below grand mean |
| twb_leader_score_short8 | Gateway | 117 | 4.741 | 4.549 | 0.192 | 1.78e-02 | 2.54e-02 | Above grand mean |
| union_view_score | ABC | 782 | 4.492 | 4.805 | -0.314 | 0.00e+00 | 0.00e+00 | Below grand mean |
| union_view_score | Dakota | 904 | 4.953 | 4.805 | 0.148 | 1.00e-07 | 3.00e-07 | Above grand mean |
| union_view_score | Salinas | 308 | 5.073 | 4.805 | 0.267 | 1.80e-06 | 6.00e-06 | Above grand mean |
| union_view_score | Rahway | 182 | 5.007 | 4.805 | 0.202 | 1.54e-03 | 3.84e-03 | Above grand mean |
| union_view_score | Burlington | 372 | 4.683 | 4.805 | -0.122 | 1.87e-02 | 3.74e-02 | Below grand mean |
# ---------------------------
# 4) One plot: district deviations (forest-style) for a chosen variable
# ---------------------------
plot_district_deviation <- function(sig_tbl, var_label_pick) {
d <- sig_tbl %>% filter(variable_label == var_label_pick)
if (nrow(d) == 0) return(NULL)
# order by diff
d <- d %>% arrange(diff) %>%
mutate(district = factor(district, levels = district))
ggplot(d, aes(x = diff, y = district)) +
geom_vline(xintercept = 0, linewidth = 0.7) +
geom_point(size = 2.2) +
labs(
title = paste0("District deviations from grand mean: ", var_label_pick),
subtitle = "Points are district mean minus overall mean (only districts significantly different after BH adjustment).",
x = "Difference from grand mean",
y = NULL
) +
theme_pub(11)
}
# Example usage (uncomment and set the label you want):
# plot_district_deviation(district_sig_tbl, "Teacher Well-Being (Composite)")
# plot_district_deviation(district_sig_tbl, "LMC Composite")“District differences exist across all outcomes, but the largest separations between districts show up in foundational supports, NPS, and overall well-being.”
“Districts differ meaningfully in foundational supports, NPS, and teacher well-being.”
“Some districts are reliably above/below the overall average—these are candidate exemplars and candidate support targets.”
“We’re seeing patterns that align: districts with higher foundational supports tend to show higher TWB and job satisfaction.” thi si snot pairwise comparisons or direct district to district tests, they are using the grand population mean
district_means <- df %>%
group_by(district) %>%
summarise(
n = n(),
lmc = mean(lmc_composite_score, na.rm = TRUE),
twb = mean(twb_composite_score, na.rm = TRUE),
job_sat = mean(job_sat, na.rm = TRUE),
job_stress = mean(job_stress, na.rm = TRUE),
retention = mean(retention_score, na.rm = TRUE),
nps = mean(nps_score, na.rm = TRUE),
.groups = "drop"
)
district_long <- district_means %>%
pivot_longer(cols = c(lmc, twb, job_sat, job_stress, retention, nps),
names_to = "metric", values_to = "mean_value")
metric_labels <- c(
lmc = "LMC Composite",
twb = "TWB Composite",
job_sat = "Job Satisfaction",
job_stress = "Job Stress",
retention = "Retention",
nps = "NPS"
)
district_long <- district_long %>%
mutate(metric_lab = metric_labels[metric])
p_district <- ggplot(district_long, aes(x = reorder(as.character(district), mean_value), y = mean_value)) +
geom_col() +
facet_wrap(~metric_lab, scales = "free_y", ncol = 3) +
coord_flip() +
labs(
title = "District Differences: Collaboration and Outcomes (Means)",
x = "District",
y = "Mean"
) +
theme_pub(10)
p_district
save_pub(p_district, "Fig4_District_Means_Dashboard.png", width = 12, height = 7.5)# Choose demographic predictors that exist in df
demo_vars <- c("role", "grade", "gender", "race", "union_interaction_ord")
# Which outcomes/constructs to test
vars_demo <- c(outcomes, "twb_composite_score", "twb_leader_score_short8", "foundational_support_score",
"lmc_composite_score", "school_lmc_score", "district_lmc_score")
demo_tests <- list()
for (dv in vars_demo) {
for (g in demo_vars) {
if (g %in% names(df)) {
f <- as.formula(paste0(dv, " ~ ", g))
m <- lm(f, data = df)
a <- anova(m)
demo_tests[[paste(dv, g, sep="__")]] <- tibble(
outcome = dv,
group = g,
F = a$`F value`[1],
df1 = a$Df[1],
df2 = a$Df[2],
p = a$`Pr(>F)`[1],
r2 = summary(m)$r.squared
)
}
}
}
demo_tests_tbl <- bind_rows(demo_tests) %>%
mutate(
outcome_label = pretty_names[outcome],
p_adj = p.adjust(p, method = "BH")
) %>%
arrange(p_adj)
demo_tests_tbl %>% mutate(across(c(F,p,p_adj,r2), ~round(.x, 4))) %>% head(30)## # A tibble: 30 × 9
## outcome group F df1 df2 p r2 outcome_label p_adj
## <chr> <chr> <dbl> <int> <int> <dbl> <dbl> <chr> <dbl>
## 1 lmc_composite_score unio… 242. 3 5445 0 0.118 LMC Composite 0
## 2 district_lmc_score unio… 113. 3 5282 0 0.0601 District LMC 0
## 3 job_stress role 79.7 4 5267 0 0.0571 Job Stress 0
## 4 school_lmc_score unio… 88.3 3 5345 0 0.0472 School LMC 0
## 5 nps_score role 37.8 4 5267 0 0.0279 NPS (School+… 0
## 6 nps_score gend… 37.4 4 5271 0 0.0276 NPS (School+… 0
## 7 twb_composite_score gend… 36.8 4 5289 0 0.0271 TWB Composite 0
## 8 twb_leader_score_s… gend… 31.6 4 5281 0 0.0234 TWB Leadersh… 0
## 9 nps_score race 5.83 36 5195 0 0.0388 NPS (School+… 0
## 10 lmc_composite_score race 5.80 36 5200 0 0.0386 LMC Composite 0
## # ℹ 20 more rows# Example: job satisfaction
m_js_mod_district <- lm(job_sat ~ lmc_composite_score * district, data = df)
anova(m_js_mod_district)## Analysis of Variance Table
##
## Response: job_sat
## Df Sum Sq Mean Sq F value Pr(>F)
## lmc_composite_score 1 3413.3 3413.3 1065.780 < 2.2e-16 ***
## district 9 488.5 54.3 16.947 < 2.2e-16 ***
## lmc_composite_score:district 9 203.1 22.6 7.047 3.381e-10 ***
## Residuals 5416 17345.4 3.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Extract just the interaction block p-value (overall test)
# (This is a blunt omnibus check; it answers "any moderation exists?")m_stress_mod_district <- lm(job_stress ~ lmc_composite_score * district, data = df)
m_ret_mod_district <- lm(retention_score ~ lmc_composite_score * district, data = df)
m_nps_mod_district <- lm(nps_score ~ lmc_composite_score * district, data = df)
anova(m_stress_mod_district)## Analysis of Variance Table
##
## Response: job_stress
## Df Sum Sq Mean Sq F value Pr(>F)
## lmc_composite_score 1 284.8 284.763 63.4682 1.971e-15 ***
## district 9 589.5 65.505 14.5998 < 2.2e-16 ***
## lmc_composite_score:district 9 180.7 20.074 4.4741 7.155e-06 ***
## Residuals 5418 24308.9 4.487
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m_ret_mod_district)## Analysis of Variance Table
##
## Response: retention_score
## Df Sum Sq Mean Sq F value Pr(>F)
## lmc_composite_score 1 522.0 521.98 653.4525 < 2.2e-16 ***
## district 9 112.9 12.55 15.7078 < 2.2e-16 ***
## lmc_composite_score:district 9 32.1 3.57 4.4704 7.259e-06 ***
## Residuals 5359 4280.8 0.80
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m_nps_mod_district)## Analysis of Variance Table
##
## Response: nps_score
## Df Sum Sq Mean Sq F value Pr(>F)
## lmc_composite_score 1 7468.5 7468.5 1421.8600 < 2.2e-16 ***
## district 9 1757.7 195.3 37.1819 < 2.2e-16 ***
## lmc_composite_score:district 9 335.8 37.3 7.1033 2.703e-10 ***
## Residuals 5419 28464.1 5.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Role moderation
m_js_mod_role <- lm(job_sat ~ lmc_composite_score * role, data = df)
anova(m_js_mod_role)## Analysis of Variance Table
##
## Response: job_sat
## Df Sum Sq Mean Sq F value Pr(>F)
## lmc_composite_score 1 3441.3 3441.3 1082.0844 <2e-16 ***
## role 4 621.8 155.5 48.8826 <2e-16 ***
## lmc_composite_score:role 4 9.6 2.4 0.7527 0.5561
## Residuals 5252 16702.7 3.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Union interaction moderation (ordered)
m_js_mod_union <- lm(job_sat ~ lmc_composite_score * union_interaction_ord, data = df)
anova(m_js_mod_union)## Analysis of Variance Table
##
## Response: job_sat
## Df Sum Sq Mean Sq F value
## lmc_composite_score 1 3411.7 3411.7 1064.5611
## union_interaction_ord 3 565.1 188.4 58.7790
## lmc_composite_score:union_interaction_ord 3 78.7 26.2 8.1814
## Residuals 5424 17382.9 3.2
## Pr(>F)
## lmc_composite_score < 2.2e-16 ***
## union_interaction_ord < 2.2e-16 ***
## lmc_composite_score:union_interaction_ord 1.972e-05 ***
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1p_role_slope <- ggplot(df, aes(x = lmc_composite_score, y = job_sat)) +
geom_point(alpha = 0.08, size = 1) +
geom_smooth(aes(group = role), method = "lm", se = FALSE) +
facet_wrap(~role) +
labs(
title = "Does LMC Relate to Job Satisfaction Similarly Across Roles?",
x = "LMC Composite",
y = "Job Satisfaction"
) +
theme_pub(11)
p_role_slope
save_pub(p_role_slope, "Fig5_LMC_Slope_by_Role_JobSat.png", width = 12, height = 7)df_plot <- df %>%
select(lmc_composite_score, job_sat, role) %>%
filter(!is.na(lmc_composite_score), !is.na(job_sat), !is.na(role))
p_role_slope <- ggplot(df_plot, aes(x = lmc_composite_score, y = job_sat)) +
geom_point(alpha = 0.06, size = 1) +
geom_smooth(method = "lm", se = TRUE, linewidth = 0.9) +
facet_wrap(~ role, scales = "free_y") +
labs(
title = "LMC → Job Satisfaction by Role",
subtitle = "Each panel shows a role-specific regression line (95% CI).",
x = "LMC Composite",
y = "Job Satisfaction"
) +
theme_pub(11)
p_role_slope
# save_pub(p_role_slope, "Fig5_LMC_Slope_by_Role_JobSat.png", width = 12, height = 7)df_m <- df %>%
select(job_sat, lmc_composite_score, role) %>%
filter(!is.na(job_sat), !is.na(lmc_composite_score), !is.na(role)) %>%
mutate(role = droplevels(as.factor(role)))
m_main <- lm(job_sat ~ lmc_composite_score + role, data = df_m)
m_int <- lm(job_sat ~ lmc_composite_score * role, data = df_m)
anova(m_main, m_int) # omnibus test: do interactions improve fit?## Analysis of Variance Table
##
## Model 1: job_sat ~ lmc_composite_score + role
## Model 2: job_sat ~ lmc_composite_score * role
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 5256 16712
## 2 5252 16703 4 9.5749 0.7527 0.5561summary(m_int) # role-specific slope differences##
## Call:
## lm(formula = job_sat ~ lmc_composite_score * role, data = df_m)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.9642 -0.9529 0.2657 1.2092 6.7744
##
## Coefficients:
## Estimate
## (Intercept) 3.96625
## lmc_composite_score 0.91867
## roleCertificated -1.79405
## roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.) -1.96029
## roleClassified -0.58427
## roleUncertificated (paraprofessionals, clerical, maintenance, etc.) -2.92581
## lmc_composite_score:roleCertificated 0.13470
## lmc_composite_score:roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.) -0.02234
## lmc_composite_score:roleClassified 0.02378
## lmc_composite_score:roleUncertificated (paraprofessionals, clerical, maintenance, etc.) 0.24646
## Std. Error
## (Intercept) 0.93772
## lmc_composite_score 0.21044
## roleCertificated 0.95138
## roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.) 1.27536
## roleClassified 0.98238
## roleUncertificated (paraprofessionals, clerical, maintenance, etc.) 2.48951
## lmc_composite_score:roleCertificated 0.21327
## lmc_composite_score:roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.) 0.28122
## lmc_composite_score:roleClassified 0.22079
## lmc_composite_score:roleUncertificated (paraprofessionals, clerical, maintenance, etc.) 0.53303
## t value
## (Intercept) 4.230
## lmc_composite_score 4.366
## roleCertificated -1.886
## roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.) -1.537
## roleClassified -0.595
## roleUncertificated (paraprofessionals, clerical, maintenance, etc.) -1.175
## lmc_composite_score:roleCertificated 0.632
## lmc_composite_score:roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.) -0.079
## lmc_composite_score:roleClassified 0.108
## lmc_composite_score:roleUncertificated (paraprofessionals, clerical, maintenance, etc.) 0.462
## Pr(>|t|)
## (Intercept) 2.38e-05
## lmc_composite_score 1.29e-05
## roleCertificated 0.0594
## roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.) 0.1243
## roleClassified 0.5520
## roleUncertificated (paraprofessionals, clerical, maintenance, etc.) 0.2399
## lmc_composite_score:roleCertificated 0.5277
## lmc_composite_score:roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.) 0.9367
## lmc_composite_score:roleClassified 0.9142
## lmc_composite_score:roleUncertificated (paraprofessionals, clerical, maintenance, etc.) 0.6438
##
## (Intercept) ***
## lmc_composite_score ***
## roleCertificated .
## roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.)
## roleClassified
## roleUncertificated (paraprofessionals, clerical, maintenance, etc.)
## lmc_composite_score:roleCertificated
## lmc_composite_score:roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.)
## lmc_composite_score:roleClassified
## lmc_composite_score:roleUncertificated (paraprofessionals, clerical, maintenance, etc.)
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.783 on 5252 degrees of freedom
## Multiple R-squared: 0.196, Adjusted R-squared: 0.1947
## F-statistic: 142.3 on 9 and 5252 DF, p-value: < 2.2e-16df_m <- df %>%
select(job_sat, lmc_composite_score, role) %>%
filter(!is.na(job_sat), !is.na(lmc_composite_score), !is.na(role)) %>%
mutate(role = droplevels(as.factor(role)))
m_main <- lm(job_sat ~ lmc_composite_score + role, data = df_m)
m_int <- lm(job_sat ~ lmc_composite_score * role, data = df_m)
anova(m_main, m_int) # omnibus test: do interactions improve fit?## Analysis of Variance Table
##
## Model 1: job_sat ~ lmc_composite_score + role
## Model 2: job_sat ~ lmc_composite_score * role
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 5256 16712
## 2 5252 16703 4 9.5749 0.7527 0.5561summary(m_int) # role-specific slope differences##
## Call:
## lm(formula = job_sat ~ lmc_composite_score * role, data = df_m)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.9642 -0.9529 0.2657 1.2092 6.7744
##
## Coefficients:
## Estimate
## (Intercept) 3.96625
## lmc_composite_score 0.91867
## roleCertificated -1.79405
## roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.) -1.96029
## roleClassified -0.58427
## roleUncertificated (paraprofessionals, clerical, maintenance, etc.) -2.92581
## lmc_composite_score:roleCertificated 0.13470
## lmc_composite_score:roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.) -0.02234
## lmc_composite_score:roleClassified 0.02378
## lmc_composite_score:roleUncertificated (paraprofessionals, clerical, maintenance, etc.) 0.24646
## Std. Error
## (Intercept) 0.93772
## lmc_composite_score 0.21044
## roleCertificated 0.95138
## roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.) 1.27536
## roleClassified 0.98238
## roleUncertificated (paraprofessionals, clerical, maintenance, etc.) 2.48951
## lmc_composite_score:roleCertificated 0.21327
## lmc_composite_score:roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.) 0.28122
## lmc_composite_score:roleClassified 0.22079
## lmc_composite_score:roleUncertificated (paraprofessionals, clerical, maintenance, etc.) 0.53303
## t value
## (Intercept) 4.230
## lmc_composite_score 4.366
## roleCertificated -1.886
## roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.) -1.537
## roleClassified -0.595
## roleUncertificated (paraprofessionals, clerical, maintenance, etc.) -1.175
## lmc_composite_score:roleCertificated 0.632
## lmc_composite_score:roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.) -0.079
## lmc_composite_score:roleClassified 0.108
## lmc_composite_score:roleUncertificated (paraprofessionals, clerical, maintenance, etc.) 0.462
## Pr(>|t|)
## (Intercept) 2.38e-05
## lmc_composite_score 1.29e-05
## roleCertificated 0.0594
## roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.) 0.1243
## roleClassified 0.5520
## roleUncertificated (paraprofessionals, clerical, maintenance, etc.) 0.2399
## lmc_composite_score:roleCertificated 0.5277
## lmc_composite_score:roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.) 0.9367
## lmc_composite_score:roleClassified 0.9142
## lmc_composite_score:roleUncertificated (paraprofessionals, clerical, maintenance, etc.) 0.6438
##
## (Intercept) ***
## lmc_composite_score ***
## roleCertificated .
## roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.)
## roleClassified
## roleUncertificated (paraprofessionals, clerical, maintenance, etc.)
## lmc_composite_score:roleCertificated
## lmc_composite_score:roleCertificated (teachers, nurses, speech therapists, counselors, social workers, etc.)
## lmc_composite_score:roleClassified
## lmc_composite_score:roleUncertificated (paraprofessionals, clerical, maintenance, etc.)
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.783 on 5252 degrees of freedom
## Multiple R-squared: 0.196, Adjusted R-squared: 0.1947
## F-statistic: 142.3 on 9 and 5252 DF, p-value: < 2.2e-16The association between labor–management collaboration and job satisfaction is strong and positive, and does not significantly differ across role categories.
Partner-facing wording
Collaboration matters for everyone. Regardless of role — certificated staff, classified staff, or paraprofessionals — higher labor–management collaboration is associated with higher job satisfaction.
Strategic implication
LMC appears to function as a system-wide lever rather than a role-specific intervention.
library(broom)
outcomes_mod <- c("job_sat","job_stress","twb_composite_score","nps_score","retention_score")
mods <- c("role","race") # add gender, grade, etc. if you want
run_modtest <- function(df, outcome, moderator) {
dat <- df %>%
select(all_of(c(outcome, "lmc_composite_score", moderator))) %>%
filter(!is.na(.data[[outcome]]),
!is.na(lmc_composite_score),
!is.na(.data[[moderator]])) %>%
mutate("{moderator}" := droplevels(as.factor(.data[[moderator]])))
# Guard: if only 1 level, skip
if (nlevels(dat[[moderator]]) < 2) return(NULL)
m_main <- lm(reformulate(c("lmc_composite_score", moderator), response = outcome), data = dat)
m_int <- lm(reformulate(c("lmc_composite_score", moderator, paste0("lmc_composite_score:", moderator)),
response = outcome), data = dat)
a <- anova(m_main, m_int)
tibble(
outcome = outcome,
moderator = moderator,
n = nrow(dat),
p_interaction = a$`Pr(>F)`[2],
delta_r2 = summary(m_int)$r.squared - summary(m_main)$r.squared,
r2_main = summary(m_main)$r.squared,
r2_int = summary(m_int)$r.squared
)
}
mod_results <- purrr::map_dfr(outcomes_mod, \(y) {
purrr::map_dfr(mods, \(m) run_modtest(df, y, m))
}) %>%
mutate(
p_adj = p.adjust(p_interaction, method = "BH")
) %>%
arrange(p_adj)
mod_results## # A tibble: 10 × 8
## outcome moderator n p_interaction delta_r2 r2_main r2_int p_adj
## <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 nps_score role 5264 0.00222 0.00240 0.242 0.245 0.0222
## 2 job_sat race 5222 0.0144 0.00674 0.184 0.191 0.0721
## 3 nps_score race 5224 0.0290 0.00599 0.222 0.228 0.0845
## 4 retention_score race 5170 0.0338 0.00675 0.119 0.125 0.0845
## 5 twb_composite_s… role 5276 0.0791 0.00114 0.283 0.284 0.158
## 6 twb_composite_s… race 5236 0.0957 0.00481 0.277 0.282 0.159
## 7 retention_score role 5207 0.152 0.00113 0.122 0.123 0.218
## 8 job_stress race 5225 0.229 0.00555 0.0365 0.0421 0.286
## 9 job_sat role 5262 0.556 0.000461 0.196 0.196 0.556
## 10 job_stress role 5264 0.539 0.000547 0.0765 0.0770 0.556df_d <- df %>%
select(job_sat, lmc_composite_score, district) %>%
filter(!is.na(job_sat), !is.na(lmc_composite_score), !is.na(district)) %>%
mutate(district = droplevels(as.factor(district)))
m_main <- lm(job_sat ~ lmc_composite_score + district, data = df_d)
m_int <- lm(job_sat ~ lmc_composite_score * district, data = df_d)
anova(m_main, m_int) # does slope differ by district?## Analysis of Variance Table
##
## Model 1: job_sat ~ lmc_composite_score + district
## Model 2: job_sat ~ lmc_composite_score * district
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 5425 17548
## 2 5416 17345 9 203.12 7.047 3.381e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(lme4)
df_d <- df %>%
select(job_sat, lmc_composite_score, district) %>%
filter(!is.na(job_sat), !is.na(lmc_composite_score), !is.na(district)) %>%
mutate(district = droplevels(as.factor(district)))
m0 <- lmer(job_sat ~ lmc_composite_score + (1 | district), data = df_d, REML = FALSE)
m1 <- lmer(job_sat ~ lmc_composite_score + (1 + lmc_composite_score | district), data = df_d, REML = FALSE)
anova(m0, m1) # if significant: slopes vary by district (moderation)## Data: df_d
## Models:
## m0: job_sat ~ lmc_composite_score + (1 | district)
## m1: job_sat ~ lmc_composite_score + (1 + lmc_composite_score | district)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## m0 4 21846 21872 -10919 21838
## m1 6 21814 21853 -10901 21802 36.195 2 1.381e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(emmeans)
df_m <- df %>%
select(job_sat, lmc_composite_score, role) %>%
filter(!is.na(job_sat), !is.na(lmc_composite_score), !is.na(role)) %>%
mutate(role = droplevels(as.factor(role)))
m <- lm(job_sat ~ lmc_composite_score * role, data = df_m)
slopes <- emtrends(m, ~ role, var = "lmc_composite_score") %>%
as.data.frame()
p_slopes <- ggplot(slopes, aes(x = role, y = lmc_composite_score.trend)) +
geom_point(size = 2) +
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.15) +
coord_flip() +
labs(
title = "Estimated LMC → Job Satisfaction Slope by Role",
subtitle = "Points are slopes; bars are 95% CIs. Non-overlap suggests meaningful slope differences.",
x = NULL,
y = "Slope (Δ Job Sat per +1 LMC)"
) +
theme_pub(11)
p_slopes
# Standardize predictors/outcomes for comparable coefficients
df_std <- df %>%
mutate(across(
c(lmc_composite_score, twb_leader_score_short8, foundational_support_score,
job_sat, job_stress, retention_score, nps_score),
~ as.numeric(scale(.x))
))
# Fit core models
m_lead <- lm(twb_leader_score_short8 ~ lmc_composite_score, data = df_std)
m_fsup <- lm(foundational_support_score ~ lmc_composite_score, data = df_std)
m_js <- lm(job_sat ~ lmc_composite_score + twb_leader_score_short8 + foundational_support_score, data = df_std)
m_st <- lm(job_stress ~ lmc_composite_score + twb_leader_score_short8 + foundational_support_score, data = df_std)
m_ret <- lm(retention_score ~ lmc_composite_score + twb_leader_score_short8 + foundational_support_score, data = df_std)
m_nps <- lm(nps_score ~ lmc_composite_score + twb_leader_score_short8 + foundational_support_score, data = df_std)
# Extract standardized betas
get_beta <- function(model, term) {
coef(summary(model))[term, "Estimate"]
}
beta_tbl <- tibble(
from = c("LMC","LMC","LMC","Leadership","Leadership","Leadership","Leadership","Foundational","Foundational","Foundational","Foundational"),
to = c("Leadership","Foundational","Job Sat (direct)","Job Sat","Job Stress","Retention","NPS","Job Sat","Job Stress","Retention","NPS"),
beta = c(
get_beta(m_lead, "lmc_composite_score"),
get_beta(m_fsup, "lmc_composite_score"),
get_beta(m_js, "lmc_composite_score"),
get_beta(m_js, "twb_leader_score_short8"),
get_beta(m_st, "twb_leader_score_short8"),
get_beta(m_ret, "twb_leader_score_short8"),
get_beta(m_nps, "twb_leader_score_short8"),
get_beta(m_js, "foundational_support_score"),
get_beta(m_st, "foundational_support_score"),
get_beta(m_ret, "foundational_support_score"),
get_beta(m_nps, "foundational_support_score")
)
)
beta_tbl## # A tibble: 11 × 3
## from to beta
## <chr> <chr> <dbl>
## 1 LMC Leadership 0.442
## 2 LMC Foundational 0.452
## 3 LMC Job Sat (direct) 0.0871
## 4 Leadership Job Sat 0.247
## 5 Leadership Job Stress -0.00352
## 6 Leadership Retention 0.0921
## 7 Leadership NPS 0.288
## 8 Foundational Job Sat 0.448
## 9 Foundational Job Stress -0.361
## 10 Foundational Retention 0.386
## 11 Foundational NPS 0.394p_path <- ggplot(beta_tbl, aes(x = beta, y = paste(from, "→", to))) +
geom_col() +
labs(
title = "Conceptual Model: Standardized Associations",
subtitle = "Coefficients from standardized regressions (higher magnitude = stronger relationship)",
x = "Standardized Beta",
y = NULL
) +
theme_pub(11)
p_path
save_pub(p_path, "Fig6_PathStrength_ModelMap.png", width = 10, height = 6.5)# ---- Labels ----
outcome_labels <- c(
job_sat = "Job Satisfaction",
job_stress = "Job Stress",
retention_score = "Retention Intentions",
nps_score = "Net Promoter (School+District)"
)
lmc_labels <- c(
union_view_score = "Union View",
school_lmc_score = "School-Level LMC",
district_lmc_score = "District-Level LMC",
lmc_composite_score = "LMC Composite"
)
# ---- Plot theme ----
theme_pub <- function(base_size = 12){
theme_minimal(base_size = base_size) +
theme(
panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold"),
plot.subtitle = element_text(color = "grey30"),
axis.title = element_text(face = "bold"),
strip.text = element_text(face = "bold")
)
}
# ---- Output folder ----
fig_dir <- file.path(data_dir, "figures_systemwide")
if (!dir.exists(fig_dir)) dir.create(fig_dir, recursive = TRUE)
save_pub <- function(p, fname, w=12, h=9){
ggsave(file.path(fig_dir, fname), p, width=w, height=h, dpi=300, bg="white")
}
outcomes <- c("job_sat","job_stress","retention_score","nps_score")
lmc_levels <- c("union_view_score","school_lmc_score","district_lmc_score","lmc_composite_score")
df_fig <- df %>%
select(all_of(outcomes), all_of(lmc_levels)) %>%
pivot_longer(cols = all_of(outcomes), names_to="outcome", values_to="y") %>%
pivot_longer(cols = all_of(lmc_levels), names_to="lmc_level", values_to="lmc") %>%
mutate(
outcome_lab = outcome_labels[outcome],
lmc_lab = lmc_labels[lmc_level]
)
p1 <- ggplot(df_fig, aes(x=lmc, y=y)) +
geom_point(alpha=.10, size=0.9) +
geom_smooth(method="lm", se=TRUE) +
facet_grid(outcome_lab ~ lmc_lab, scales="free_y") +
labs(
title="Systemwide: LMC and Outcomes",
subtitle="Each panel shows the association between a specific LMC level and an outcome (linear fit shown)",
x="LMC score", y=NULL
) +
theme_pub(11)
p1
save_pub(p1, "Fig1_LMC_byLevel_byOutcome.png", 12, 9)pretty_names <- c(
job_sat="Job Satisfaction", job_stress="Job Stress",
retention_score="Retention Intentions", nps_score="NPS (School+District)",
union_view_score="Union View", school_lmc_score="School LMC",
district_lmc_score="District LMC", lmc_composite_score="LMC Composite",
twb_leader_score_short8="TWB Leadership (8-item)", growth_score="TWB Growth",
wellbeing_score="TWB Well-being", acceptance_score="TWB Acceptance",
adaptability_score="TWB Adaptability", depletion_score="TWB Depletion",
foundational_support_score="TWB Foundational Supports",
twb_composite_score="TWB Composite"
)
desc_vars <- c(outcomes, lmc_levels,
"twb_composite_score","twb_leader_score_short8","foundational_support_score",
"growth_score","acceptance_score","adaptability_score","wellbeing_score","depletion_score")
desc_tbl <- df %>%
summarise(across(all_of(desc_vars),
list(
n = ~sum(!is.na(.x)),
mean = ~mean(.x, na.rm=TRUE),
sd = ~sd(.x, na.rm=TRUE)
))) %>%
pivot_longer(everything(),
names_to=c("var","stat"),
names_sep="_(?=[^_]+$)") %>%
pivot_wider(names_from=stat, values_from=value) %>%
mutate(
var_label = pretty_names[var],
mean = round(mean, 3),
sd = round(sd, 3)
) %>%
select(var, var_label, n, mean, sd) %>%
arrange(desc(mean))
knitr::kable(desc_tbl, caption = "Table 1. Systemwide Descriptives (N, Mean, SD)")| var | var_label | n | mean | sd |
|---|---|---|---|---|
| job_stress | Job Stress | 5480 | 7.128 | 2.160 |
| job_sat | Job Satisfaction | 5477 | 7.051 | 1.986 |
| nps_score | NPS (School+District) | 5469 | 6.604 | 2.642 |
| growth_score | TWB Growth | 5745 | 5.195 | 0.733 |
| adaptability_score | TWB Adaptability | 5666 | 5.131 | 0.650 |
| wellbeing_score | TWB Well-being | 5737 | 5.047 | 0.941 |
| union_view_score | Union View | 5447 | 4.805 | 1.012 |
| acceptance_score | TWB Acceptance | 5666 | 4.768 | 0.831 |
| twb_leader_score_short8 | TWB Leadership (8-item) | 5881 | 4.549 | 1.212 |
| lmc_composite_score | LMC Composite | 5454 | 4.516 | 0.831 |
| school_lmc_score | School LMC | 5353 | 4.459 | 0.986 |
| twb_composite_score | TWB Composite | 5987 | 4.432 | 0.762 |
| district_lmc_score | District LMC | 5290 | 4.420 | 0.906 |
| retention_score | Retention Intentions | 5435 | 3.862 | 0.960 |
| depletion_score | TWB Depletion | 5620 | 3.725 | 1.236 |
| foundational_support_score | TWB Foundational Supports | 5597 | 3.467 | 1.088 |
desc_ci <- df %>%
summarise(across(all_of(desc_vars),
list(
n = ~sum(!is.na(.x)),
mean = ~mean(.x, na.rm=TRUE),
sd = ~sd(.x, na.rm=TRUE)
))) %>%
pivot_longer(everything(),
names_to=c("var","stat"),
names_sep="_(?=[^_]+$)") %>%
pivot_wider(names_from=stat, values_from=value) %>%
mutate(
se = sd/sqrt(n),
ci_low = mean - 1.96*se,
ci_high = mean + 1.96*se,
var_label = pretty_names[var]
)
p2 <- ggplot(desc_ci, aes(x=mean, y=reorder(var_label, mean))) +
geom_point(size=2) +
geom_errorbarh(aes(xmin=ci_low, xmax=ci_high), height=.2) +
labs(title="Systemwide Descriptives", subtitle="Means with 95% confidence intervals", x="Mean", y=NULL) +
theme_pub(11)
p2
save_pub(p2, "Fig2_Descriptives_MeansCI.png", 9, 7)# Helper to run a simple model and extract slope + R2
extract_lm <- function(formula_str, data){
m <- lm(as.formula(formula_str), data=data)
s <- summary(m)
co <- broom::tidy(m)
# slope term is always the first predictor after intercept
pred <- co$term[co$term != "(Intercept)"][1]
slope <- co %>% filter(term == pred) %>% select(estimate, std.error, statistic, p.value)
tibble(
model = formula_str,
predictor = pred,
beta = slope$estimate,
se = slope$std.error,
t = slope$statistic,
p = slope$p.value,
r2 = s$r.squared,
n = nobs(m)
)
}
rows <- list()
for (lmc in lmc_levels){
for (y in outcomes){
f <- paste0(y, " ~ ", lmc)
rows[[paste(y,lmc,sep="__")]] <- extract_lm(f, df) %>%
mutate(
outcome = y,
outcome_label = outcome_labels[y],
lmc_label = lmc_labels[lmc]
)
}
}
tbl_lmc_out <- bind_rows(rows) %>%
mutate(
beta = round(beta, 3),
se = round(se, 3),
r2 = round(r2, 3),
p = signif(p, 3)
) %>%
select(outcome_label, lmc_label, beta, se, r2, p, n) %>%
arrange(outcome_label, lmc_label)
knitr::kable(tbl_lmc_out,
caption="Table 2. Systemwide Associations: LMC Level Predicting Outcomes (bivariate OLS)")| outcome_label | lmc_label | beta | se | r2 | p | n |
|---|---|---|---|---|---|---|
| Job Satisfaction | District-Level LMC | 0.860 | 0.028 | 0.154 | 0.000 | 5278 |
| Job Satisfaction | LMC Composite | 0.954 | 0.030 | 0.159 | 0.000 | 5436 |
| Job Satisfaction | School-Level LMC | 0.800 | 0.025 | 0.158 | 0.000 | 5341 |
| Job Satisfaction | Union View | 0.426 | 0.026 | 0.047 | 0.000 | 5430 |
| Job Stress | District-Level LMC | -0.321 | 0.032 | 0.018 | 0.000 | 5280 |
| Job Stress | LMC Composite | -0.275 | 0.035 | 0.011 | 0.000 | 5438 |
| Job Stress | School-Level LMC | -0.221 | 0.030 | 0.010 | 0.000 | 5342 |
| Job Stress | Union View | -0.050 | 0.029 | 0.001 | 0.086 | 5432 |
| Net Promoter (School+District) | District-Level LMC | 1.284 | 0.036 | 0.193 | 0.000 | 5279 |
| Net Promoter (School+District) | LMC Composite | 1.411 | 0.039 | 0.196 | 0.000 | 5439 |
| Net Promoter (School+District) | School-Level LMC | 1.126 | 0.033 | 0.176 | 0.000 | 5342 |
| Net Promoter (School+District) | Union View | 0.681 | 0.034 | 0.068 | 0.000 | 5433 |
| Retention Intentions | District-Level LMC | 0.343 | 0.014 | 0.105 | 0.000 | 5226 |
| Retention Intentions | LMC Composite | 0.376 | 0.015 | 0.105 | 0.000 | 5379 |
| Retention Intentions | School-Level LMC | 0.276 | 0.013 | 0.081 | 0.000 | 5286 |
| Retention Intentions | Union View | 0.195 | 0.013 | 0.042 | 0.000 | 5373 |
This is the clean partner story: “Even controlling for well-being and leadership climate, LMC still matters / does not matter.”
Two versions: - Model A: outcome ~
lmc_composite + twb_composite
- Model B: outcome ~ lmc_composite + TWB subscales
(leadership + foundational + etc.)
library(dplyr)
library(purrr)
library(broom)
library(stringr)
library(tidyr)
# -----------------------------
# Inputs (edit as needed)
# -----------------------------
outcomes <- c("job_sat","job_stress","retention_score","nps_score")
term_labels <- c(
lmc_composite_score = "LMC Composite",
twb_composite_score = "TWB Composite",
twb_leader_score_short8 = "TWB Leadership (8-item)",
foundational_support_score = "TWB Foundational Supports",
growth_score = "TWB Growth Orientation",
wellbeing_score = "TWB Personal Well-being",
acceptance_score = "TWB Acceptance",
adaptability_score = "TWB Adaptability",
depletion_score = "TWB Depletion"
)
outcome_labels <- c(
job_sat = "Job Satisfaction",
job_stress = "Job Stress",
retention_score = "Retention Intentions",
nps_score = "Net Promoter (School+District)"
)
sig_stars <- function(p){
case_when(
is.na(p) ~ "",
p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .10 ~ "†",
TRUE ~ ""
)
}
# -----------------------------
# Core function: fit two models + tidy
# -----------------------------
run_two_models <- function(outcome_var){
# Model A: LMC + TWB composite
fA <- reformulate(c("lmc_composite_score","twb_composite_score"), response = outcome_var)
mA <- lm(fA, data = df)
# Model B: LMC + TWB subscales
predsB <- c("lmc_composite_score",
"twb_leader_score_short8","foundational_support_score",
"growth_score","wellbeing_score","acceptance_score","adaptability_score","depletion_score")
fB <- reformulate(predsB, response = outcome_var)
mB <- lm(fB, data = df)
tidy_keep <- function(m, model_name){
broom::tidy(m) %>%
filter(term %in% c("lmc_composite_score","twb_composite_score",
"twb_leader_score_short8","foundational_support_score",
"growth_score","wellbeing_score","acceptance_score","adaptability_score","depletion_score")) %>%
mutate(
model = model_name,
r2 = summary(m)$r.squared,
n = stats::nobs(m),
outcome = outcome_var
) %>%
select(outcome, model, term, estimate, std.error, statistic, p.value, r2, n)
}
bind_rows(
tidy_keep(mA, "A: LMC + TWB composite"),
tidy_keep(mB, "B: LMC + TWB subscales")
)
}
# -----------------------------
# Build Table 3
# -----------------------------
tbl3 <- map_dfr(outcomes, run_two_models) %>%
mutate(
outcome_label = recode(outcome, !!!outcome_labels),
term_label = recode(term, !!!term_labels),
stars = sig_stars(p.value),
est_se = sprintf("%.3f (%.3f)%s", estimate, std.error, stars),
p_value_fmt = ifelse(is.na(p.value), NA_character_,
ifelse(p.value < .001, "<.001", sprintf("%.3f", p.value))),
r2 = round(r2, 3)
) %>%
select(outcome_label, model, term_label, est_se, statistic, p_value_fmt, r2, n) %>%
arrange(outcome_label, model, term_label)
tbl3## # A tibble: 40 × 8
## outcome_label model term_label est_se statistic p_value_fmt r2 n
## <chr> <chr> <chr> <chr> <dbl> <chr> <dbl> <int>
## 1 Job Satisfaction A: LMC … LMC Compo… 0.278… 9.54 <.001 0.395 5436
## 2 Job Satisfaction A: LMC … TWB Compo… 1.511… 46.0 <.001 0.395 5436
## 3 Job Satisfaction B: LMC … LMC Compo… 0.191… 6.95 <.001 0.482 5391
## 4 Job Satisfaction B: LMC … TWB Accep… 0.085… 2.65 0.008 0.482 5391
## 5 Job Satisfaction B: LMC … TWB Adapt… -0.00… -0.0201 0.984 0.482 5391
## 6 Job Satisfaction B: LMC … TWB Deple… -0.25… -14.1 <.001 0.482 5391
## 7 Job Satisfaction B: LMC … TWB Found… 0.680… 29.2 <.001 0.482 5391
## 8 Job Satisfaction B: LMC … TWB Growt… 0.248… 7.86 <.001 0.482 5391
## 9 Job Satisfaction B: LMC … TWB Leade… 0.380… 19.0 <.001 0.482 5391
## 10 Job Satisfaction B: LMC … TWB Perso… 0.035… 1.47 0.141 0.482 5391
## # ℹ 30 more rowsWhat this table is doing (one sentence)
You are testing whether Labor–Management Collaboration (LMC) predicts key outcomes on its own (Model A) and whether it still matters after accounting for specific dimensions of Teacher Well-Being (TWB) (Model B).
- Job Satisfaction Model A: LMC + TWB Composite
LMC Composite: 0.278*
TWB Composite: 1.511*
R² = 0.395
Interpretation: Both collaboration and overall well-being strongly predict job satisfaction. TWB is the dominant predictor, but LMC independently contributes. Together they explain ~40% of variance.
Model B: LMC + TWB Subscales
LMC Composite: 0.191* (still strong)
Largest TWB drivers:
Foundational Supports: 0.680* (largest effect)
Leadership (8-item): 0.380*
Growth Orientation: 0.248*
Depletion: –0.253* (strong negative)
Acceptance: small positive
Adaptability & Personal Well-being: not unique predictors
R² = 0.482
Interpretation: Job satisfaction is primarily driven by systems and supports (foundational resources, leadership quality, low depletion). LMC remains a significant, independent contributor, even after accounting for all TWB components.
Plain-language takeaway:
Collaboration boosts job satisfaction above and beyond individual well-being—especially when foundational supports and leadership are strong.
- Job Stress Model A: LMC + TWB Composite
LMC Composite: not significant
TWB Composite: –0.706*
R² = 0.055
Interpretation: Stress is weakly explained at this level; overall well-being matters, but collaboration alone does not explain much stress variance.
Model B: LMC + TWB Subscales
LMC Composite: 0.110 (positive, significant)
Primary drivers:
Depletion: 0.715* (dominant)
Foundational Supports: –0.445*
Acceptance: –0.179*
Growth Orientation: 0.171*
Leadership & Personal Well-being: not significant
R² = 0.286
Interpretation (important nuance):
Stress is overwhelmingly a depletion + supports story
Once depletion and supports are accounted for, LMC’s coefficient turns positive due to shared variance / suppression
This does not mean collaboration causes stress
Correct interpretation:
Collaboration operates upstream through supports and depletion; stress reduction depends on whether collaboration translates into real resources and workload relief.
- Net Promoter Score (NPS) Model A: LMC + TWB Composite
LMC Composite: 0.532*
TWB Composite: 1.963*
R² = 0.422
Interpretation: Both collaboration and well-being strongly predict advocacy and recommendation behavior.
Model B: LMC + TWB Subscales
LMC Composite: 0.434* (very strong)
Largest predictors:
Foundational Supports: 0.875* (largest)
Leadership (8-item): 0.594*
Growth Orientation: 0.357*
Depletion: –0.192*
Personal Well-being: –0.157* (unique negative effect)
R² = 0.475
Interpretation: NPS is the clearest systems signal:
People recommend districts where collaboration, leadership, and foundational supports are strong.
LMC remains one of the strongest independent predictors.
Partner framing:
Collaboration isn’t just internal—it shapes how willing educators are to publicly stand behind their district.
- Retention Intentions Model A: LMC + TWB Composite
LMC Composite: 0.163*
TWB Composite: 0.472*
R² = 0.203
Interpretation: Both collaboration and well-being predict intent to stay, with TWB stronger overall.
Model B: LMC + TWB Subscales
LMC Composite: 0.114* (still meaningful)
Strong predictors:
Foundational Supports: 0.312*
Growth Orientation: 0.154*
Leadership (8-item): 0.059*
Depletion: –0.043*
Adaptability: –0.070 (small but meaningful)
R² = 0.266
Interpretation: Retention is shaped by organizational conditions, not just personal resilience. LMC independently predicts retention even after accounting for leadership and supports.
Cross-Outcome Synthesis (the real insight) 1. LMC is not redundant with TWB
Across all four outcomes, LMC remains significant in the full models.
- Foundational supports are the strongest universal lever
They dominate satisfaction, stress, retention, and NPS.
- Depletion is the main driver of stress—and a drag everywhere else
Reducing burnout matters more than boosting generic “well-being.”
- Leadership is especially critical for satisfaction and advocacy
Leadership shows up strongest for job satisfaction and NPS.
One-sentence executive takeaway
Labor–Management Collaboration consistently predicts better satisfaction, retention, and advocacy—and its benefits are strongest when collaboration translates into real foundational supports, effective leadership, and reduced depletion.
Systems most strongly shape whether educators feel supported, valued, and would recommend or remain — less so their personal growth orientation. Demographics shape experience — but systems explain more variance than identity alone. This reframes equity conversations:
Not “Which groups are struggling?”
But “Which system conditions create differential exposure to strain or support?”
Labor–Management Collaboration operates as a structural resource in school systems. When collaboration is stronger, educators report:
better leadership experiences,
stronger foundational supports,
higher well-being,
lower stress,
and greater willingness to stay and recommend their school or district.
These relationships persist across roles, genders, grades, and racial groups, suggesting LMC is not merely a perception issue — it is a system condition that shapes experience. Where labor and management collaborate effectively, educators experience stronger leadership, greater well-being, and higher retention — regardless of who they are or where they sit in the system.
# ============================================================
# CORE LMC FIGURES (publication-ready ggplots)
# - Robust: handles missing, trims extremes if desired
# - Works across outcomes + LMC levels
# ============================================================
library(tidyverse)
library(scales)
# ---------------------------
# 0) Plot theme (clean, journal-ish)
# ---------------------------
theme_pub <- function(base_size = 12) {
theme_minimal(base_size = base_size) +
theme(
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
axis.title = element_text(face = "bold"),
plot.title = element_text(face = "bold", size = base_size + 2),
plot.subtitle = element_text(size = base_size),
legend.title = element_text(face = "bold"),
legend.position = "top"
)
}
# ------------------------------------------------------------
# Update helper functions so NA LMC / NA outcomes never appear
# in plots (prevents "NA" quartile bin and removes ggplot warnings)
# ------------------------------------------------------------
make_long_outcomes <- function(df, outcomes) {
df %>%
select(pid, district, school, role, gender, race, grade,
union_interaction_ord,
lmc_composite_score, school_lmc_score, district_lmc_score, union_view_score,
all_of(outcomes)) %>%
pivot_longer(cols = all_of(outcomes),
names_to = "outcome",
values_to = "value") %>%
# Drop missing outcome values early
filter(!is.na(value))
}
add_lmc_bins <- function(df_long, lmc_var = "lmc_composite_score", n_bins = 4) {
df_long %>%
mutate(lmc_value = .data[[lmc_var]]) %>%
# Drop missing LMC values early (prevents NA bin)
filter(!is.na(lmc_value)) %>%
mutate(
lmc_bin = ntile(lmc_value, n_bins),
lmc_bin = factor(lmc_bin, levels = 1:n_bins,
labels = paste0("Q", 1:n_bins))
)
}
# ---------------------------
# 3) Helper: summary stats for binned plots
# ---------------------------
summarize_bins <- function(df_binned) {
df_binned %>%
group_by(outcome, lmc_bin) %>%
summarise(
n = sum(!is.na(value) & !is.na(lmc_value)),
mean = mean(value, na.rm = TRUE),
se = sd(value, na.rm = TRUE) / sqrt(sum(!is.na(value))),
ci_low = mean - 1.96 * se,
ci_high = mean + 1.96 * se,
.groups = "drop"
)
}
# ---------------------------
# 4) Figure 1: Dose–response (LMC quartiles -> outcome means)
# ---------------------------
plot_lmc_quartiles <- function(df, outcomes,
lmc_var = "lmc_composite_score",
lmc_label = "LMC Composite",
outcome_labels = NULL) {
df_long <- make_long_outcomes(df, outcomes) %>%
add_lmc_bins(lmc_var = lmc_var, n_bins = 4)
summ <- summarize_bins(df_long)
if (!is.null(outcome_labels)) {
summ <- summ %>%
mutate(outcome_lab = recode(outcome, !!!outcome_labels))
} else {
summ <- summ %>% mutate(outcome_lab = outcome)
}
ggplot(summ, aes(x = lmc_bin, y = mean, group = 1)) +
geom_line(linewidth = 0.8) +
geom_point(size = 2.2) +
geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.12) +
facet_wrap(~ outcome_lab, scales = "free_y", ncol = 2) +
labs(
title = paste0(lmc_label, " Quartiles and Outcomes"),
subtitle = "Points = mean; error bars = 95% CI. Quartiles computed within full sample.",
x = paste0(lmc_label, " (Quartiles)"),
y = "Outcome"
) +
theme_pub()
}
# ---------------------------
# 5) Figure 2: Scatter + smooth (continuous LMC -> outcome)
# ---------------------------
plot_lmc_scatter <- function(df, outcomes,
lmc_var = "lmc_composite_score",
lmc_label = "LMC Composite",
outcome_labels = NULL,
alpha_points = 0.08) {
df_long <- make_long_outcomes(df, outcomes) %>%
mutate(lmc_value = .data[[lmc_var]])
if (!is.null(outcome_labels)) {
df_long <- df_long %>%
mutate(outcome_lab = recode(outcome, !!!outcome_labels))
} else {
df_long <- df_long %>% mutate(outcome_lab = outcome)
}
ggplot(df_long, aes(x = lmc_value, y = value)) +
geom_point(alpha = alpha_points) +
geom_smooth(method = "lm", se = TRUE, linewidth = 0.9) +
facet_wrap(~ outcome_lab, scales = "free_y", ncol = 2) +
labs(
title = paste0(lmc_label, " and Outcomes (Continuous)"),
subtitle = "Linear fit with 95% CI. Points shown with transparency for large-N readability.",
x = lmc_label,
y = "Outcome"
) +
theme_pub()
}
# ---------------------------
# 6) Figure 3: District-level variation (district means)
# ---------------------------
plot_district_means <- function(df, outcomes,
lmc_var = "lmc_composite_score",
lmc_label = "LMC Composite",
outcome_labels = NULL,
min_n = 30) {
# district summaries
dist_summ <- df %>%
group_by(district) %>%
summarise(
n = n(),
lmc_mean = mean(.data[[lmc_var]], na.rm = TRUE),
across(all_of(outcomes), ~ mean(.x, na.rm = TRUE), .names = "mean_{.col}"),
.groups = "drop"
) %>%
filter(n >= min_n)
dist_long <- dist_summ %>%
pivot_longer(cols = starts_with("mean_"),
names_to = "outcome",
values_to = "outcome_mean") %>%
mutate(outcome = str_remove(outcome, "^mean_"))
if (!is.null(outcome_labels)) {
dist_long <- dist_long %>%
mutate(outcome_lab = recode(outcome, !!!outcome_labels))
} else {
dist_long <- dist_long %>% mutate(outcome_lab = outcome)
}
ggplot(dist_long, aes(x = lmc_mean, y = outcome_mean)) +
geom_point(size = 2.2, alpha = 0.85) +
geom_smooth(method = "lm", se = TRUE, linewidth = 0.9) +
facet_wrap(~ outcome_lab, scales = "free_y", ncol = 2) +
labs(
title = "District-Level Means: LMC and Outcomes",
subtitle = paste0("Each point is a district mean (districts with n ≥ ", min_n, ")."),
x = paste0(lmc_label, " (District Mean)"),
y = "Outcome (District Mean)"
) +
theme_pub()
}
# ---------------------------
# 7) Figure 4: Stratified dose–response (by role or gender)
# ---------------------------
plot_lmc_quartiles_by_group <- function(df, outcomes,
group_var = "role",
lmc_var = "lmc_composite_score",
lmc_label = "LMC Composite",
outcome_labels = NULL) {
df_long <- make_long_outcomes(df, outcomes) %>%
add_lmc_bins(lmc_var = lmc_var, n_bins = 4) %>%
mutate(group = .data[[group_var]])
summ <- df_long %>%
group_by(outcome, group, lmc_bin) %>%
summarise(
n = sum(!is.na(value) & !is.na(lmc_value)),
mean = mean(value, na.rm = TRUE),
se = sd(value, na.rm = TRUE) / sqrt(sum(!is.na(value))),
ci_low = mean - 1.96 * se,
ci_high = mean + 1.96 * se,
.groups = "drop"
)
if (!is.null(outcome_labels)) {
summ <- summ %>% mutate(outcome_lab = recode(outcome, !!!outcome_labels))
} else {
summ <- summ %>% mutate(outcome_lab = outcome)
}
ggplot(summ, aes(x = lmc_bin, y = mean, group = group, color = group)) +
geom_line(linewidth = 0.8) +
geom_point(size = 2) +
facet_wrap(~ outcome_lab, scales = "free_y", ncol = 2) +
labs(
title = paste0(lmc_label, " Quartiles and Outcomes by ", str_to_title(group_var)),
subtitle = "Group-specific means across LMC quartiles (95% CI not shown for readability).",
x = paste0(lmc_label, " (Quartiles)"),
y = "Outcome",
color = str_to_title(group_var)
) +
theme_pub()
}
# ============================================================
# USAGE (drop into your Rmd after scoring)
# ============================================================
# 1) Choose outcomes to display together
core_outcomes <- c("twb_composite_score", "job_sat", "job_stress", "retention_score", "nps_score")
# Pretty labels
outcome_labels <- c(
twb_composite_score = "Teacher Well-Being (Composite)",
job_sat = "Job Satisfaction",
job_stress = "Job Stress",
retention_score = "Retention Intentions",
nps_score = "Net Promoter Score (School+District)"
)
# ---- Figure set A: LMC Composite ----
p1 <- plot_lmc_quartiles(df, core_outcomes,
lmc_var = "lmc_composite_score",
lmc_label = "LMC Composite",
outcome_labels = outcome_labels)
p2 <- plot_lmc_scatter(df, core_outcomes,
lmc_var = "lmc_composite_score",
lmc_label = "LMC Composite",
outcome_labels = outcome_labels)
p3 <- plot_district_means(df, core_outcomes,
lmc_var = "lmc_composite_score",
lmc_label = "LMC Composite",
outcome_labels = outcome_labels,
min_n = 30)
p4_role <- plot_lmc_quartiles_by_group(df, core_outcomes,
group_var = "role",
lmc_var = "lmc_composite_score",
lmc_label = "LMC Composite",
outcome_labels = outcome_labels)
# Print to knit output
p1
p2
p3
p4_role
# ---- Figure set B: School LMC vs District LMC (repeat) ----
p_school <- plot_lmc_quartiles(df, core_outcomes,
lmc_var = "school_lmc_score",
lmc_label = "School-Level LMC",
outcome_labels = outcome_labels)
p_district <- plot_lmc_quartiles(df, core_outcomes,
lmc_var = "district_lmc_score",
lmc_label = "District-Level LMC",
outcome_labels = outcome_labels)
p_school
p_district
# ============================================================
# OPTIONAL: Save publication-ready figures
# ============================================================
# ggsave("Fig_LMC_Quartiles_CoreOutcomes.png", p1, width = 9, height = 7, dpi = 300)
# ggsave("Fig_LMC_Scatter_CoreOutcomes.png", p2, width = 9, height = 7, dpi = 300)
# ggsave("Fig_LMC_DistrictMeans_CoreOutcomes.png", p3, width = 9, height = 7, dpi = 300)
# ggsave("Fig_LMC_Quartiles_ByRole.png", p4_role, width = 9, height = 7, dpi = 300)“Across the sample, educators reporting stronger labor–management collaboration also report meaningfully higher satisfaction, well-being, retention, and recommendation — and somewhat lower stress.”
# ============================================================
# Single combined figure:
# Composite vs School-level vs District-level LMC (quartiles)
# Facets = outcomes; Columns = LMC level; Points=mean, bars=95% CI
# ============================================================
plot_lmc_quartiles_threepanel <- function(df, outcomes, outcome_labels = NULL) {
# Build one long dataset repeated for each LMC "level"
df_long_base <- make_long_outcomes(df, outcomes)
lmc_specs <- tibble::tibble(
lmc_var = c("lmc_composite_score", "school_lmc_score", "district_lmc_score"),
lmc_panel = c("LMC Composite", "School-Level LMC", "District-Level LMC")
)
# Create binned summaries for each LMC spec
summ_all <- purrr::map2_dfr(
lmc_specs$lmc_var, lmc_specs$lmc_panel,
~{
df_long_base %>%
add_lmc_bins(lmc_var = .x, n_bins = 4) %>%
summarize_bins() %>%
mutate(lmc_panel = .y)
}
)
# Apply pretty outcome labels
if (!is.null(outcome_labels)) {
summ_all <- summ_all %>%
mutate(outcome_lab = dplyr::recode(outcome, !!!outcome_labels))
} else {
summ_all <- summ_all %>% mutate(outcome_lab = outcome)
}
ggplot(summ_all, aes(x = lmc_bin, y = mean, group = 1)) +
geom_line(linewidth = 0.8) +
geom_point(size = 2.2) +
geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.12) +
facet_grid(outcome_lab ~ lmc_panel, scales = "free_y") +
labs(
title = "LMC Quartiles and Outcomes Across Levels",
subtitle = "Points = mean; error bars = 95% CI. Quartiles computed within full sample; missing values excluded.",
x = "LMC Quartile",
y = "Outcome"
) +
theme_pub()
}
# ---- Build and print the combined figure ----
p_lmc_3panel <- plot_lmc_quartiles_threepanel(
df,
outcomes = core_outcomes,
outcome_labels = outcome_labels
)
p_lmc_3panel
# Optional save
# ggsave("Fig_LMC_ThreePanel_Quartiles.png", p_lmc_3panel, width = 12, height = 9, dpi = 300)library(tidyverse)
theme_pub <- function(base_size = 12) {
theme_minimal(base_size = base_size) +
theme(
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
axis.title = element_text(face = "bold"),
plot.title = element_text(face = "bold", size = base_size + 2),
plot.subtitle = element_text(size = base_size),
legend.position = "top"
)
}
outcomes <- c("job_sat","job_stress","retention_score","nps_score","twb_composite_score")
outcome_labels <- c(
job_sat = "Job Satisfaction",
job_stress = "Job Stress",
retention_score = "Retention Intentions",
nps_score = "Net Promoter (School+District)",
twb_composite_score = "Teacher Well-Being (Composite)"
)
df_plot <- df %>%
select(lmc_composite_score, all_of(outcomes)) %>%
mutate(across(everything(), ~ suppressWarnings(as.numeric(as.character(.))))) %>%
pivot_longer(cols = all_of(outcomes), names_to = "outcome", values_to = "y") %>%
filter(!is.na(lmc_composite_score), !is.na(y)) %>%
mutate(outcome_lab = recode(outcome, !!!outcome_labels))
ggplot(df_plot, aes(x = lmc_composite_score, y = y)) +
geom_point(alpha = 0.06) +
geom_smooth(method = "lm", se = TRUE, linewidth = 0.9) +
facet_wrap(~ outcome_lab, scales = "free_y", ncol = 2) +
labs(
title = "Labor–Management Collaboration and Key Outcomes",
subtitle = "Each panel shows a linear association with 95% CI; rows with missing values for that panel are removed.",
x = "Labor–Management Collaboration (Composite)",
y = NULL
) +
theme_pub()
df_long <- df %>%
select(job_sat, union_view_score, school_lmc_score, district_lmc_score) %>%
mutate(across(everything(), ~ suppressWarnings(as.numeric(as.character(.))))) %>%
pivot_longer(
cols = c(union_view_score, school_lmc_score, district_lmc_score),
names_to = "lmc_level",
values_to = "lmc_score"
) %>%
filter(!is.na(job_sat), !is.na(lmc_score)) %>%
mutate(
lmc_level = recode(
lmc_level,
union_view_score = "Union View",
school_lmc_score = "School-Level LMC",
district_lmc_score = "District-Level LMC"
)
)
ggplot(df_long, aes(x = lmc_score, y = job_sat)) +
geom_point(alpha = 0.07) +
geom_smooth(method = "lm", se = TRUE, linewidth = 0.9) +
facet_wrap(~ lmc_level, scales = "free_x", ncol = 3) +
labs(
title = "Job Satisfaction by Level of Collaboration",
subtitle = "Each panel uses that LMC level’s original scale (free x-axis). Linear fit with 95% CI.",
x = "LMC Score",
y = "Job Satisfaction"
) +
theme_pub()
z_score <- function(x) {
m <- mean(x, na.rm = TRUE)
s <- sd(x, na.rm = TRUE)
if (is.na(s) || s == 0) return(rep(NA_real_, length(x)))
(x - m) / s
}
df_long_z <- df %>%
select(job_sat, union_view_score, school_lmc_score, district_lmc_score) %>%
mutate(across(everything(), ~ suppressWarnings(as.numeric(as.character(.))))) %>%
mutate(
union_view_z = z_score(union_view_score),
school_lmc_z = z_score(school_lmc_score),
district_lmc_z = z_score(district_lmc_score)
) %>%
select(job_sat, union_view_z, school_lmc_z, district_lmc_z) %>%
pivot_longer(
cols = c(union_view_z, school_lmc_z, district_lmc_z),
names_to = "lmc_level",
values_to = "lmc_z"
) %>%
filter(!is.na(job_sat), !is.na(lmc_z)) %>%
mutate(
lmc_level = recode(
lmc_level,
union_view_z = "Union View (z)",
school_lmc_z = "School-Level LMC (z)",
district_lmc_z = "District-Level LMC (z)"
)
)
ggplot(df_long_z, aes(x = lmc_z, y = job_sat)) +
geom_point(alpha = 0.07) +
geom_smooth(method = "lm", se = TRUE, linewidth = 0.9) +
facet_wrap(~ lmc_level, ncol = 3) +
labs(
title = "Job Satisfaction by Collaboration Level (Standardized)",
subtitle = "LMC predictors are z-scored; slopes are directly comparable across panels. Linear fit with 95% CI.",
x = "LMC (z-score)",
y = "Job Satisfaction"
) +
theme_pub()
library(tidyverse)
# ---------------------------
# 0) Publication-ish theme
# ---------------------------
theme_pub <- function(base_size = 12) {
theme_minimal(base_size = base_size) +
theme(
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
axis.title = element_text(face = "bold"),
plot.title = element_text(face = "bold", size = base_size + 2),
plot.subtitle = element_text(size = base_size),
strip.text = element_text(face = "bold"),
legend.position = "top"
)
}
# ---------------------------
# 1) Helpers
# ---------------------------
zscore <- function(x) {
(x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE)
}
# Labels (edit if you want different phrasing)
outcome_labels <- c(
job_stress = "Job Stress",
retention_score = "Retention Intentions",
nps_score = "Net Promoter (School+District)",
twb_composite_score = "Teacher Well-Being (Composite)"
)
lmc_labels <- c(
lmc_composite_score = "LMC Composite",
district_lmc_score = "District-Level LMC",
school_lmc_score = "School-Level LMC",
union_view_score = "Union View"
)
# ---------------------------
# 2) Choose outcomes + LMC levels (4 levels)
# ---------------------------
outcomes <- c("job_stress", "retention_score", "nps_score", "twb_composite_score")
lmc_vars <- c("lmc_composite_score", "district_lmc_score", "school_lmc_score", "union_view_score")
# ---------------------------
# 3) Long dataframe (drops missing pairs automatically)
# - We standardize x and y within their variable for comparable slopes visually
# ---------------------------
df_long <- df %>%
select(pid, all_of(outcomes), all_of(lmc_vars)) %>%
pivot_longer(
cols = all_of(outcomes),
names_to = "outcome",
values_to = "y"
) %>%
pivot_longer(
cols = all_of(lmc_vars),
names_to = "lmc_level",
values_to = "x"
) %>%
# Drop missing x or y for plotting (this is the "remove NAs" you wanted)
filter(!is.na(x) & !is.na(y)) %>%
group_by(outcome) %>% mutate(y_z = zscore(y)) %>% ungroup() %>%
group_by(lmc_level) %>% mutate(x_z = zscore(x)) %>% ungroup() %>%
mutate(
outcome_lab = recode(outcome, !!!outcome_labels),
lmc_lab = recode(lmc_level, !!!lmc_labels)
)
# Optional sanity check (counts)
df_long %>%
count(outcome_lab, lmc_lab) %>%
arrange(outcome_lab, lmc_lab)## # A tibble: 16 × 3
## outcome_lab lmc_lab n
## <chr> <chr> <int>
## 1 Job Stress District-Level LMC 5280
## 2 Job Stress LMC Composite 5438
## 3 Job Stress School-Level LMC 5342
## 4 Job Stress Union View 5432
## 5 Net Promoter (School+District) District-Level LMC 5279
## 6 Net Promoter (School+District) LMC Composite 5439
## 7 Net Promoter (School+District) School-Level LMC 5342
## 8 Net Promoter (School+District) Union View 5433
## 9 Retention Intentions District-Level LMC 5226
## 10 Retention Intentions LMC Composite 5379
## 11 Retention Intentions School-Level LMC 5286
## 12 Retention Intentions Union View 5373
## 13 Teacher Well-Being (Composite) District-Level LMC 5289
## 14 Teacher Well-Being (Composite) LMC Composite 5453
## 15 Teacher Well-Being (Composite) School-Level LMC 5352
## 16 Teacher Well-Being (Composite) Union View 5446# ---------------------------
# 4) Figure A: Raw scale scatter + linear fit (one outcome x 4 panels)
# ---------------------------
p_raw <- ggplot(df_long, aes(x = x, y = y)) +
geom_point(alpha = 0.08) +
geom_smooth(method = "lm", se = TRUE, linewidth = 0.9) +
facet_grid(outcome_lab ~ lmc_lab, scales = "free") +
labs(
title = "Outcomes vs. LMC by Level (Raw Scales)",
subtitle = "Each panel shows a linear fit (95% CI). Points are transparent for large-N readability.",
x = "LMC score",
y = "Outcome"
) +
theme_pub()
p_raw
# ---------------------------
# 5) Figure B: Standardized scatter + linear fit (comparable across panels)
# - y-axis is in SD units; x-axis is in SD units
# - This makes slopes visually comparable across LMC levels/outcomes
# ---------------------------
p_std <- ggplot(df_long, aes(x = x_z, y = y_z)) +
geom_point(alpha = 0.08) +
geom_smooth(method = "lm", se = TRUE, linewidth = 0.9) +
facet_grid(outcome_lab ~ lmc_lab) +
labs(
title = "Outcomes vs. LMC by Level (Standardized)",
subtitle = "Both axes standardized within variable (SD units) so slopes are comparable across panels.",
x = "LMC (z)",
y = "Outcome (z)"
) +
theme_pub()
p_std
# ---------------------------
# 6) OPTIONAL: save
# ---------------------------
# ggsave("Fig_LMC_levels_outcomes_raw.png", p_raw, width = 14, height = 10, dpi = 300)
# ggsave("Fig_LMC_levels_outcomes_std.png", p_std, width = 14, height = 10, dpi = 300)library(tidyverse)
library(broom)
# ----------------------------
# Outcomes + predictors
# ----------------------------
outcomes <- c(
"job_sat",
"job_stress",
"retention_score",
"nps_score"
)
predictors <- c(
"lmc_composite_score",
"twb_leader_score_short8",
"foundational_support_score",
"growth_score",
"wellbeing_score",
"acceptance_score",
"adaptability_score",
"depletion_score"
)
pretty_terms <- c(
lmc_composite_score = "LMC Composite",
twb_leader_score_short8 = "TWB Leadership",
foundational_support_score = "TWB Foundational Supports",
growth_score = "TWB Growth Orientation",
wellbeing_score = "TWB Personal Well-being",
acceptance_score = "TWB Acceptance",
adaptability_score = "TWB Adaptability",
depletion_score = "TWB Depletion"
)
pretty_outcomes <- c(
job_sat = "Job Satisfaction",
job_stress = "Job Stress",
retention_score = "Retention Intentions",
nps_score = "Net Promoter Score"
)
# ----------------------------
# Helper: standardize variables
# ----------------------------
z <- function(x) as.numeric(scale(x))
df_z <- df %>%
select(all_of(outcomes), all_of(predictors)) %>%
drop_na() %>%
mutate(across(everything(), z))
# ----------------------------
# Fit standardized models
# ----------------------------
std_betas <- map_dfr(outcomes, function(out) {
f <- as.formula(
paste(out, "~", paste(predictors, collapse = " + "))
)
lm(f, data = df_z) %>%
tidy(conf.int = TRUE) %>%
filter(term != "(Intercept)") %>%
mutate(
outcome = out,
term_label = pretty_terms[term],
outcome_label = pretty_outcomes[out]
)
})
# ----------------------------
# Plot: standardized betas
# ----------------------------
p_std <- ggplot(
std_betas,
aes(
x = estimate,
y = fct_reorder(term_label, estimate),
xmin = conf.low,
xmax = conf.high,
color = term_label
)
) +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
geom_errorbarh(height = 0.25, linewidth = 0.8) +
geom_point(size = 2.5) +
facet_wrap(~ outcome_label, scales = "free_x") +
scale_color_brewer(palette = "Dark2", guide = "none") +
labs(
title = "Standardized Associations Between LMC, Teacher Well-Being, and Outcomes",
subtitle = "Points show standardized β coefficients; bars indicate 95% confidence intervals",
x = "Standardized Effect Size (β)",
y = NULL
) +
theme_pub(11)
p_std
# save_pub(p_std, "Fig_LMC_TWB_Standardized_Betas.png", width = 11, height = 7)Across the full system, labor–management collaboration (LMC) is consistently associated with better outcomes for educators, including higher job satisfaction, stronger retention intentions, and greater advocacy for schools and districts. However, standardized analyses reveal that LMC operates primarily as a contextual enabler rather than a direct driver of well-being. The strongest and most consistent predictors of positive outcomes are teachers’ foundational well-being supports—such as feeling supported, resourced, and psychologically safe—and lower levels of depletion. Even in highly collaborative districts, outcomes suffer when educators are exhausted or lack foundational support. In contrast, districts that pair strong LMC with robust well-being supports show markedly better outcomes across satisfaction, stress, retention, and Net Promoter scores. These findings suggest that collaboration creates the conditions for success, but investments in teacher well-being determine whether those conditions translate into sustained impact.
library(tidyverse)
library(broom)
# ----------------------------
# Inputs (match your prior chunk)
# ----------------------------
outcomes <- c("job_sat","job_stress","retention_score","nps_score")
predictors <- c(
"lmc_composite_score",
"twb_leader_score_short8",
"foundational_support_score",
"growth_score",
"wellbeing_score",
"acceptance_score",
"adaptability_score",
"depletion_score"
)
pretty_terms <- c(
lmc_composite_score = "LMC Composite",
twb_leader_score_short8 = "TWB Leadership",
foundational_support_score = "TWB Foundational Supports",
growth_score = "TWB Growth Orientation",
wellbeing_score = "TWB Personal Well-being",
acceptance_score = "TWB Acceptance",
adaptability_score = "TWB Adaptability",
depletion_score = "TWB Depletion"
)
pretty_outcomes <- c(
job_sat = "Job Satisfaction",
job_stress = "Job Stress",
retention_score = "Retention Intentions",
nps_score = "Net Promoter Score"
)
# ----------------------------
# Standardize (complete-case across all vars used)
# ----------------------------
z <- function(x) as.numeric(scale(x))
df_z <- df %>%
select(all_of(outcomes), all_of(predictors)) %>%
drop_na() %>%
mutate(across(everything(), z))
# ----------------------------
# Fit standardized models + extract with CIs
# ----------------------------
std_betas <- map_dfr(outcomes, function(out) {
f <- as.formula(paste(out, "~", paste(predictors, collapse = " + ")))
lm(f, data = df_z) %>%
tidy(conf.int = TRUE) %>%
filter(term != "(Intercept)") %>%
mutate(
outcome = out,
outcome_label = pretty_outcomes[out],
term_label = pretty_terms[term],
sig = case_when(
p.value < .001 ~ "***",
p.value < .01 ~ "**",
p.value < .05 ~ "*",
p.value < .10 ~ "†",
TRUE ~ ""
),
beta_lbl = paste0(
sprintf("%.2f", estimate),
sig
)
)
})
# ----------------------------
# Plot with numeric beta labels ON the points
# ----------------------------
p_std_labeled <- ggplot(
std_betas,
aes(
x = estimate,
y = fct_reorder(term_label, estimate),
xmin = conf.low,
xmax = conf.high
)
) +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
geom_errorbarh(height = 0.25, linewidth = 0.8) +
geom_point(size = 2.5) +
geom_text(
aes(label = beta_lbl),
nudge_x = 0.02,
size = 3.2
) +
facet_wrap(~ outcome_label, scales = "free_x") +
labs(
title = "Standardized Associations (β) Between LMC, TWB, and Outcomes",
subtitle = "Numbers show standardized β (with significance). Bars indicate 95% confidence intervals.",
x = "Standardized Effect Size (β)",
y = NULL
) +
theme_pub(11)
p_std_labeled
# save_pub(p_std_labeled, "Fig_LMC_TWB_Standardized_Betas_Labeled.png", width = 11, height = 7)
# ----------------------------
# Optional: print a clean table version too
# ----------------------------
std_betas %>%
select(outcome_label, term_label, estimate, conf.low, conf.high, p.value) %>%
mutate(
across(c(estimate, conf.low, conf.high), ~ round(.x, 3)),
p_value_fmt = case_when(
p.value < .001 ~ "<.001",
TRUE ~ sprintf("%.3f", p.value)
)
) %>%
arrange(outcome_label, desc(abs(estimate)))## # A tibble: 32 × 7
## outcome_label term_label estimate conf.low conf.high p.value p_value_fmt
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Job Satisfaction TWB Found… 0.37 0.345 0.395 3.33e-171 <.001
## 2 Job Satisfaction TWB Leade… 0.232 0.208 0.256 3.86e- 76 <.001
## 3 Job Satisfaction TWB Deple… -0.158 -0.18 -0.136 5.58e- 44 <.001
## 4 Job Satisfaction TWB Growt… 0.089 0.066 0.112 3.11e- 14 <.001
## 5 Job Satisfaction LMC Compo… 0.082 0.059 0.104 2.53e- 12 <.001
## 6 Job Satisfaction TWB Accep… 0.035 0.009 0.062 8.69e- 3 0.009
## 7 Job Satisfaction TWB Perso… 0.016 -0.007 0.038 1.74e- 1 0.174
## 8 Job Satisfaction TWB Adapt… 0.001 -0.026 0.027 9.61e- 1 0.961
## 9 Job Stress TWB Deple… 0.408 0.382 0.434 2.18e-192 <.001
## 10 Job Stress TWB Found… -0.221 -0.25 -0.191 3.61e- 48 <.001
## # ℹ 22 more rows