# Load and prepare data 
suppressPackageStartupMessages(library(dplyr))
library(dplyr)
library(haven)
library(knitr)
library(ggplot2)
library(tidyr)
data_full <- read_dta("C:/Users/JH/OneDrive - Kitces.com/Desktop/JH/Research_Ticket/Wellbeing Senior Advisors/Stata_WellBeing_2025/Stata/Master_Well_Being_2025.dta")

# Filter
senior_advisors <- data_full %>% filter(rspbizrole == 2)

# Table 1: Distribution of Compensation Mix
senior_advisors %>%
 mutate(rsppaystruc = as_factor(rsppaystruc)) %>%   
 count(rsppaystruc, name = "Count") %>%
 mutate(Percent = round(Count / sum(Count) * 100, 1)
 ) %>%
 arrange(rsppaystruc) %>%
 kable(
 caption = "Distribution of Compensation Mix - Senior Advisors",
 col.names = c("Compensation Mix", "Count", "Percent (%)"),
 align = c("l", "r", "r")
 )
Distribution of Compensation Mix - Senior Advisors
Compensation Mix Count Percent (%)
Fixed salary only 32 4.9
Salary as draw 13 2.0
Salary with revenue incentive 109 16.8
Salary with non-revenue incentive 37 5.7
Revenue-based only 124 19.1
Net profits only 301 46.5
Other 31 4.8
NA 1 0.2 
knitr::asis_output("Senior advisor compensation is concentrated in variable pay structures. Nearly half of senior advisors are compensated solely through net profits, with another one-fifth relying on revenue-based pay only, indicating substantial income variability at the senior level. Hybrid salary-plus-incentive models are present but less common, while purely fixed compensation is relatively rare.")

Senior advisor compensation is concentrated in variable pay structures. Nearly half of senior advisors are compensated solely through net profits, with another one-fifth relying on revenue-based pay only, indicating substantial income variability at the senior level. Hybrid salary-plus-incentive models are present but less common, while purely fixed compensation is relatively rare.

# Table 2: Well-Being by Compensation Mix for Senior Advisors
senior_advisors %>%
 mutate(rsppaystruc = as_factor(rsppaystruc)) %>%   
 group_by(rsppaystruc) %>%
 summarise(Mean_Wellbeing = round(mean(cantril_prsnt, na.rm = TRUE), 2),
 N = n(),
 .groups = "drop"
 ) %>%
 arrange(rsppaystruc) %>%
 kable(
 col.names = c("Compensation Mix", "Mean Wellbeing", "Number of Advisors"),
 caption = "Average Present Cantril Ladder Score by Compensation Structure - Senior Advisors",
 align = c("l", "r", "r")
 )
Average Present Cantril Ladder Score by Compensation Structure - Senior Advisors
Compensation Mix Mean Wellbeing Number of Advisors
Fixed salary only 7.12 32
Salary as draw 6.69 13
Salary with revenue incentive 7.43 109
Salary with non-revenue incentive 7.11 37
Revenue-based only 7.31 124
Net profits only 7.64 301
Other 7.29 31
NA 8.00 1 
knitr::asis_output("Average wellbeing scores are relatively high and tightly clustered across compensation structures, ranging roughly from 6.7 to 7.6 on the Cantril ladder. Senior advisors with profit-based and revenue-based compensation report well-being levels comparable to, or slightly higher than, those with salary-based arrangements. Notably, there is no clear well-being penalty associated with more variable compensation structures.")

Average wellbeing scores are relatively high and tightly clustered across compensation structures, ranging roughly from 6.7 to 7.6 on the Cantril ladder. Senior advisors with profit-based and revenue-based compensation report well-being levels comparable to, or slightly higher than, those with salary-based arrangements. Notably, there is no clear well-being penalty associated with more variable compensation structures.

# Table 3: Well-Being by Client-Facing Experience
senior_advisors %>%
 mutate(grpclexp = as_factor(grpclexp)) %>%    
 group_by(grpclexp) %>%
 summarise(
 Mean_Wellbeing = round(mean(cantril_prsnt, na.rm = TRUE), 2),
 N = n(),
 .groups = "drop"
 ) %>%
 arrange(grpclexp) %>%
 kable(
 col.names = c("Client-Facing Experience", "Mean Wellbeing", "Number of Advisors"),
 caption = "Wellbeing by Client-Facing Experience - Senior Advisors",
 align = c("l", "r", "r")
 )
Wellbeing by Client-Facing Experience - Senior Advisors
Client-Facing Experience Mean Wellbeing Number of Advisors
<5 yrs 7.27 65
5 to 9 6.99 112
10 to 19 7.35 189
20+ yrs 7.75 274
NA 7.50 8 
knitr::asis_output("Average wellbeing increases modestly with client-facing experience, with the lowest scores among advisors with 5-9 years of experience and the highest among those with 20+ years. Advisors earlier in their careers report slightly lower wellbeing, while more experienced advisors appear to benefit from greater stability or role fit.")

Average wellbeing increases modestly with client-facing experience, with the lowest scores among advisors with 5-9 years of experience and the highest among those with 20+ years. Advisors earlier in their careers report slightly lower wellbeing, while more experienced advisors appear to benefit from greater stability or role fit.

# Table 4: Well-Being by Team Structure / Support Level
senior_advisors %>%
 mutate(pstructure = as_factor(pstructure)) %>%  
 group_by(pstructure) %>%
 summarise(
  Mean_Wellbeing = round(mean(cantril_prsnt, na.rm = TRUE), 2),
  N = n(),
  .groups = "drop"
  ) %>%
  arrange(pstructure) %>%
  kable(
 col.names = c("Team Structure / Support Level", "Mean Wellbeing", "Number of Advisors"),
  caption = "Wellbeing by Team Structure - Senior Advisors",
 align = c("l", "r", "r")
 )
Wellbeing by Team Structure - Senior Advisors
Team Structure / Support Level Mean Wellbeing Number of Advisors
Unsupported Solo 7.29 180
Supported Solo 7.45 169
Silo 7.42 88
Ensemble 7.60 211 
knitr::asis_output("Wellbeing is lowest among unsupported solo advisors and highest among those in ensemble structures. Supported solo advisors report wellbeing levels comparable to ensemble firms, suggesting that access to support matters more than firm size alone. These patterns indicate that team context and structural support are meaningfully associated with senior advisor wellbeing.")

Wellbeing is lowest among unsupported solo advisors and highest among those in ensemble structures. Supported solo advisors report wellbeing levels comparable to ensemble firms, suggesting that access to support matters more than firm size alone. These patterns indicate that team context and structural support are meaningfully associated with senior advisor wellbeing.

# Table 5: Well-Being by Ownership Status (Partner/Owner vs Employee only) 
senior_advisors %>%
 filter(rspstatus %in% c(1, 2)) %>%  
 mutate(rspstatus = as_factor(rspstatus)) %>%
 group_by(rspstatus) %>%
 summarise(
 Mean_Wellbeing = round(mean(cantril_prsnt, na.rm = TRUE), 2),
 N = n(),
 .groups = "drop"
 ) %>%
 arrange(rspstatus) %>%
 kable(
 col.names = c("Ownership Status", "Mean Well-Being", "Number of Advisors"),
 caption = "Wellbeing by Ownership Status - Senior Advisors",
 align = c("l", "r", "r")
 )
Wellbeing by Ownership Status - Senior Advisors
Ownership Status Mean Well-Being Number of Advisors
Partner/Owner 7.56 517
Employee 7.03 102 
knitr::asis_output("Owner/Partner senior advisors report higher well-being compared to Employees. This suggests ownership may contribute to greater satisfaction, purpose, or control among senior advisors.")

Owner/Partner senior advisors report higher well-being compared to Employees. This suggests ownership may contribute to greater satisfaction, purpose, or control among senior advisors.

# Chart: Well-Being by Ownership Status and Experience Group  
 senior_advisors %>%
 filter(rspstatus %in% c(1, 2)) %>%
 mutate(
rspstatus = as_factor(rspstatus),
exp_group = case_when(
rspfsexp <= 5 ~ "0–5 years",
rspfsexp <= 10 ~ "6–10 years",
rspfsexp <= 15 ~ "11–15 years",
TRUE ~ ">15 years"
 ),
exp_group = factor(exp_group, levels = c("0–5 years", "6–10 years", "11–15 years", ">15 years"))
 ) %>%
 group_by(rspstatus, exp_group) %>%
 summarise(
Mean_Wellbeing = round(mean(cantril_prsnt, na.rm = TRUE), 2),
N = n(),
.groups = "drop"
) %>%
ggplot(aes(x = exp_group, y = Mean_Wellbeing, fill = rspstatus)) +
geom_bar(stat = "identity", position = "dodge", width = 0.65) +
geom_text(aes(label = Mean_Wellbeing),
position = position_dodge(width = 0.65),
vjust = 1.2,
size = 2.8,
color = "white",
fontface = "bold") +
labs(
title = "Wellbeing by Ownership Status - Senior Advisors",
x = "Years of Experience",
y = "Mean Cantril Ladder Score (Present)",
fill = ""
) +
scale_fill_manual(values = c("Partner/Owner" = "#001F3F",
"Employee" = "#007BFF")) +
theme_minimal(base_size = 10) +
theme(
plot.title = element_text(hjust = 0, 
colour = "grey50", 
size = 10, 
face = "plain"),
axis.title = element_text(face = "bold"),
axis.text = element_text(size = 10),
legend.position = "bottom",
legend.title = element_text(size = 10, face = "bold"),
legend.text = element_text(size = 10),
legend.key.size = unit(1.3, "lines"),
legend.background = element_rect(fill = "white", colour = "grey80"),
plot.margin = margin(10, 10, 10, 10)
)

knitr::asis_output("Partners/Owners consistently report higher wellbeing than Employees across all experience groups. The advantage is most pronounced in the early career stage and remains stable in later years. This suggests ownership may provide greater purpose, autonomy, or financial security that supports life satisfaction from the start of a senior advisor's career.")

Partners/Owners consistently report higher wellbeing than Employees across all experience groups. The advantage is most pronounced in the early career stage and remains stable in later years. This suggests ownership may provide greater purpose, autonomy, or financial security that supports life satisfaction from the start of a senior advisor’s career.

# Table 6: Wellbeing by Compensation Structure – by Ownership
senior_advisors %>%
filter(rspstatus %in% c(1, 2)) %>%
mutate(
rsppaystruc = as_factor(rsppaystruc),
rspstatus = as_factor(rspstatus)
 ) %>%
 group_by(rsppaystruc, rspstatus) %>%
 summarise(Mean_Wellbeing = round(mean(cantril_prsnt, na.rm = TRUE), 2), .groups = "drop") %>%
 pivot_wider(
 names_from = rspstatus,
 values_from = Mean_Wellbeing,
 names_prefix = "Mean_Wellbeing_"
 ) %>%
 arrange(rsppaystruc) %>%
 kable(
 col.names = c("Compensation Mix", "Partner/Owner", "Employee"),
 caption = "Mean Wellbeing by Compensation Structure and Ownership Status - Senior Advisors",
 align  = c("l", "r", "r"),
 digits = 2
 )
Mean Wellbeing by Compensation Structure and Ownership Status - Senior Advisors
Compensation Mix Partner/Owner Employee
Fixed salary only 7.37 6.77
Salary as draw 6.80 6.33
Salary with revenue incentive 7.64 7.03
Salary with non-revenue incentive 7.55 6.59
Revenue-based only 7.41 7.44
Net profits only 7.64 NA
Other 7.34 NA
NA 8.00 NA 
knitr::asis_output("Owners generally report higher wellbeing than employees across most compensation structures, with the largest gaps in 'Salary with non-revenue incentive'. This suggests ownership may enhance satisfaction in variable-pay environments.")

Owners generally report higher wellbeing than employees across most compensation structures, with the largest gaps in ‘Salary with non-revenue incentive’. This suggests ownership may enhance satisfaction in variable-pay environments.

# Wellbeing by Total Compensation Level and Ownership Status
senior_advisors %>%
 filter(rspstatus %in% c(1, 2), !is.na(rspinc)) %>%
 mutate(
 comp_level = case_when(
 rspinc < 150000 ~ "< $150,000",
 rspinc < 250000 ~ "$150,000 – $249,999",
 rspinc < 400000 ~ "$250,000 – $399,999",
 rspinc >= 400000  ~ "$400,000+"
 ),
 comp_level = factor(comp_level, levels = c("< $150,000", "$150,000 – $249,999", "$250,000 – $399,999", "$400,000+")),
rspstatus = as_factor(rspstatus)
 ) %>%
 group_by(comp_level, rspstatus) %>%
 summarise(Mean_Wellbeing = mean(cantril_prsnt, na.rm = TRUE), .groups = "drop") %>%
 ggplot(aes(x = comp_level, y = Mean_Wellbeing, fill = rspstatus)) +
 geom_bar(stat = "identity", position = "dodge", width = 0.7) +
 geom_text(aes(label = round(Mean_Wellbeing, 2)),
 position = position_dodge(0.7), vjust = -0.4, size = 3) +
 labs(
 title = "Wellbeing by Total Compensation Level and Ownership Status - Senior Advisors",
 x  = "Total Compensation (rspinc)",
 y = "Mean Cantril Ladder Score (Present)",
 fill = "Ownership Status"
 ) +
 scale_fill_manual(values = c("Partner/Owner" = "#001F3F", "Employee" = "#007BFF")) +
 theme_minimal(base_size = 10) +
 theme(
 plot.title = element_text(hjust = 0, colour = "grey50", size = 10, face = "plain"),
 axis.title = element_text(face = "bold"),
 axis.text.x = element_text(angle = 15, hjust = 1)
 )

knitr::asis_output("Across all compensation levels, wellbeing rises consistently with higher total compensation for both ownership statuses. Partner/Owners tend to score marginally higher than Employees at each compensation band, though the gap remains small.")

Across all compensation levels, wellbeing rises consistently with higher total compensation for both ownership statuses. Partner/Owners tend to score marginally higher than Employees at each compensation band, though the gap remains small.

# Wellbeing by Comp Structure vs. Total Comp Level (Owners/Partners only – 2x2 comparison)
owners_only <- senior_advisors %>%
filter(rspstatus == 1) %>% 
 mutate(
 structure_group = case_when(
 rsppaystruc %in% c(1, 2) ~ "Exclusively Variable\n(Profits or Revenue Only)",
!is.na(rsppaystruc) ~ "At Least Partially Fixed",
 TRUE ~ "Missing"
 ),
 structure_group = factor(structure_group, levels = c(
 "Exclusively Variable\n(Profits or Revenue Only)",
 "At Least Partially Fixed"
 ))
 ) %>%
 filter(!is.na(structure_group), !is.na(rspinc), !is.na(cantril_prsnt))

# Compute median among this sample
median_comp <- median(owners_only$rspinc, na.rm = TRUE)
median_label <- paste0("$", round(median_comp / 1000), "k")  

owners_plot_data <- owners_only %>%
 mutate(
 comp_group = if_else(rspinc >= median_comp, "Above Median", "Below Median"),
 comp_group = factor(comp_group, levels = c("Below Median", "Above Median"))
 ) %>%
 group_by(structure_group, comp_group) %>%
 summarise(
 Mean_Wellbeing = mean(cantril_prsnt, na.rm = TRUE),
 N = n(),
 .groups = "drop"
 )

# The chart
owners_plot_data %>%
 ggplot(aes(x = structure_group, y = Mean_Wellbeing, fill = comp_group)) +
 geom_bar(stat = "identity", position = "dodge", width = 0.65) +
 geom_text(aes(label = round(Mean_Wellbeing, 2)),
 position = position_dodge(width = 0.65),
 vjust = -0.4, size = 3.2) +
 labs(
 title = paste0("Wellbeing by Compensation Structure vs. Total Comp Level\n(Owners/Partners Only; Median Comp = ", median_label, ")"),
 x = "Compensation Structure Group",
 y = "Mean Cantril Ladder Score (Present)",
 fill = ""
 ) +
 scale_fill_manual(values = c("Below Median" = "#007BFF", "Above Median" = "#001F3F")) +
 theme_minimal(base_size = 10) +
 theme(
 plot.title = element_text(hjust = 0, colour = "grey50", size = 10, face = "plain"),
 axis.title = element_text(face = "bold"),
 axis.text.x = element_text(angle = 10, hjust = 1, size = 9),
 legend.position = "bottom"
 )

knitr::asis_output("Among owner/partners, higher total compensation is associated with modestly higher wellbeing in both structure groups. Absolute compensation level seems to matter at least as much as, if not more than, purely variable vs. mixed/fixed payout structure for this subgroup.")

Among owner/partners, higher total compensation is associated with modestly higher wellbeing in both structure groups. Absolute compensation level seems to matter at least as much as, if not more than, purely variable vs. mixed/fixed payout structure for this subgroup.

# Chart: Wellbeing by Comp Structure vs. Total Comp Level
# (Non-Owner / Employee Senior Advisors Only)

nonowners_only <- data_full %>%
 filter(rspbizrole == 2, rspstatus == 2) %>% 
 mutate(
 structure_group = case_when(
 rsppaystruc %in% c(1, 2) ~ "Exclusively Variable\n(Profits or Revenue Only)",
!is.na(rsppaystruc) ~ "At Least Partially Fixed",
 TRUE ~ "Missing"
 ),
 structure_group = factor(structure_group, levels = c(
 "Exclusively Variable\n(Profits or Revenue Only)",
 "At Least Partially Fixed"
 ))
 ) %>%
 filter(!is.na(structure_group), !is.na(rspinc), !is.na(cantril_prsnt))

if (nrow(nonowners_only) == 0) {
 print("No data after filtering — check rspstatus coding or missing values")
} else {
 median_comp_nonowners <- median(nonowners_only$rspinc, na.rm = TRUE)
 median_label_nonowners <- paste0("$", round(median_comp_nonowners / 1000), "k")

 nonowners_plot_data <- nonowners_only %>%
 mutate(
 comp_group = if_else(rspinc >= median_comp_nonowners, "Above Median", "Below Median"),
 comp_group = factor(comp_group, levels = c("Below Median", "Above Median"))
 ) %>%
 group_by(structure_group, comp_group) %>%
 summarise(
 Mean_Wellbeing = mean(cantril_prsnt, na.rm = TRUE),
N = n(),
 .groups = "drop"
 )

 # The chart
 nonowners_plot_data %>%
 ggplot(aes(x = structure_group, y = Mean_Wellbeing, fill = comp_group)) +
 geom_bar(stat = "identity", position = "dodge", width = 0.65) +
 geom_text(aes(label = round(Mean_Wellbeing, 2)),
 position = position_dodge(width = 0.65),
 vjust = -0.4, size = 3.2) +
 labs(
title = paste0("Wellbeing by Compensation Structure vs. Total Comp Level\n(Non-Owner / Employee Senior Advisors Only; Median Comp = ", median_label_nonowners, ")"),
 x = "Compensation Structure Group",
 y = "Mean Cantril Ladder Score (Present)",
 fill = ""
 ) +
 scale_fill_manual(values = c("Below Median" = "#007BFF", "Above Median" = "#001F3F")) +
 theme_minimal(base_size = 10) +
 theme(
 plot.title = element_text(hjust = 0, colour = "grey50", size = 10, face = "plain"),
 axis.title = element_text(face = "bold"),
 axis.text.x = element_text(angle = 10, hjust = 1, size = 9),
 legend.position = "bottom"
 )

}

knitr::asis_output("For employee senior advisors, wellbeing increases with total compensation above the median in both structure groups. The gain is more pronounced in exclusively variable pay (6.00 to 7.14) than in partially fixed structures (6.76 to 7.40), suggesting that for non-owners, the amount of compensation may matter more than whether it is fixed or variable.")

For employee senior advisors, wellbeing increases with total compensation above the median in both structure groups. The gain is more pronounced in exclusively variable pay (6.00 to 7.14) than in partially fixed structures (6.76 to 7.40), suggesting that for non-owners, the amount of compensation may matter more than whether it is fixed or variable.