Who’s Helping Whom? How Organizational Social Support Hierarchies Shape Mental Health Outcomes

Author

Joe Quinn

Published

August 4, 2025

Program Overview

Goals of this program:

  • Explore the attribute-based data

  • Make and document analysis sample decisions

  • Univariate descriptives for DV

  • Univariate descriptives for the nodal attribute that goes into our IV

  • Bivariate descriptives for things that might be correlated with both our IV and DV

  • Bivariate models to test the idea (sans ties themselves)

Up next (before I leave today, add):

  • Reconstruct dyadic measures from full network data

  • Construct nodal measures from full network data

  • Run permutation tests

  • Append measures to analysis sample file

  • Bivariate network regression models

Read in and prepare the CHI data

Done - section’s code hidden from reader for convenience.

Filter the Sample

We have talked about some exclusion criteria for nodes. We want to keep:

  1. Non-retired clergy
  2. Full-time clergy
  3. Clergy only in the NC Annual and Western NC Annual conferences.

I take care of that below. I capture all discussion about this by creating an analysis sample based on those decisions. Changes to the analysis frame (and different files we build for robustness checks) should appear here too.

Note: as of 8/4/25, we only analyze 2021 observations (see chunk above). We should bring in 2023 data as well. Comeback.

#|echo: false

# pre-cleaning peep
table(n.attrs$ax1, n.attrs$a2c, deparse.level = 2)
                                          n.attrs$a2c
n.attrs$ax1                                1/4 time 1/2 time 3/4 time Full time
  Retired, under episcopal appointment           54       61       12        12
  Retired, not under episcopal appointment        0        0        0         0
  No, not retired                                79       97       42       880
# keep only full-time and not retired clergy
n.attrs.ansamp <- n.attrs |> 
  dplyr::filter(ax1 == "No, not retired") |>
  dplyr::filter(a2c == "Full time") |>
  dplyr::filter(ax11 %in% c("NC conf.","Western NC conf."))

# remove empty factor levels
n.attrs.ansamp <- n.attrs.ansamp |> 
  dplyr::mutate(across(c(ax1, a2c, ax11), ~ droplevels(.)))

# post-cleaning peep
table(n.attrs.ansamp$ax1, n.attrs.ansamp$a2c, deparse.level = 2)
                  n.attrs.ansamp$a2c
n.attrs.ansamp$ax1 Full time
   No, not retired       859
table(n.attrs.ansamp$ax1, n.attrs.ansamp$ax11, deparse.level = 2)
                  n.attrs.ansamp$ax11
n.attrs.ansamp$ax1 NC conf. Western NC conf.
   No, not retired      321              538

Our survey analysis sample includes 859 full-time, non-retired UMC clergy who do their work in one of the two annual conferences in NC.

A few notes:

First, by removing non-FTE clergy, we exclude many provisional and student clergy. We also make ansamp estimates for such provisional / student clergy weird. For example, full-time local pastors who are also students are not “normal” for the LP category (or the student LP category). Because this filtering leads to non-representative-ness of clergy within their status groups, I think the survey analysis sample should exclude provisional / student clergy.

Second, when we create measures related to support-seeking ties via the network data (by contrast to #1), we should include UMC clergy who are not in the survey analysis sample. From the POV of the ego, the amount of support one seeks / provides is not constrained to the analysis sample, and determining whether an alter provides one with support is ego-level perception-based data. By keeping all nominees / nominations when computing support network measures, we are not keeping non-representative cases like we would in #1; we are just counting support-seeking volume and support-providing load correctly. Especially because it is likely that non-trainees will provide more support to trainees (removing them would underestimate support-giving volume of the clergy in the analysis sample).

In other words, we should use all of the 2021 nodes when computing node-level measures about connectivity and edge info, but we should only analyze the analysis sample respondents in regression-style analyses we run after creating those full-graph measures.

So let’s add another exclusion criterion here:

  1. No student and provisional clergy

Remove non-terminal roles from the ansamp file here (we will still keep the edge attrs file for assigning connectivity features to each of these provisional clergy).

#|echo: false

# keep only non-training statuses
n.attrs.ansamp <- n.attrs.ansamp |>
  filter(ax10 %in% c("Elder", "LP", "Deacon")) |>
   droplevels()

# check the filter
table(n.attrs.ansamp$ax10, deparse.level = 2)
n.attrs.ansamp$ax10
 Elder     LP Deacon 
   590    140     57 
# peep n
nrow(n.attrs.ansamp)
[1] 787
# visualize breakdown
n.attrs.ansamp |> 
  ggplot(aes(x = ax10)) + geom_bar() + coord_flip()

We now have an analysis file with 787 full-time, non-retired, non-trainee clergy. This includes 590 elders, 140 local pastors, and 57 deacons.

Notes:

  1. Discuss with RJP: is this the right set of choices for dealing with non-FTE, non-traininee, and non-retired actors analytically? They will show up in the network even when they are not in our analysis sample. We still need info on their position in the hierarchy in order to estimate the effects of these hierarchical ties on any outcome, not merely their existence as a tie (i.e., we need their nodal attributes but they’re not in the analysis sample).
  2. These are currently 2021 respondents only - #comeback.
  3. Most clergy are elders, so we may not have a ton of contrast. We will add district superintendent info to facilitate (also - we should consider an “ever district super” flag, as the position likely alters one’s network even after they rotate out of the role).

#comeback when we get a clean DS info file with UIDs and years.

#|echo: false
#|include: false

# comeback.

Cross Clergy Status with Other Traits

Here I explore clergy status by covariates of interest so we can model relationships in the data responsibly (avoid/note small cell sizes; consider other vars associated with our predictors; etc.)

Gender

#|echo: false

# Cross-tab of clergy type vs. gender
ctable(
  x       = n.attrs.ansamp$ax10,
  y       = n.attrs.ansamp$female,
  prop    = "r", 
  useNA   = "no",
  totals  = TRUE,
  chisq   = TRUE,
  headings= TRUE
)
Cross-Tabulation, Row Proportions  
ax10 * female  
Data Frame: n.attrs.ansamp  


-------- -------- ------------- ------------- --------------
           female        Female          Male          Total
    ax10                                                    
   Elder            198 (33.6%)   392 (66.4%)   590 (100.0%)
      LP             40 (28.6%)   100 (71.4%)   140 (100.0%)
  Deacon             44 (77.2%)    13 (22.8%)    57 (100.0%)
   Total            282 (35.8%)   505 (64.2%)   787 (100.0%)
-------- -------- ------------- ------------- --------------

----------------------------
 Chi.squared   df   p.value 
------------- ---- ---------
   46.9449     2       0    
----------------------------

Female R’s make up only 35.8% of clergy in the data.

Gender and clergy status are not independent. Most elders and local pastors are men (66.4% and 71.4%, respectively. Most deacons are women (77.2%).

Conference

#|echo: false

# Cross-tab of clergy type vs. conference
ctable(
  x       = n.attrs.ansamp$ax11,
  y       = n.attrs.ansamp$ax10,
  prop    = "r",
  useNA   = "no",
  totals  = F,
  chisq   = TRUE,
  headings= TRUE
)
Cross-Tabulation, Row Proportions  
ax11 * ax10  
Data Frame: n.attrs.ansamp  


------------------ ------ ------------- ------------ -----------
                     ax10         Elder           LP      Deacon
              ax11                                              
          NC conf.          220 (73.3%)   65 (21.7%)   15 (5.0%)
  Western NC conf.          370 (76.0%)   75 (15.4%)   42 (8.6%)
------------------ ------ ------------- ------------ -----------

----------------------------
 Chi.squared   df   p.value 
------------- ---- ---------
   7.6373      2     0.022  
----------------------------

More clergy in the ansamp are from the Western NC Annual Conference (n=487) than the NC Annual Conference (n=300). There is some evidence that clergy status and conference are not independent: a larger proportion of clergy in the WNCC conference are deacons, and a smaller proportion are local pastors, compared to the proportions of the NCC.

Ask Logan, Rae Jean, and David for requests on other ways to cut the data that may be consequential (things that are correlated with clergy status that we will need to account for, as we construct our key explanatory variable(s) from clergy status and need to address alternative explanations for our outcomes. So: #comeback.

Inspect Outcomes

The three outcomes of interest are distributed as such in the ansamp:

#|echo: false

summary(n.attrs.ansamp$codi_mean)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  0.000   1.000   1.400   1.395   1.800   3.000      38 
summary(n.attrs.ansamp$min_satis_mean)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 0.4167  1.9167  2.2500  2.1861  2.5000  3.0000     141 
summary(n.attrs.ansamp$flourisher)
Not flourishing     Flourishing            NA's 
            307             462              18 
n.attrs.ansamp |> 
  ggplot(aes(x = codi_mean)) + 
  geom_histogram(binwidth = 0.2, fill = "lightblue")
Warning: Removed 38 rows containing non-finite outside the scale range
(`stat_bin()`).

n.attrs.ansamp |> 
  ggplot(aes(x = min_satis_mean)) + 
  geom_histogram()
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Warning: Removed 141 rows containing non-finite outside the scale range
(`stat_bin()`).

n.attrs.ansamp |> 
  ggplot(aes(x = flourisher)) + 
  geom_bar()

  • codi_mean is normally distributed and ranges (empirically) from 0 to 3; 38 obs are missing.

  • min_satis_mean is left-skewed and ranges (empirically) from ~0.42 to 3; 141 obs are missing.

  • flourisher is a skewed proportion; 18 obs are missing

Note: ask RJP and LT about the missing-ness in min_satisf_mean. Could this be a balloting thing, or is the scale composed of items that people are more likely to skip?

Bivariate - CODI x Clergy Status

Check out the bivariate relationship between occupational distress and clergy type. Assuming that higher-status clergy provide support to more people (which we can measure materially, not merely the vibe like other work in this space), the orgs literature predicts that Elders should have the highest CODI, followed by Deacons, followed by LPs.

Why? Because hierarchically higher roles tend to be accompanied with more support-giving, which can benefit individuals who receive the support but harm those who provide it. In this line of thinking, if the above is true (i.e., if CODI scores are Elders > Deacons > LPs), then it is true in part because support burden also decreases in the same order.

Likewise, it also means that within a category, those with more support-giving nominations will have a lower CODI (e.g., a Deacon nominated by one person as a support tie will have a lower CODI than a Deacon nominated by four people as a support tie). This deductive statement is central to the question we explore in our paper.

m_codi <- aov(codi_mean ~ ax10, data = n.attrs.ansamp)
Anova(m_codi, type = "II") 
Anova Table (Type II tests)

Response: codi_mean
           Sum Sq  Df F value   Pr(>F)   
ax10        4.517   2  5.8742 0.002943 **
Residuals 286.826 746                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
v_codi <- n.attrs.ansamp |>
  group_by(ax10) |>
  summarise(
    mean_val = mean(codi_mean, na.rm = TRUE),
    se = sd(codi_mean, na.rm = TRUE) / sqrt(n()))

ggplot(v_codi, aes(x = ax10, y = mean_val)) +
  geom_col(fill = "steelblue") +
  geom_errorbar(aes(ymin = mean_val - se, ymax = mean_val + se), width = 0.2) +
  labs(
    x = "Clergy Status Groups",
    y = "codi_mean",
    title = "codi_mean estimates by clergy status"
  ) +
  scale_y_continuous(limits = c(0, 2)) +
  theme_minimal(base_size = 14)

A clergy-member’s CODI is dependent upon their occupational status. Elders have the highest CODI (clearly higher than LPs and Deacons). Deacons have the next highest CODI - clearly lower than elders, and higher (but uncertainly so) than local pastors.

While some people on our study team shared thoughts that local pastors may feel more distressed because of their lower levels of connectivity to the rest of the clergy support network, it appears that the opposite might be true: LP well-being may be protected by their minimal connectivity to this support network (TBD - we need to re-establish these differences in connectivity; see the CHI support ties memo).

If the LP prediction was based on some indicator from the literature, this could be presented as a competing hypothesis (a “do support ties help or hurt?” question).

Outcome 2: Satisfaction

m_satis <- aov(min_satis_mean ~ ax10, data = n.attrs.ansamp)
Anova(m_satis, type = "II") 
Anova Table (Type II tests)

Response: min_satis_mean
           Sum Sq  Df F value  Pr(>F)  
ax10        1.552   2  4.0333 0.01817 *
Residuals 123.751 643                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
v_satis <- n.attrs.ansamp |>
  group_by(ax10) |>
  summarise(
    mean_val = mean(min_satis_mean, na.rm = TRUE),
    se = sd(min_satis_mean, na.rm = TRUE) / sqrt(n()))

ggplot(v_satis, aes(x = ax10, y = mean_val)) +
  geom_col(fill = "steelblue") +
  geom_errorbar(aes(ymin = mean_val - se, ymax = mean_val + se), width = 0.2) +
  labs(
    x = "Clergy Status Groups",
    y = "min_satis_mean",
    title = "min_satis_mean estimates by clergy status"
  ) +
  scale_y_continuous(limits = c(0, 3)) +   # Set y-axis limits
  theme_minimal(base_size = 14) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),  # Keep angled labels for readability
    plot.margin = margin(20, 20, 40, 20)
  )

A clergy-member’s min_satisf_mean score is dependent upon their occupational status. Local pastors have the highest min_satisf_mean (clearly higher than LPs and Deacons). Elders and deacons have comparatively lower min_satisf_mean scores (and they are not clearly distinct from each other).

This result complements the bi-variate results for clergy status x codi_mean above.

While some people on our study team shared thoughts that local pastors may feel less satisfied because of their lower levels of connectivity to the rest of the clergy support network, the opposite might be true: LP satisfaction may be protected by their minimal connectivity to this support network (e.g., less less support-giving).

Again, if the LP prediction was based on some indicator from the literature, this could be presented as a competing hypothesis (a “do support ties help or hurt?” question).

Note: 141 of the 787 observations in the analysis sample lack a min_satisf_mean score. If these are not missing at random, we should avoid analyzing this outcome. If they are, we should impute or we’ll probably have power issues.

Outcome 3: Flourishing

ctable(
  x       = n.attrs.ansamp$ax10,
  y       = n.attrs.ansamp$flourisher,
  prop    = "r", #r = row proportions
  useNA   = "no",
  totals  = FALSE,
  chisq   = TRUE,
  headings= TRUE
)
Cross-Tabulation, Row Proportions  
ax10 * flourisher  
Data Frame: n.attrs.ansamp  


-------- ------------ ----------------- -------------
           flourisher   Not flourishing   Flourishing
    ax10                                             
   Elder                    240 (41.8%)   334 (58.2%)
      LP                     46 (33.3%)    92 (66.7%)
  Deacon                     21 (36.8%)    36 (63.2%)
-------- ------------ ----------------- -------------

----------------------------
 Chi.squared   df   p.value 
------------- ---- ---------
   3.5779      2    0.1671  
----------------------------
# Compute proportion of flourishers by ax10
v_flourish <- n.attrs.ansamp |>
  filter(!is.na(flourisher)) |>
  group_by(ax10) |>
  summarise(
    n = n(),
    prop_flourish = mean(flourisher == "Flourishing", na.rm = TRUE)
  ) |>
  mutate(
    se = sqrt((prop_flourish * (1 - prop_flourish)) / n)
  )

ggplot(v_flourish, aes(x = ax10, y = prop_flourish)) +
  geom_col(fill = "steelblue") +
  geom_errorbar(aes(ymin = prop_flourish - se, ymax = prop_flourish + se),
                width = 0.2) +
  labs(
    x = "Clergy Status Groups",
    y = "Proportion Flourishing",
    title = "Proportion of Flourishers by Clergy Status"
  ) +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0, 1)) +
  theme_minimal(base_size = 14) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.margin = margin(20, 20, 40, 20)
  )

We lack strongevidence to suggest that clergy status and flourisher status are not independent (or that clergy status may have an “effect” on flourishing) in this super cursory analysis). This bivariate analysis reflects the patterns from the clergy status x codi_mean and clergy status x min_satisf_mean analyses: the plotted point estimates show that local pastors might be somewhat more likely to be flourishing than elders. And they reflect the general trend of flourishing proportions that are the inverse of the codi_mean trend (as they should be): local pastors > deacons > elders. But a chi-square test of these two categorical variables suggests that there is not sufficient evidence to reject the null hypothesis that they are independent (p = .16).

We should continue to explore this outcome in different ways, given its difference as a binary outcome at the individual level. For example:

#|echo: false

n.attrs.ansamp <- n.attrs.ansamp |> 
  mutate(
    elder  = if_else(ax10 == "Elder",  1L, 0L, missing = 0L),
    lp     = if_else(ax10 == "LP",     1L, 0L, missing = 0L),
    deacon = if_else(ax10 == "Deacon", 1L, 0L, missing = 0L)
  )

m1 <- glm(flourisher ~ lp + deacon, 
          data = n.attrs.ansamp, 
          family = binomial(link = "logit"))

summary(m1)

Call:
glm(formula = flourisher ~ lp + deacon, family = binomial(link = "logit"), 
    data = n.attrs.ansamp)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.33050    0.08462   3.906  9.4e-05 ***
lp           0.36265    0.19942   1.818    0.069 .  
deacon       0.20849    0.28733   0.726    0.468    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1034.6  on 768  degrees of freedom
Residual deviance: 1031.0  on 766  degrees of freedom
  (18 observations deleted due to missingness)
AIC: 1037

Number of Fisher Scoring iterations: 4
tab <- tidy(m1) |>
  mutate(
    odds_ratio  = exp(estimate), # e^β
    log_odds    = if_else(term == "(Intercept)",
                          estimate, # 'elders' is reference category
                          estimate + coef(m1)["(Intercept)"]),
    pred_prob   = plogis(log_odds), # logistic transformation
    comparison  = case_when(# intuitive labels
      term == "(Intercept)" ~ "Elders (reference)",
      term == "lp"          ~ "Local Pastors vs. Elders",
      term == "deacon"      ~ "Deacons vs. Elders"
    )
  ) |>
  select(comparison, estimate, odds_ratio, pred_prob, p.value) |> 
  rename(
    log_odds_coef = estimate,
    p_value       = p.value
  ) |>
  mutate(across(c(log_odds_coef, odds_ratio, pred_prob, p_value),
                ~ round(.x, 3)))

print(tab)
# A tibble: 3 × 5
  comparison               log_odds_coef odds_ratio pred_prob p_value
  <chr>                            <dbl>      <dbl>     <dbl>   <dbl>
1 Elders (reference)               0.331       1.39     0.582   0    
2 Local Pastors vs. Elders         0.363       1.44     0.667   0.069
3 Deacons vs. Elders               0.208       1.23     0.632   0.468

Ok - we have some evidence to suggest that local pastors might be more likely to get classified as flourishers than elders. This would be tricky to power if we add additional variables. And longitudinal data won’t really help us out with this sort of thing because these categories are largely time-invariant (especially as we’ve removed trainee categories).

Note: before excluding trainees / provisional clergy, the lowest point estimates for proportion of flourisher observations are clergy in the provisional deacon and provisional elder categories. The provisional categories are likely associated with support-seeking and social network ties; this might be worth exploring, but pause on this for now.