ct_selection_untitled-B77F94F7.R

gg1 <- gss  %>% 
  filter(year <= 2018) %>% 
  mutate(happy = frcode(happy == 1 ~ "Very", 
                        happy == 2 ~ "Pretty", 
                        happy == 3 ~ "Not too")) %>% 
  group_by(year) %>% 
  ct(happy, wt = wtssall, show_na = FALSE)

gg2 <- gss  %>% 
  filter(year > 2018) %>% 
  mutate(happy = frcode(happy == 1 ~ "Very", 
                        happy == 2 ~ "Pretty", 
                        happy == 3 ~ "Not too")) %>% 
  group_by(year) %>% 
  ct(happy, wt = wtssnrps, show_na = FALSE)

both <- bind_rows(gg1, gg2)

ends <- both %>% filter(year == 1972 | year == 2024)



loess_ends <- both %>%
  group_by(happy) %>%
  arrange(year) %>%
  summarise(
    y_1972 = predict(loess(pct ~ year, span = 0.75), newdata = data.frame(year = 1972)),
    y_2024 = predict(loess(pct ~ year, span = 0.75), newdata = data.frame(year = 2024))
  ) %>%
  pivot_longer(cols = c(y_1972, y_2024), names_to = "year", values_to = "pct") %>%
  mutate(year = as.numeric(str_remove(year, "y_")))

ends <- loess_ends %>%
  mutate(
    nx = ifelse(year == 1972, -1.5, 1.5),
    ny = case_when(
      year == 1972 & happy == "Very"    ~  0.00,
      year == 1972 & happy == "Pretty"  ~ -0.00,
      year == 1972 & happy == "Not too" ~ -0.00,
      TRUE ~ 0
    )
  )

# Then use loess_ends instead of ends in geom_label

both %>% 
  ggplot(., aes(x = year, y = pct, color = happy, group = happy, color = happy)) +
  geom_point(stroke = .5, shape = 21, alpha = .7) +
  geom_labelsmooth(aes(label = happy), method = "loess", formula = y ~ x, family = "bold", 
                   linewidth = 1, text_smoothing = 30, size = 6, boxlinewidth = 0.3) +
  geom_label(data = ends, aes(label = scales::percent(pct, accuracy = 1),
                              nudge_x = nx, nudge_y = ny),  # won't work in aes
             family = "bold", size = 3.5, show.legend = FALSE,
             nudge_x = ends$nx, nudge_y = ends$ny,
             label.size = 0.2) +
  theme_rb() +
  scale_color_manual(values = c("Not too" = "#2166ac", "Pretty" = "#1b7837", "Very" = "#d6191b")) +
  scale_y_continuous(labels = percent, limits = c(0, .6)) +
  labs(x = "", y = "", 
       title = "Taken all together, how would you say things are these days--\nwould you say that you are very happy, pretty happy, or not too happy?", 
       caption = "@ryanburge + @religiondata | Data: General Social Survey, 1972-2024"
  )

save("happy_three_lines_gss24.png")



gg1 <- gss  %>% 
  filter(year <= 2018) %>% 
  mutate(hap = case_when(happy == 3 ~ 1, 
                         happy <= 2 ~ 0)) %>% 
  mutate(trad2 = frcode(reltrad != 7 ~ "Religious", 
                        reltrad == 7 ~ "Non-Religious")) %>% 
  group_by(year, trad2) %>% 
  mean_ci(hap, wt = wtssall, ci = .84)

gg2 <- gss  %>% 
  filter(year > 2018) %>% 
  mutate(hap = case_when(happy == 3 ~ 1, 
                         happy <= 2 ~ 0)) %>% 
  mutate(trad2 = frcode(reltrad != 7 ~ "Religious", 
                        reltrad == 7 ~ "Non-Religious")) %>% 
  group_by(year, trad2) %>% 
  mean_ci(hap, wt = wtssnrps, ci = .84)

both <- bind_rows(gg1, gg2) %>% filter(trad2 != "NA")



loess_ends <- both %>%
  group_by(trad2) %>%
  arrange(year) %>%
  summarise(
    y_start = predict(loess(mean ~ year, span = 0.75), newdata = data.frame(year = min(both$year))),
    y_end   = predict(loess(mean ~ year, span = 0.75), newdata = data.frame(year = max(both$year)))
  ) %>%
  pivot_longer(cols = c(y_start, y_end), names_to = "side", values_to = "mean") %>%
  mutate(
    year = ifelse(side == "y_start", min(both$year), max(both$year)),
    nx   = ifelse(side == "y_start", -1.5, 1.5),
    ny   = case_when(
      side == "y_end" & trad2 == "Non-Religious" ~ -0.02,
      TRUE ~ 0
    )
  )

both %>%
  ggplot(aes(x = year, y = mean, color = trad2, group = trad2)) +
  geom_labelsmooth(aes(label = trad2), method = "loess", formula = y ~ x, family = "bold", 
                   linewidth = 1, text_smoothing = 30, size = 6, boxlinewidth = 0.3) +
  geom_label(data = loess_ends,
             aes(label = scales::percent(mean, accuracy = 1)),
             nudge_x = loess_ends$nx,
             nudge_y = loess_ends$ny,
             family = "bold", size = 3.5, show.legend = FALSE,
             label.size = 0.2) +
  scale_color_calc() + 
  scale_y_continuous(labels = percent, limits = c(0, .27)) +
  labs(x = "", y = "",
       title = "Share Reporting 'Not Too Happy' by Religious Affiliation",
       caption = "@ryanburge + @religiondata | Data: General Social Survey") +
  theme_rb(legend = FALSE)

save("not_happy_gss2lines.png")



gss_cohort <- gss %>%
  mutate(
    birthyr = year - age,
    gen = frcode(
      birthyr <= 1927 ~ "Greatest",
      birthyr >= 1928 & birthyr <= 1945 ~ "Silent",
      birthyr >= 1946 & birthyr <= 1964 ~ "Boomer",
      birthyr >= 1965 & birthyr <= 1980 ~ "Gen X",
      birthyr >= 1981 & birthyr <= 1996 ~ "Millennial",
      birthyr >= 1997 ~ "Gen Z"
    ),
    lifestage = frcode(
      age >= 18 & age <= 29 ~ "18–29",
      age >= 30 & age <= 44 ~ "30–44",
      age >= 45 & age <= 59 ~ "45–59",
      age >= 60 & age <= 74 ~ "60–74",
      age >= 75             ~ "75+"
    ),
    hap = case_when(happy == 3 ~ 1, happy <= 2 ~ 0),
    wt = ifelse(year <= 2018, wtssall, wtssnrps)
  ) %>%
  filter(!is.na(gen), !is.na(lifestage), !is.na(hap))

gg <- gss_cohort %>%
  group_by(gen, lifestage) %>%
  mean_ci(hap, wt = wt, ci = .84)

gg %>%
  mutate(sparse = n_eff < 300) %>%
  ggplot(aes(x = lifestage, y = mean, color = gen, group = gen)) +
  geom_line(linewidth = 1) +
  geom_point(aes(shape = sparse), size = 3) +
  scale_shape_manual(values = c("FALSE" = 19, "TRUE" = 1), guide = "none") +
  scale_y_continuous(labels = percent, limits = c(0, .32)) +
  scale_color_manual(values = c(
    "Greatest"   = "#e41a1c",
    "Silent"     = "#377eb8",
    "Boomer"     = "#ff7f00",
    "Gen X"      = "#4daf4a",
    "Millennial" = "#984ea3",
    "Gen Z"      = "#000000"
  )) +
  guides(color = guide_legend(nrow = 1)) +
  labs(
    x = "Life Stage", y = "",
    title = "Share Reporting 'Not Too Happy' by Generation and Life Stage",
    subtitle = "Each line traces one generation across age groups; open circles indicate sparse cells (N < 300)",
    caption = "@ryanburge + @religiondata | Data: General Social Survey"
  ) +
  theme_rb(legend = TRUE) +
  theme(legend.position = "bottom")
save("not_happy_age_life_stage.png")





gg1 <- gss %>%
  filter(year <= 2018) %>%
  mutate(
    hap = case_when(happy == 3 ~ 1, happy <= 2 ~ 0),
    attend_grp = frcode(
      attend <= 1 ~ "Never/Seldom",
      attend >= 6 ~ "Weekly+"
    )
  ) %>%
  filter(!is.na(attend_grp)) %>%
  group_by(year, attend_grp) %>%
  mean_ci(hap, wt = wtssall, ci = .84)

gg2 <- gss %>%
  filter(year > 2018) %>%
  mutate(
    hap = case_when(happy == 3 ~ 1, happy <= 2 ~ 0),
    attend_grp = frcode(
      attend <= 1 ~ "Never/Seldom",
      attend >= 6 ~ "Weekly+"
    )
  ) %>%
  filter(!is.na(attend_grp)) %>%
  group_by(year, attend_grp) %>%
  mean_ci(hap, wt = wtssnrps, ci = .84)

both <- bind_rows(gg1, gg2)

loess_ends <- both %>%
  group_by(attend_grp) %>%
  arrange(year) %>%
  summarise(
    y_start = predict(loess(mean ~ year, span = 0.75), newdata = data.frame(year = min(both$year))),
    y_end   = predict(loess(mean ~ year, span = 0.75), newdata = data.frame(year = max(both$year)))
  ) %>%
  pivot_longer(cols = c(y_start, y_end), names_to = "side", values_to = "mean") %>%
  mutate(
    year = ifelse(side == "y_start", min(both$year), max(both$year)),
    nx   = ifelse(side == "y_start", -1.5, 1.5),
    ny   = 0
  )

both %>%
  ggplot(aes(x = year, y = mean, color = attend_grp, group = attend_grp)) +
  geom_labelsmooth(aes(label = attend_grp), method = "loess", formula = y ~ x, family = "bold",
                   linewidth = 1, text_smoothing = 30, size = 6, boxlinewidth = 0.3) +
  geom_label(data = loess_ends,
             aes(label = scales::percent(mean, accuracy = 1)),
             nudge_x = loess_ends$nx,
             nudge_y = loess_ends$ny,
             family = "bold", size = 3.5, show.legend = FALSE,
             label.size = 0.2) +
  scale_color_manual(values = c("Never/Seldom" = "#2166ac", "Weekly+" = "#d6191b")) +
  scale_y_continuous(labels = percent, limits = c(0, .27)) +
  labs(x = "", y = "",
       title = "Share Reporting 'Not Too Happy' by Church Attendance",
       caption = "@ryanburge + @religiondata | Data: General Social Survey") +
  theme_rb(legend = FALSE)
save("not_happy_attend_lines.png")



model_df <- gss %>%
  filter(year >= 2021) %>%
  mutate(
    hap = case_when(happy == 3 ~ 1, happy <= 2 ~ 0),
    female = case_when(sex == 2 ~ 1, sex == 1 ~ 0),
    educ_cat = frcode(
      educ < 12  ~ "No HS",
      educ == 12 ~ "HS Diploma",
      educ <= 15 ~ "Some College",
      educ == 16 ~ "BA",
      educ > 16  ~ "Grad Degree"
    ),
    inc_cat = frcode(
      realinc < 20000                        ~ "Low",
      realinc >= 20000 & realinc < 60000     ~ "Middle Income",
      realinc >= 60000                       ~ "High Income"
    ),
    relig_cat = frcode(
      reltrad == 1 ~ "Evangelical",
      reltrad == 2 ~ "Mainline",
      reltrad == 3 ~ "Black Prot",
      reltrad == 5 ~ "Catholic",
      reltrad == 7 ~ "None",
      !is.na(reltrad) ~ "Other"
    ),
    attend_cat = frcode(
      attend <= 1  ~ "Never/Seldom",
      attend %in% 2:5 ~ "Sometimes",
      attend >= 6  ~ "Weekly+"
    ),
    pid3 = frcode(
      partyid %in% 0:2 ~ "Democrat",
      partyid == 3     ~ "Independent",
      partyid %in% 4:6 ~ "Republican"
    ),
    wt = wtssnrps
  ) %>%
  filter(!is.na(hap))

model <- glm(hap ~ female + age + educ_cat + inc_cat + relig_cat + attend_cat + pid3,
             data = model_df,
             weights = wt,
             family = binomial(link = "logit"))

summary(model)


library(tidytext)

tidy(model, conf.int = TRUE, conf.level = 0.84) %>%
  filter(term != "(Intercept)") %>%
  mutate(term = case_when(
    term == "female"                ~ "Female",
    term == "age"                   ~ "Age",
    term == "educ_catHS Diploma"    ~ "HS Diploma",
    term == "educ_catSome College"  ~ "Some College",
    term == "educ_catBA"            ~ "BA",
    term == "educ_catGrad Degree"   ~ "Grad Degree",
    term == "inc_catMiddle Income"  ~ "Middle Income",
    term == "inc_catHigh Income"    ~ "High Income",
    term == "attend_catSometimes"   ~ "Sometimes",
    term == "attend_catWeekly+"     ~ "Weekly+",
    term == "relig_catMainline"     ~ "Mainline Prot",
    term == "relig_catBlack Prot"   ~ "Black Prot",
    term == "relig_catCatholic"     ~ "Catholic",
    term == "relig_catNone"         ~ "None",
    term == "relig_catOther"        ~ "Other Religion",
    term == "pid3Independent"       ~ "Independent",
    term == "pid3Republican"        ~ "Republican"
  )) %>%
  mutate(
    group = case_when(
      term %in% c("Female", "Age")                                    ~ "Demographics",
      term %in% c("HS Diploma", "Some College", "BA", "Grad Degree")  ~ "Education (ref: No HS)",
      term %in% c("Middle Income", "High Income")                     ~ "Income (ref: Low)",
      term %in% c("Sometimes", "Weekly+")                             ~ "Attendance (ref: Never/Seldom)",
      term %in% c("Mainline Prot", "Black Prot", "Catholic",
                  "None", "Other Religion")                           ~ "Religion (ref: Evangelical)",
      term %in% c("Independent", "Republican")                        ~ "Party ID (ref: Democrat)"
    ),
    group = factor(group, levels = c("Demographics",
                                     "Education (ref: No HS)",
                                     "Income (ref: Low)",
                                     "Attendance (ref: Never/Seldom)",
                                     "Religion (ref: Evangelical)",
                                     "Party ID (ref: Democrat)")),
    significant = sign(conf.low) == sign(conf.high)
  ) %>%
  ggplot(aes(x = estimate, y = reorder_within(term, estimate, group), color = significant)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.2, linewidth = 0.7) +
  geom_point(size = 3) +
  facet_wrap(~ group, scales = "free_y", ncol = 1) +
  scale_y_reordered() +
  scale_color_manual(values = c("TRUE" = "#d6191b", "FALSE" = "gray60")) +
  labs(
    x = "Log-Odds (positive = more likely 'Not Too Happy')", y = "",
    title = "Predictors of Reporting 'Not Too Happy'",
    caption = "@ryanburge + @religiondata | Data: General Social Survey, 2021-2024"
  ) +
  theme_rb(legend = FALSE)
save("not_happy_coefplot.png", ht = 10)



# Option 1: Simple predicted probabilities at each attendance level
library(marginaleffects)

plot_predictions(model, condition = "attend_cat") +
  scale_y_continuous(labels = percent, limits = c(0, .40)) +
  scale_x_discrete(limits = c("Never/Seldom", "Sometimes", "Weekly+")) +
  labs(
    x = "", y = "",
    title = "Church Attendance and Unhappiness",
    subtitle = "Predicted probability of reporting 'Not Too Happy'\nholding all other predictors at their means",
    caption = "@ryanburge + @religiondata | Data: General Social Survey, 2021-2024"
  ) +
  theme_rb(legend = FALSE)



library(marginaleffects)

predictions(model, 
            newdata = datagrid(attend_cat = c("Never/Seldom", "Sometimes", "Weekly+"))) %>%
  as_tibble() %>%
  mutate(attend_cat = factor(attend_cat, levels = c("Never/Seldom", "Sometimes", "Weekly+"))) %>%
  ggplot(aes(x = attend_cat, y = estimate)) +
  geom_col(fill = "#2166ac", alpha = 0.8, width = 0.75, color = "black") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.1, linewidth = 0.7) +
  geom_text(aes(y = .03, label = scales::percent(estimate, accuracy = 1)),
            family = "bold", size = 11) +
  scale_y_continuous(labels = percent, limits = c(0, .37)) +
  labs(
    x = "", y = "",
    title = "The More You Attend Church, the Happier You Are",
    subtitle = "Predicted probability of 'Not Too Happy' by attendance\nControlling for income, education, religion, age, gender, and party",
    caption = "@ryanburge + @religiondata | Data: General Social Survey, 2021-2024"
  ) +
  theme_rb(legend = FALSE)
save("not_happy_attend_bars.png", wd = 6)


library(marginaleffects)

inc   <- predictions(model, newdata = datagrid(inc_cat    = c("Low", "High Income"))) %>%
  as_tibble() %>% mutate(predictor = "Income", label = inc_cat)
educ  <- predictions(model, newdata = datagrid(educ_cat   = c("No HS", "Grad Degree"))) %>%
  as_tibble() %>% mutate(predictor = "Education", label = educ_cat)
att   <- predictions(model, newdata = datagrid(attend_cat = c("Never/Seldom", "Weekly+"))) %>%
  as_tibble() %>% mutate(predictor = "Attendance", label = attend_cat)
rel   <- predictions(model, newdata = datagrid(relig_cat  = c("Evangelical", "None"))) %>%
  as_tibble() %>% mutate(predictor = "Religion", label = relig_cat)
pid   <- predictions(model, newdata = datagrid(pid3       = c("Democrat", "Republican"))) %>%
  as_tibble() %>% mutate(predictor = "Party ID", label = pid3)
gen   <- predictions(model, newdata = datagrid(female     = c(0, 1))) %>%
  as_tibble() %>% mutate(predictor = "Gender", label = ifelse(female == 0, "Male", "Female"))



bind_rows(
  predictions(model, newdata = datagrid(inc_cat    = c("Low", "High Income"))) %>%
    as_tibble() %>% mutate(predictor = "Income", label = inc_cat),
  predictions(model, newdata = datagrid(educ_cat   = c("No HS", "Grad Degree"))) %>%
    as_tibble() %>% mutate(predictor = "Education", label = educ_cat),
  predictions(model, newdata = datagrid(attend_cat = c("Never/Seldom", "Weekly+"))) %>%
    as_tibble() %>% mutate(predictor = "Attendance", label = attend_cat),
  predictions(model, newdata = datagrid(relig_cat  = c("Evangelical", "None"))) %>%
    as_tibble() %>% mutate(predictor = "Religion", label = relig_cat),
  predictions(model, newdata = datagrid(pid3       = c("Democrat", "Republican"))) %>%
    as_tibble() %>% mutate(predictor = "Party ID", label = pid3),
  predictions(model, newdata = datagrid(female     = c(0, 1))) %>%
    as_tibble() %>% mutate(predictor = "Gender", label = ifelse(female == 0, "Male", "Female"))
) %>%
  mutate(predictor = factor(predictor, levels = c("Income", "Attendance", "Religion",
                                                  "Education", "Party ID", "Gender"))) %>%
  ggplot(aes(x = estimate, y = predictor)) +
  geom_line(aes(group = predictor), color = "gray70", linewidth = 2) +
  geom_point(size = 5, color = "#2166ac") +
  geom_text(aes(label = paste0(label, "\n", scales::percent(estimate, accuracy = 1))),
            vjust = -0.7, family = "bold", size = 3.5) +
  scale_x_continuous(labels = percent, limits = c(0, .40)) +
  labs(
    x = "Predicted Probability of 'Not Too Happy'", y = "",
    title = "What Predicts Unhappiness? Low vs. High End of Each Factor",
    subtitle = "Predicted probabilities from logistic regression | All other variables held at means",
    caption = "@ryanburge + @religiondata | Data: General Social Survey, 2021-2024"
  ) +
  theme_rb(legend = FALSE)
save("not_happy_dumbbell_all.png")



gg1 <- gss %>%
  filter(year <= 2018) %>%
  mutate(
    hap = case_when(happy == 3 ~ 1, happy <= 2 ~ 0),
ideo = frcode(
  polviews %in% 1:3 ~ "Liberal",
  polviews == 4     ~ "Moderate",
  polviews %in% 5:7 ~ "Conservative"
)
  ) %>%
  filter(!is.na(ideo)) %>%
  group_by(year, ideo) %>%
  mean_ci(hap, wt = wtssall, ci = .84)

gg2 <- gss %>%
  filter(year > 2018) %>%
  mutate(
    hap = case_when(happy == 3 ~ 1, happy <= 2 ~ 0),
    ideo = frcode(
      polviews %in% 1:2 ~ "Liberal",
      polviews == 4     ~ "Moderate",
      polviews %in% 6:7 ~ "Conservative"
    )
  ) %>%
  filter(!is.na(ideo)) %>%
  group_by(year, ideo) %>%
  mean_ci(hap, wt = wtssnrps, ci = .84)

both <- bind_rows(gg1, gg2)

loess_ends <- both %>%
  group_by(ideo) %>%
  arrange(year) %>%
  summarise(
    y_start = predict(loess(mean ~ year, span = 0.75), newdata = data.frame(year = min(both$year))),
    y_end   = predict(loess(mean ~ year, span = 0.75), newdata = data.frame(year = max(both$year)))
  ) %>%
  pivot_longer(cols = c(y_start, y_end), names_to = "side", values_to = "mean") %>%
  mutate(
    year = ifelse(side == "y_start", min(both$year), max(both$year)),
    nx   = ifelse(side == "y_start", -1.5, 1.5),
    ny   = case_when(
      side == "y_end" & ideo == "Liberal"      ~  -0.01,
      side == "y_end" & ideo == "Moderate"     ~ -0.01,
      TRUE ~ 0
    )
  )

both %>%
  ggplot(aes(x = year, y = mean, color = ideo, group = ideo)) +
  geom_labelsmooth(aes(label = ideo), method = "loess", formula = y ~ x, family = "bold",
                   linewidth = 1, text_smoothing = 30, size = 6, boxlinewidth = 0.3) +
  geom_label(data = loess_ends,
             aes(label = scales::percent(mean, accuracy = 1)),
             nudge_x = loess_ends$nx,
             nudge_y = loess_ends$ny,
             family = "bold", size = 3.5, show.legend = FALSE,
             label.size = 0.2) +
  scale_color_manual(values = c(
    "Liberal"      = "#2166ac",
    "Moderate"     = "#4daf4a",
    "Conservative" = "#d6191b"
  )) +
  scale_y_continuous(labels = percent, limits = c(0, .27)) +
  labs(x = "", y = "",
       title = "Share Reporting 'Not Too Happy' by Political Ideology",
       caption = "@ryanburge + @religiondata | Data: General Social Survey") +
  theme_rb(legend = FALSE)
save("not_happy_ideo_lines.png")