#|echo: false
n.attrs |>
ggplot(aes(x = ax10)) + geom_bar() + coord_flip()
Who’s Helping Whom? How Organizational Social Support Hierarchies Shape Mental Health Outcomes
Project Overview
…
Read in and prepare the CHI data
Explore Clergy Status Frequencies
Here’s the breakdown of the raw categories in ax10 - clergy status - which we will recode into our status hierarchy.
Most clergy are Full Elders, followed by Local Pastors, and then Deacons. The trainee versions of these categories proceed within the same order of relative frequency but contain far fewer R’s.
- NOTE: there are very few “Provisional Deacon” Rs. Breaking this category out may be an issue for modeling.
In the rest of this section, I explore clergy status by covariates of interest so we can avoid creating microscopic cells and get a feel for the data…
Gender
#|echo: false
# Simple 2×2 (or more) cross-tab of clergy status vs. conference
ctable(
x = n.attrs$ax10,
y = n.attrs$female,
prop = "r", #r = row proportions (use "c" or "t" for column/total)
useNA = "no",
totals = TRUE,
headings= TRUE
)Cross-Tabulation, Row Proportions
ax10 * female
Data Frame: n.attrs
----------------------------------------------- -------- ------------- ------------- ---------------
female Female Male Total
ax10
Provisional Elder 20 (35.7%) 36 (64.3%) 56 (100.0%)
Elder in Full Connection 295 (32.6%) 611 (67.4%) 906 (100.0%)
Student Local Pastor 15 (42.9%) 20 (57.1%) 35 (100.0%)
Local Pastor (but not a student local pastor) 111 (32.4%) 232 (67.6%) 343 (100.0%)
Provisional Deacon 10 (83.3%) 2 (16.7%) 12 (100.0%)
Deacon in Full Connection 69 (80.2%) 17 (19.8%) 86 (100.0%)
Associate Member 3 (17.6%) 14 (82.4%) 17 (100.0%)
Other (specify) 3 (30.0%) 7 (70.0%) 10 (100.0%)
Total 526 (35.9%) 939 (64.1%) 1465 (100.0%)
----------------------------------------------- -------- ------------- ------------- ---------------Female R’s make up 35.9% of clergy in the data. Within a category (Elder, Deacon, Local Pastor), their representation decreases from the trainee category to the full category (e.g., 32.4% of Local Pastors are female –> 42.9% of student local pastors are female). Female R’s are vastly over-represented as Deacons (83.3% –> 80.2%). They are slightly under-represented as Local Pastors (32.4% –> 42.9%) and Elders (32.6% –> 35.7%).
I can add other covariates of interest in this section. I’ll stop here for now.
Sample Filtering?
Our conversations have involved some exclusion criteria. I am capturing these exclusions here before nodal measure creation. I capture all discussion about this, but do not yet formally apply any filters to the analysis sample besides classifying “Associate Member” and “Other (specify)” ax10 responses as NA.
#|echo: false
# retired clergy
freq(n.attrs$ax1)Frequencies
n.attrs$ax1
Label: retired status
Type: Factor
Freq % Valid % Valid Cum. % Total % Total Cum.
---------------------------------------------- ------ --------- -------------- --------- --------------
Retired, under episcopal appointment 139 9.46 9.46 9.46 9.46
Retired, not under episcopal appointment 194 13.21 22.67 13.21 22.67
No, not retired 1136 77.33 100.00 77.33 100.00
<NA> 0 0.00 100.00
Total 1469 100.00 100.00 100.00 100.00Create Nodal Measures
Here we transform or create pastor-level measures for analysis based on conversations with RJP and DE. Some important notes:
- The measures built in this section currently reclassify all clergy in the “Associate Member” and “Other (specify)” categories as NA.
- We have previously talked about how years of experience may be highly correlated with aspects of this hierarchy. We have not yet explored this in more detail or identified a solution. I suspect this will more more related to item 3…
- We have thus far only discussed operationalizing the hierarchy as a single ladder (Elder-ish, Deacon-ish, Local Pastor-ish). The other way to think about this hierarchy is a set of three ladders with two categories each - a trainee subcategory and a full sub-category. This decision should be based on how RJP and DE think about the extent to which these categories are isolated or integrated, and Sohale’s review of the literature (e.g., other organizations have this sort of structure where there is a within-group hierarchy for administrative assistant and office manager that is not the same as analyst and vice president). I will create a trainee flag for now. More discussion is necessary.
#|echo: false
#|include: false
# Create three hierarchical indicators #
#--------------------------------------#
n.attrs <- n.attrs %>%
mutate(
# 1) three‐category grouping factor
clrgystatus_3cat = case_when(
ax10 %in% c("Provisional Elder", "Elder in Full Connection") ~ "Elder",
ax10 %in% c("Student Local Pastor",
"Local Pastor (but not a student local pastor)") ~ "Local Pastor",
ax10 %in% c("Provisional Deacon", "Deacon in Full Connection") ~ "Deacon",
TRUE ~ NA_character_
),
clrgystatus_3cat = factor(
clrgystatus_3cat,
levels = c("Elder", "Local Pastor", "Deacon")
),
# 2) six‐category detailed status factor
clrgystatus_6cat = case_when(
ax10 == "Provisional Elder" ~ "Elder - Provisional",
ax10 == "Elder in Full Connection" ~ "Elder - Full",
ax10 == "Student Local Pastor" ~ "Local Pastor - Student",
ax10 == "Local Pastor (but not a student local pastor)"
~ "Local Pastor - Full",
ax10 == "Provisional Deacon" ~ "Deacon - Provisional",
ax10 == "Deacon in Full Connection" ~ "Deacon - Full",
TRUE ~ NA_character_
),
clrgystatus_6cat = factor(
clrgystatus_6cat,
levels = c(
"Elder - Provisional", "Elder - Full",
"Local Pastor - Student", "Local Pastor - Full",
"Deacon - Provisional", "Deacon - Full"
)
),
# 3) trainee flag (0/1 integer)
trainee_flag = case_when(
ax10 %in% c("Provisional Elder", "Student Local Pastor",
"Provisional Deacon") ~ 1,
ax10 %in% c("Elder in Full Connection",
"Local Pastor (but not a student local pastor)",
"Deacon in Full Connection") ~ 0,
TRUE ~ NA_real_
)
)
# Check new vars
ctable(x = n.attrs$ax10,
y = n.attrs$clrgystatus_3cat,
prop = "r", #r = row proportions (use "c" or "t" for column/total)
useNA = "ifany",
totals = FALSE)Cross-Tabulation, Row Proportions
ax10 * clrgystatus_3cat
Data Frame: n.attrs
----------------------------------------------- ------------------ -------------- -------------- ------------- -------------
clrgystatus_3cat Elder Local Pastor Deacon <NA>
ax10
Provisional Elder 56 (100.0%) 0 ( 0.0%) 0 ( 0.0%) 0 ( 0.0%)
Elder in Full Connection 908 (100.0%) 0 ( 0.0%) 0 ( 0.0%) 0 ( 0.0%)
Student Local Pastor 0 ( 0.0%) 35 (100.0%) 0 ( 0.0%) 0 ( 0.0%)
Local Pastor (but not a student local pastor) 0 ( 0.0%) 345 (100.0%) 0 ( 0.0%) 0 ( 0.0%)
Provisional Deacon 0 ( 0.0%) 0 ( 0.0%) 12 (100.0%) 0 ( 0.0%)
Deacon in Full Connection 0 ( 0.0%) 0 ( 0.0%) 86 (100.0%) 0 ( 0.0%)
Associate Member 0 ( 0.0%) 0 ( 0.0%) 0 ( 0.0%) 17 (100.0%)
Other (specify) 0 ( 0.0%) 0 ( 0.0%) 0 ( 0.0%) 10 (100.0%)
----------------------------------------------- ------------------ -------------- -------------- ------------- -------------ctable(x = n.attrs$clrgystatus_6cat,
y = n.attrs$clrgystatus_3cat,
prop = "r", #r = row proportions (use "c" or "t" for column/total)
useNA = "ifany",
totals = FALSE)Cross-Tabulation, Row Proportions
clrgystatus_6cat * clrgystatus_3cat
Data Frame: n.attrs
------------------------ ------------------ -------------- -------------- ------------- -------------
clrgystatus_3cat Elder Local Pastor Deacon <NA>
clrgystatus_6cat
Elder - Provisional 56 (100.0%) 0 ( 0.0%) 0 ( 0.0%) 0 ( 0.0%)
Elder - Full 908 (100.0%) 0 ( 0.0%) 0 ( 0.0%) 0 ( 0.0%)
Local Pastor - Student 0 ( 0.0%) 35 (100.0%) 0 ( 0.0%) 0 ( 0.0%)
Local Pastor - Full 0 ( 0.0%) 345 (100.0%) 0 ( 0.0%) 0 ( 0.0%)
Deacon - Provisional 0 ( 0.0%) 0 ( 0.0%) 12 (100.0%) 0 ( 0.0%)
Deacon - Full 0 ( 0.0%) 0 ( 0.0%) 86 (100.0%) 0 ( 0.0%)
<NA> 0 ( 0.0%) 0 ( 0.0%) 0 ( 0.0%) 27 (100.0%)
------------------------ ------------------ -------------- -------------- ------------- -------------ctable(x = n.attrs$ax10,
y = n.attrs$trainee_flag,
prop = "r", #r = row proportions (use "c" or "t" for column/total)
useNA = "ifany",
totals = FALSE)Cross-Tabulation, Row Proportions
ax10 * trainee_flag
Data Frame: n.attrs
----------------------------------------------- -------------- -------------- ------------- -------------
trainee_flag 0 1 <NA>
ax10
Provisional Elder 0 ( 0.0%) 56 (100.0%) 0 ( 0.0%)
Elder in Full Connection 908 (100.0%) 0 ( 0.0%) 0 ( 0.0%)
Student Local Pastor 0 ( 0.0%) 35 (100.0%) 0 ( 0.0%)
Local Pastor (but not a student local pastor) 345 (100.0%) 0 ( 0.0%) 0 ( 0.0%)
Provisional Deacon 0 ( 0.0%) 12 (100.0%) 0 ( 0.0%)
Deacon in Full Connection 86 (100.0%) 0 ( 0.0%) 0 ( 0.0%)
Associate Member 0 ( 0.0%) 0 ( 0.0%) 17 (100.0%)
Other (specify) 0 ( 0.0%) 0 ( 0.0%) 10 (100.0%)
----------------------------------------------- -------------- -------------- ------------- -------------# Create a simplified conference variable #
#-----------------------------------------#
n.attrs <- n.attrs %>%
mutate(
conference_binary = case_when(
ax11 == "North Carolina Annual Conference" ~ "NC Annual",
ax11 == "Western North Carolina Annual Conference" ~ "Western NC Annual",
TRUE ~ NA_character_
),
conference_binary = factor(
conference_binary,
levels = c("NC Annual", "Western NC Annual")
)
)
# Check it
ctable(x = n.attrs$ax11,
y = n.attrs$conference_binary,
prop = "r", #r = row proportions (use "c" or "t" for column/total)
useNA = "no",
totals = FALSE)Cross-Tabulation, Row Proportions
ax11 * conference_binary
Data Frame: n.attrs
----------------------------------------------- ------------------- -------------- -------------------
conference_binary NC Annual Western NC Annual
ax11
North Carolina Annual Conference 580 (100.0%) 0 ( 0.0%)
Western North Carolina Annual Conference 0 ( 0.0%) 847 (100.0%)
A UMC annual conference not in North Carolina 0 ( 0.0%) 0 ( 0.0%)
Not serving in any UMC annual conference 0 ( 0.0%) 0 ( 0.0%)
----------------------------------------------- ------------------- -------------- -------------------# The outcomes... #
#-----------------#
n.attrs <- n.attrs %>%
mutate(
flourisher = factor(
flourisher,
levels = c("Not flourishing", "Flourishing")
)
)
# check
table(n.attrs$flourisher, useNA="ifany")
Not flourishing Flourishing <NA>
474 907 88 Let’s chat about the sample sizes and theoretical implications of the different grouping strategies in our next check-in. Also consider retirement status downstream - potentially excluding pastors who are flagged as retired in the data (perhaps the analysis sample is only n.attrs$ax1==0).
Build the Network
Make the iGraph object
#|include: false
#|echo: false
# First, make sure you don't have any missing fids or uids in the file...
e.attrs <- e.attrs %>%
filter(
!is.na(uid), !is.na(fid), # no missing
uid != 0, # drop the phantom 0
fid != 0
) %>%
distinct(uid, fid, .keep_all = TRUE) # also dedupe
# remove self-nominations
e.attrs <- e.attrs %>% filter(uid != fid)
# make vectors of egos and alters, dropping any missing IDs
egos <- e.attrs |> pull(uid)
alters <- e.attrs |> pull(fid)
nodes <- c(egos, alters) %>%
unique() %>%
tibble(uid = .) %>% # if I don't use old school piping I get an error.
bind_rows(n.attrs %>% select(uid)) %>%
distinct(uid) %>%
filter(!is.na(uid))
# now bring in the attributes
nodes <- nodes |>
left_join(n.attrs, by = "uid")
# build the graph
g21 <- graph_from_data_frame(
d = e.attrs,
directed = TRUE,
vertices = nodes
) |>
igraph::simplify(
remove.multiple = TRUE,
remove.loops = TRUE
)Make some net measures relevant to analysis
Include some generic network measures that allow us to explore structural dependencies in the data (degree, betweenness, etc). This is more or less a holding place for code I expect to build out for now.
Create Relational Support Measures
Make hierarchical measures
We want to use the nodal measures I just made to operationalize who is in what position within the broader status hierarchy for the church.
This section is frought with concerns about analytical decisions I’m making. Based on exchanges with RJP, I may be overstating the extent to which these roles are ordered in ways that correspond to clear social deference or status distinctions. To this really needs to be discussed further. I approach this in three ways:
- A 3-category hierarchy to accompany clrgystatus_3cat.
- A 6-category hierarchy to accompany clrgystatus_6cat.
- A 2-category hierarchy to accompany trainee_flag.
The following code chunk labels and recodes items in the edge attributes file that relate to key types of and quality of support provided by each tie type.
The flags I created for support frequency are not mutually exclusive. They allow us to explore different versions / degrees of “treatment” in terms of support frequency before we make a selection regarding a clear cutoff for any analyses in our project that include them.
#|include: false
#|echo: false
# Shared frequency labels for emotional, spiritual, and vocational support
freq_labels <- c(
"Almost every day",
"More than once a week",
"About once a week",
"More than once a month",
"About once a month",
"More than once a year",
"About once a year",
"Less than once a year",
"Never"
)
# Recode edges: factorize, create flags, and bins
e.attrs <- e.attrs %>%
# 1. preserve raw codes
mutate(
edge_emosupp_code = edge_emosupp,
edge_spirsupp_code = edge_spirsupp,
edge_vocsupp_code = edge_vocsupp
) %>%
# 2. turn into labeled factors
mutate(
edge_emosupp = factor(edge_emosupp_code,
levels = 1:9,
labels = freq_labels),
edge_spirsupp = factor(edge_spirsupp_code,
levels = 1:9,
labels = freq_labels),
edge_vocsupp = factor(edge_vocsupp_code,
levels = 1:9,
labels = freq_labels),
edge_closeness = factor(edge_closeness,
levels = 1:5,
labels = c(
"Extremely close",
"Very close",
"Moderately close",
"Somewhat close",
"Not close at all"
))
) %>%
# 3. flags for >= monthly (<=5) and >= weekly (<=3)
mutate(
emosupp_monthly = if_else(edge_emosupp_code <= 5, 1, 0),
emosupp_weekly = if_else(edge_emosupp_code <= 3, 1, 0),
spirsupp_monthly = if_else(edge_spirsupp_code <= 5, 1, 0),
spirsupp_weekly = if_else(edge_spirsupp_code <= 3, 1, 0),
vocsupp_monthly = if_else(edge_vocsupp_code <= 5, 1, 0),
vocsupp_weekly = if_else(edge_vocsupp_code <= 3, 1, 0)
) %>%
# 4. binned categories
mutate(
emosupp_bin = factor(case_when(
edge_emosupp_code <= 3 ~ "Weekly or more",
edge_emosupp_code <= 5 ~ "Monthly or more",
edge_emosupp_code <= 8 ~ "Less than monthly",
edge_emosupp_code == 9 ~ "Never",
TRUE ~ NA_character_
), levels = c("Weekly or more", "Monthly or more", "Less than monthly", "Never")),
spirsupp_bin = factor(case_when(
edge_spirsupp_code <= 3 ~ "Weekly or more",
edge_spirsupp_code <= 5 ~ "Monthly or more",
edge_spirsupp_code <= 8 ~ "Less than monthly",
edge_spirsupp_code == 9 ~ "Never",
TRUE ~ NA_character_
), levels = c("Weekly or more", "Monthly or more", "Less than monthly", "Never")),
vocsupp_bin = factor(case_when(
edge_vocsupp_code <= 3 ~ "Weekly or more",
edge_vocsupp_code <= 5 ~ "Monthly or more",
edge_vocsupp_code <= 8 ~ "Less than monthly",
edge_vocsupp_code == 9 ~ "Never",
TRUE ~ NA_character_
), levels = c("Weekly or more", "Monthly or more", "Less than monthly", "Never"))
)
# attach the “code” and factor versions
E(g21)$edge_emosupp_code <- e.attrs$edge_emosupp_code
E(g21)$edge_spirsupp_code <- e.attrs$edge_spirsupp_code
E(g21)$edge_vocsupp_code <- e.attrs$edge_vocsupp_code
E(g21)$edge_emosupp <- e.attrs$edge_emosupp
E(g21)$edge_spirsupp <- e.attrs$edge_spirsupp
E(g21)$edge_vocsupp <- e.attrs$edge_vocsupp
E(g21)$edge_closeness <- e.attrs$edge_closeness
# attach the binary flags
E(g21)$emosupp_monthly <- e.attrs$emosupp_monthly
E(g21)$emosupp_weekly <- e.attrs$emosupp_weekly
E(g21)$spirsupp_monthly <- e.attrs$spirsupp_monthly
E(g21)$spirsupp_weekly <- e.attrs$spirsupp_weekly
E(g21)$vocsupp_monthly <- e.attrs$vocsupp_monthly
E(g21)$vocsupp_weekly <- e.attrs$vocsupp_weekly
# attach the coarser bins
E(g21)$emosupp_bin <- e.attrs$emosupp_bin
E(g21)$spirsupp_bin <- e.attrs$spirsupp_bin
E(g21)$vocsupp_bin <- e.attrs$vocsupp_bin
# run a quick check
#table(E(g21)$edge_emosupp, useNA="ifany")
#table(E(g21)$emosupp_monthly)Note: I need to come back and create a set of overarching variables for support frequency that are support-type agnostic.
Note: I also need to come back and remove retired people apparently, oops (5/5/25)
Peep Outcomes
I need to talk to RJP about what these outcomes actually mean. Flourishing is intuitive and I’m familiar with the scale. The other outcomes seem to be reductions that I haven’t been able to map to the codebook at a quick glance.
#|include: false
#|echo: false
summary(V(g21)$codi_mean, na.rm = T) #if what I think it is, lower better Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.000 0.800 1.200 1.261 1.600 3.000 527 summary(V(g21)$min_satis_mean, na.rm = T) #higher probably better Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.000 2.000 2.250 2.244 2.583 3.000 634 #coerce flourisher back into a factor (idk why it broke)
V(g21)$flourisher <- factor(
V(g21)$flourisher,
levels = c("Not flourishing", "Flourishing")
)
summary(V(g21)$flourisher, na.rm=T) # 1 good, 0 badNot flourishing Flourishing NA's
474 907 310 Models :)
#|echo: false
# 1) pull out vertices & edges as tibbles
verts <- as_data_frame(g21, what = "vertices") %>%
# assume you already have a 0/1 numeric flag on the graph
mutate(flourisher_flag = as.integer(flourisher == "Flourishing"))
edges <- as_data_frame(g21, what = "edges")
# 2) compute node‐level exposures: for each node, incoming and/or outgoing
node_exposures <- edges %>%
# you can choose in, out, or both; here we do ALL incident edges:
pivot_longer(c(from, to), names_to="role", values_to="uid") %>%
group_by(uid) %>%
summarize(
# proportion of ties that are at least monthly emotional support
p_emosupp_monthly = mean(emosupp_monthly, na.rm = TRUE),
# proportion at least monthly spiritual support
p_spirsupp_monthly = mean(spirsupp_monthly, na.rm = TRUE),
# proportion at least monthly vocational support
p_vocsupp_monthly = mean(vocsupp_monthly, na.rm = TRUE),
# average closeness (1–5 scale)
mean_closeness = mean(as.integer(edge_closeness), na.rm = TRUE),
.groups = "drop"
)
# 3) bind exposures back onto verts
model_df <- verts %>%
rename(uid = name) %>% # igraph uses `name` for the vertex ID
left_join(node_exposures, by = "uid") %>%
# bring in any other nodal attributes you want, e.g.
mutate(
#age = as.numeric(age),
female = factor(female)#,
#clrgystatus_3cat = clrgystatus_3cat # factor
)
fit0a <- glm(
flourisher_flag ~ trainee_flag,
data = model_df,
family = binomial(link = "logit")
)
fit0b <- glm(
flourisher_flag ~ clrgystatus_3cat,
data = model_df,
family = binomial(link = "logit")
)
fit0c <- glm(
flourisher_flag ~ clrgystatus_6cat,
data = model_df,
family = binomial(link = "logit")
)
broom::tidy(fit0a, conf.int=TRUE, exponentiate = TRUE) %>%
select(term, estimate, p.value, std.error, conf.low, conf.high)# A tibble: 2 × 6
term estimate p.value std.error conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 1.91 1.23e-27 0.0594 1.70 2.15
2 trainee_flag 0.902 6.35e- 1 0.218 0.592 1.39AIC(fit0a, fit0b, fit0c) df AIC
fit0a 2 1750.492
fit0b 3 1730.049
fit0c 6 1734.515BIC(fit0a, fit0b, fit0c) df BIC
fit0a 2 1760.915
fit0b 3 1745.684
fit0c 6 1765.785#looks like the best variance-minimizing choice is 3cat. We'll proceed with the 3cat version for now. We should use ANOVA to highlight that the variance is meaningfully explained by differences between the 3cat groups (moreso than trainee status or other stuff).
#check this one
fit0d <- glm(
flourisher_flag ~ trainee_flag * clrgystatus_3cat,
data = model_df,
family = binomial(link = "logit")
)
BIC(fit0a, fit0b, fit0c, fit0d) df BIC
fit0a 2 1760.915
fit0b 3 1745.684
fit0c 6 1765.785
fit0d 6 1765.785#same conclusion...
# Now bring in some additive predictors of interest:
fit1a <- glm(
flourisher_flag ~ clrgystatus_3cat +
p_emosupp_monthly +
p_spirsupp_monthly +
p_vocsupp_monthly,
data = model_df,
family = binomial(link = "logit")
)
broom::tidy(fit1a, conf.int=TRUE, exponentiate = TRUE) %>%
select(term, estimate, p.value, std.error, conf.low, conf.high)# A tibble: 6 × 6
term estimate p.value std.error conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 1.90 0.0261 0.287 1.09 3.38
2 clrgystatus_3catElder 0.762 0.309 0.267 0.445 1.27
3 clrgystatus_3catLocal Pastor 1.58 0.120 0.295 0.878 2.80
4 p_emosupp_monthly 1.29 0.414 0.308 0.703 2.36
5 p_spirsupp_monthly 0.478 0.0130 0.297 0.265 0.850
6 p_vocsupp_monthly 1.97 0.0152 0.280 1.14 3.43 #hmm - looks like type of support matters a lot, we have much more to do here...but spiritual support's RR is < 1, ans vocational support's RR is almost 2.
fit1b <- glm(
flourisher_flag ~ clrgystatus_3cat*p_emosupp_monthly,
data = model_df,
family = binomial(link = "logit")
)
broom::tidy(fit1b, conf.int=TRUE, exponentiate = TRUE) %>%
select(term, estimate, p.value, std.error, conf.low, conf.high)# A tibble: 6 × 6
term estimate p.value std.error conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 2.75 0.0483 0.513 1.05 8.04
2 clrgystatus_3catElder 0.579 0.307 0.535 0.191 1.59
3 clrgystatus_3catLocal Pastor 0.918 0.880 0.570 0.285 2.72
4 p_emosupp_monthly 0.692 0.624 0.752 0.152 2.98
5 clrgystatus_3catElder:p_emosupp… 1.54 0.586 0.788 0.331 7.47
6 clrgystatus_3catLocal Pastor:p_… 2.44 0.293 0.849 0.467 13.3 fit1c <- glm(
flourisher_flag ~ clrgystatus_3cat*p_spirsupp_monthly,
data = model_df,
family = binomial(link = "logit")
)
broom::tidy(fit1c, conf.int=TRUE, exponentiate = TRUE) %>%
select(term, estimate, p.value, std.error, conf.low, conf.high)# A tibble: 6 × 6
term estimate p.value std.error conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 5.36 0.00754 0.628 1.69 20.4
2 clrgystatus_3catElder 0.345 0.0980 0.644 0.0882 1.13
3 clrgystatus_3catLocal Pastor 0.535 0.351 0.671 0.131 1.87
4 p_spirsupp_monthly 0.249 0.110 0.869 0.0416 1.29
5 clrgystatus_3catElder:p_spirsup… 3.22 0.193 0.898 0.583 20.3
6 clrgystatus_3catLocal Pastor:p_… 5.39 0.0766 0.952 0.877 37.5 # local pastors may benefit tremendously from spiritual support tie interactions
# I need to confirm what tie is here to make sure this isn't giving (I think it is giving; local pastors who send a lot of spiritual support appear to be flourishing).
fit1d <- glm(
flourisher_flag ~ clrgystatus_3cat*p_vocsupp_monthly,
data = model_df,
family = binomial(link = "logit")
)
broom::tidy(fit1d, conf.int=TRUE, exponentiate = TRUE) %>%
select(term, estimate, p.value, std.error, conf.low, conf.high)# A tibble: 6 × 6
term estimate p.value std.error conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 2.43 0.108 0.553 0.858 7.75
2 clrgystatus_3catElder 0.592 0.360 0.573 0.180 1.75
3 clrgystatus_3catLocal Pastor 0.909 0.874 0.601 0.263 2.85
4 p_vocsupp_monthly 0.814 0.787 0.760 0.174 3.55
5 clrgystatus_3catElder:p_vocsupp… 1.60 0.552 0.794 0.343 7.97
6 clrgystatus_3catLocal Pastor:p_… 2.86 0.222 0.860 0.538 16.1 fit2 <- glm(
flourisher_flag ~ clrgystatus_3cat +
p_emosupp_monthly *
p_spirsupp_monthly *
p_vocsupp_monthly,
data = model_df,
family = binomial(link = "logit")
)
broom::tidy(fit2, conf.int=TRUE, exponentiate = TRUE) %>%
select(term, estimate, p.value, std.error, conf.low, conf.high)# A tibble: 10 × 6
term estimate p.value std.error conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 1.46 0.240 0.322 0.783 2.78
2 clrgystatus_3catElder 0.746 0.276 0.269 0.433 1.25
3 clrgystatus_3catLocal Pastor 1.59 0.118 0.298 0.878 2.84
4 p_emosupp_monthly 2.45 0.142 0.610 0.772 8.55
5 p_spirsupp_monthly 0.417 0.193 0.671 0.106 1.52
6 p_vocsupp_monthly 7.33 0.00265 0.663 2.14 28.9
7 p_emosupp_monthly:p_spirsupp_m… 1.44 0.745 1.13 0.161 13.6
8 p_emosupp_monthly:p_vocsupp_mo… 0.110 0.0330 1.04 0.0141 0.824
9 p_spirsupp_monthly:p_vocsupp_m… 0.439 0.473 1.15 0.0462 4.24
10 p_emosupp_monthly:p_spirsupp_m… 2.67 0.519 1.52 0.132 53.0 #People who seek spiritual support monthly AND vocsupp are a messUp next, I need to rebuild the relational / hierarchical bits so I can speak to the “who’s helping whom” question of it all…
…then I need to do some more EDA and think through what the output means.