ct_selection_untitled-329C0CD0.R

# Never attends = 1; all other valid attend values = 0; keep NA as NA
mk_never <- function(x) ifelse(x == 0, 1L,
                               ifelse(x %in% 1:8, 0L, NA_integer_))

# 1972–2002
pre2004 <- gss %>%
  filter(year <= 2002) %>%
  mutate(never = mk_never(attend)) %>%
  group_by(year) %>%
  mean_ci(never, wt = wtssps)

# 2004–2024
post2004 <- gss %>%
  filter(year >= 2004) %>%
  mutate(never = mk_never(attend)) %>%
  group_by(year) %>%
  mean_ci(never, wt = wtssnrps)

never <- bind_rows(pre2004, post2004) %>% arrange(year)

never %>% 
  ggplot(., aes(x = year, y = mean)) +
  geom_point() +
  geom_line() +
  scale_y_continuous(labels =  percent, limits = c(0, .35)) +
  theme_rb() + 
  labs(x = "", y = "", title = "Share Who Never Attend Religious Services, 1972-2024", caption = "@ryanburge + @religiondata | Data: General Social Survey, 1972-2024")
save("gss_never_attends.png")


never %>% 
  ggplot(aes(x = year, y = mean)) +
  geom_area(fill = "#a2cffe", alpha = 0.6) +                                # soft blue area
  geom_line(color = "#0057b7", linewidth = 1.3) +                           # bold line
  geom_point(color = "#003f88", size = 2) +                                 # darker dots
  scale_y_continuous(
    labels = scales::percent_format(accuracy = 1),
    limits = c(0, 0.35),
    breaks = seq(0, 0.35, 0.05)
  ) +
  scale_x_continuous(
    breaks = seq(1972, 2024, 4),
    expand = expansion(mult = c(0.01, 0.01))
  ) +
  theme_rb() +
  labs(
    title = "Americans Who Never Attend Religious Services, 1972–2024",
    caption = "@ryanburge + @religiondata | Data: GSS 1972–2024",
    x = "", y = ""
  ) +
  theme(
    plot.title = element_text(size = 18, face = "bold"),
    plot.subtitle = element_text(size = 13, color = "gray40"),
    panel.grid.minor = element_blank()
  ) +
  add_text(x = 1974, y = .09, word = "10%", sz = 6) +
  add_text(x = 2000, y = .225, word = "21%", sz = 6) +
  add_text(x = 2022, y = .345, word = "33%", sz = 6) +
  add_text(x = 2022.5, y = .280, word = "29%", sz = 6)



save("gss_never_attends_pretty.png", wd = 9, ht = 6)


# --- Non-attenders who are "Nones" among total non-attenders ---

# 1988–2002: use wtssps
none_pre2004 <- gss %>% 
  filter(year >= 1972, year <= 2002, attend == 0) %>% 
  mutate(none = case_when(reltrad == 7 ~ 1,
                          reltrad <= 6 ~ 0,
                          TRUE ~ NA_real_)) %>% 
  group_by(year) %>% 
  mean_ci(none, wt = wtssps)

# 2004–2024: use wtssnrps
none_post2004 <- gss %>% 
  filter(year >= 2004, attend == 0) %>% 
  mutate(none = case_when(reltrad == 7 ~ 1,
                          reltrad <= 6 ~ 0,
                          TRUE ~ NA_real_)) %>% 
  group_by(year) %>% 
  mean_ci(none, wt = wtssnrps)

# Combine and clean
gg <- bind_rows(none_pre2004, none_post2004) %>% arrange(year) %>% filter(year != 2021)

gg

gg %>% 
  ggplot(aes(x = year, y = mean)) +
  geom_area(fill = "#FFD6A5", alpha = 0.6) +             # soft orange fill
  geom_line(color = "#FF7A00", linewidth = 1.3) +        # bold orange line
  geom_point(color = "#C2410C", size = 2) +              # darker accent points

  scale_x_continuous(
    breaks = seq(1972, 2024, 4),
    expand = expansion(mult = c(0.01, 0.01))
  ) +
  theme_rb() +
  y_pct() + 
  labs(
    title = "Among Americans Who Never Attend, Share with No Religion, 1972–2024",
    caption = "@ryanburge + @religiondata | Data: General Social Survey 1972–2024",
    x = "", y = ""
  ) +
  theme(
    plot.title = element_text(size = 18, face = "bold"),
    panel.grid.minor = element_blank()
  ) +
  add_text(x = 1973.5, y = 0.325, word = "35%", sz = 6) +
  
  add_text(x = 1989.5, y = 0.345, word = "38%", sz = 6) +
  add_text(x = 2004, y = 0.585, word = "56%", sz = 6) +
  add_text(x = 2023, y = 0.62, word = "65%", sz = 6)

save("gss_never_attend_share_none_orange.png", wd = 9, ht = 6)


gd_levels <- c("Don't Believe","No Way to Find Out","Some Higher Power",
               "Believe Sometimes","Believe But Doubts","Believe No Doubts")

# 1988 non-attenders
gd_1988 <- gss %>%
  filter(year == 1988, attend == 0) %>%
  mutate(gd = frcode(god == 1 ~ "Don't Believe",
                     god == 2 ~ "No Way to Find Out",
                     god == 3 ~ "Some Higher Power",
                     god == 4 ~ "Believe Sometimes",
                     god == 5 ~ "Believe But Doubts",
                     god == 6 ~ "Believe No Doubts")) %>%
  group_by(year) %>%
  ct(gd, wt = wtssps, show_na = FALSE) %>%
  mutate(gd = factor(gd, levels = gd_levels))

gd_1998 <- gss %>%
  filter(year == 1998, attend == 0) %>%
  mutate(gd = frcode(god == 1 ~ "Don't Believe",
                     god == 2 ~ "No Way to Find Out",
                     god == 3 ~ "Some Higher Power",
                     god == 4 ~ "Believe Sometimes",
                     god == 5 ~ "Believe But Doubts",
                     god == 6 ~ "Believe No Doubts")) %>%
  group_by(year) %>%
  ct(gd, wt = wtssps, show_na = FALSE) %>%
  mutate(gd = factor(gd, levels = gd_levels))

gd_2008 <- gss %>%
  filter(year == 2008, attend == 0) %>%
  mutate(gd = frcode(god == 1 ~ "Don't Believe",
                     god == 2 ~ "No Way to Find Out",
                     god == 3 ~ "Some Higher Power",
                     god == 4 ~ "Believe Sometimes",
                     god == 5 ~ "Believe But Doubts",
                     god == 6 ~ "Believe No Doubts")) %>%
  group_by(year) %>%
  ct(gd, wt = wtssps, show_na = FALSE) %>%
  mutate(gd = factor(gd, levels = gd_levels))


gd_2018 <- gss %>%
  filter(year == 2018, attend == 0) %>%
  mutate(gd = frcode(god == 1 ~ "Don't Believe",
                     god == 2 ~ "No Way to Find Out",
                     god == 3 ~ "Some Higher Power",
                     god == 4 ~ "Believe Sometimes",
                     god == 5 ~ "Believe But Doubts",
                     god == 6 ~ "Believe No Doubts")) %>%
  group_by(year) %>%
  ct(gd, wt = wtssps, show_na = FALSE) %>%
  mutate(gd = factor(gd, levels = gd_levels))

# 2024 non-attenders
gd_2024 <- gss %>%
  filter(year == 2024, attend == 0) %>%
  mutate(gd = frcode(god == 1 ~ "Don't Believe",
                     god == 2 ~ "No Way to Find Out",
                     god == 3 ~ "Some Higher Power",
                     god == 4 ~ "Believe Sometimes",
                     god == 5 ~ "Believe But Doubts",
                     god == 6 ~ "Believe No Doubts")) %>%
  group_by(year) %>%
  ct(gd, wt = wtssnrps, show_na = FALSE) %>%
  mutate(gd = factor(gd, levels = gd_levels))

gg1 <- bind_rows(gd_1988, gd_1998, gd_2008, gd_2018, gd_2024) %>% arrange(year, gd)
gg1



gg1 %>% 
  mutate(lab = round(pct, 2)) %>% 
  ggplot(., aes(x = 1, y = pct, fill = fct_rev(gd))) +
  geom_col(color = "black") + 
  coord_flip() +
  facet_wrap(~ year, ncol =1, strip.position = "left") +
  theme_rb(legend = TRUE) + 
  scale_fill_manual(values = c(
    "Don't Believe"        = "#2F4858",
    "No Way to Find Out"   = "#33658A",
    "Some Higher Power"    = "#86BBD8",
    "Believe Sometimes"    = "#F6AE2D",
    "Believe But Doubts"   = "#F26419",
    "Believe No Doubts"    = "#8AC926"
  )) +
  theme(legend.position = "bottom") +
  scale_y_continuous(labels = percent) +
  theme(strip.text.y.left = element_text(angle = 0)) +
  guides(fill = guide_legend(reverse=T, nrow = 1)) +
  theme(axis.title.y=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) +
  theme(panel.grid.minor.y=element_blank(), panel.grid.major.y=element_blank()) +
  geom_text(aes(label = ifelse(pct >.08, paste0(lab*100, '%'), '')), position = position_stack(vjust = 0.5), size = 8, family = "font", color = "white") +
  labs(x = "", y = "", title = "Belief in God Among Never Attenders",  caption = "@ryanburge + @religiondata | Data: General Social Survey, 1988-2024") 
save("never_belief_compare.png", wd = 11, ht = 4.5)


spirit_lvls <- c("Very", "Moderately", "Slightly", "Not at all")

spirit_pre2004 <- gss %>%
  filter(year >= 1988, year <= 2002, attend == 0, sprtprsn %in% 1:4) %>%
  mutate(spirit = frcode(
    sprtprsn == 1 ~ "Very",
    sprtprsn == 2 ~ "Moderately",
    sprtprsn == 3 ~ "Slightly",
    sprtprsn == 4 ~ "Not at all"
  )) %>%
  group_by(year) %>%
  ct(spirit, wt = wtssps, show_na = FALSE) %>%
  mutate(spirit = factor(spirit, levels = spirit_lvls))

# 2004–2024
spirit_post2004 <- gss %>%
  filter(year >= 2004, attend == 0, sprtprsn %in% 1:4) %>%
  mutate(spirit = frcode(
    sprtprsn == 1 ~ "Very",
    sprtprsn == 2 ~ "Moderately",
    sprtprsn == 3 ~ "Slightly",
    sprtprsn == 4 ~ "Not at all"
  )) %>%
  group_by(year) %>%
  ct(spirit, wt = wtssnrps, show_na = FALSE) %>%
  mutate(spirit = factor(spirit, levels = spirit_lvls))

# Combine
spirit_by_year <- bind_rows(spirit_pre2004, spirit_post2004) %>% arrange(year, spirit)
spirit_by_year


spirit_by_year %>% 
  ggplot(., aes(x = year, y = pct, color = spirit, group = spirit)) +
  geom_point(stroke = 1, shape = 21, fill = "white") +
  geom_labelsmooth(aes(label = spirit), method = "loess", formula = y ~ x, 
                   family = "font", linewidth = 1, text_smoothing = 30, 
                   size = 6, linewidth = 1, boxlinewidth = 0.3, hjust = .15) +
  theme_rb() + 
  scale_y_continuous(labels = percent, limits = c(0, .4)) +
  scale_color_manual(values = c("Very" = "#E41A1C", 
                                "Moderately" = "#377EB8", 
                                "Slightly" = "#4DAF4A", 
                                "Not at all" = "#FF7F00")) +
  labs(x = "", y = "", 
       title = "To what extent do you consider yourself a spiritual person?\nAmong those who never attend religious services", 
       caption = "@ryanburge + @religiondata\nData: General Social Survey, 1998-2024") 

save("spirit_person_nevers_gss24.png")




gss_sec <- gss %>%
  filter(attend == 0) %>% 
  mutate(
    # 1. Religious affiliation: none = 1, else 0
    none_affil = if_else(reltrad == 7, 1, 0),
    
    # 2. Belief in God: code from most to least belief
    god_score = case_when(
      god == 6 ~ 0,  # Believe, no doubts
      god == 5 ~ 0.2,
      god == 4 ~ 0.4,
      god == 3 ~ 0.6,
      god == 2 ~ 0.8,
      god == 1 ~ 1.0,  # Don’t believe in God
      TRUE ~ NA_real_
    ),
    
    # 3. Spirituality (sprtprsn): not at all = 1, very spiritual = 0
    spirit_score = case_when(
      sprtprsn == 4 ~ 1.0,
      sprtprsn == 3 ~ 0.66,
      sprtprsn == 2 ~ 0.33,
      sprtprsn == 1 ~ 0,
      TRUE ~ NA_real_
    )
  )


gss_sec <- gss_sec %>%
  mutate(
    secularity_index = rowMeans(select(., none_affil, god_score, spirit_score), na.rm = TRUE)
  )


# ensure the index exists only where inputs exist (1988+), then give it a simple name
gss_sec_idx <- gss_sec %>%
  filter(year >= 1988) %>%                 # spirituality starts in 1988
  mutate(
    wt  = if_else(year <= 2002, wtssps, wtssnrps),
    idx = secularity_index                 # <- plain column for mean_ci()
  ) %>%
  filter(!is.na(idx), !is.na(wt))




mean_ci <- function(df, var, wt = NULL, ci = 0.95, dist = c("t", "normal"), na.rm = TRUE) {
  stopifnot(is.numeric(ci), length(ci) == 1, is.finite(ci), ci > 0, ci < 1)
  dist <- match.arg(dist)
  
  var <- rlang::enquo(var)
  wt_missing <- missing(wt)  # Check BEFORE enquo
  wt  <- rlang::enquo(wt)
  
  # Validate: exactly one column selected and numeric
  sel <- tidyselect::eval_select(var, df)
  if (length(sel) != 1L) stop("`var` must select exactly one numeric column.", call. = FALSE)
  if (!is.numeric(df[[sel]])) stop("`var` must refer to a numeric column.", call. = FALSE)
  
  compute_stats <- function(x, w = NULL, ci = 0.95, dist = "t") {
    if (!is.numeric(x)) stop("`var` must be numeric.", call. = FALSE)
    
    if (is.null(w)) {
      if (na.rm) x <- x[!is.na(x)]
      if (!length(x)) {
        return(tibble::tibble(mean = NA_real_, sd = NA_real_, n = 0L, n_eff = 0,
                              se = NA_real_, lower = NA_real_, upper = NA_real_, ci = ci))
      }
      N_raw <- length(x)
      m  <- mean(x)
      sd <- stats::sd(x)
      df_ <- max(N_raw - 1, 1)
      alpha <- 1 - ci
      crit <- if (dist == "t") stats::qt(1 - alpha/2, df = df_) else stats::qnorm(1 - alpha/2)
      se <- sd / sqrt(N_raw)
      return(tibble::tibble(
        mean = m, sd = sd, n = N_raw, n_eff = N_raw,
        se = se, lower = m - crit * se, upper = m + crit * se, ci = ci
      ))
    } else {
      if (!is.numeric(w)) stop("`wt` must be numeric.", call. = FALSE)
      if (length(w) != length(x)) stop("`wt` must be same length as `var`.", call. = FALSE)
      
      keep <- is.finite(w) & (w > 0)
      if (na.rm) keep <- keep & !is.na(x) & !is.na(w)
      x <- x[keep]; w <- w[keep]
      if (!length(x)) {
        return(tibble::tibble(mean = NA_real_, sd = NA_real_, n = 0L, n_eff = 0,
                              se = NA_real_, lower = NA_real_, upper = NA_real_, ci = ci))
      }
      w_sum <- sum(w)
      if (!is.finite(w_sum) || w_sum <= 0) stop("Sum of weights must be positive.", call. = FALSE)
      
      N_raw <- length(x)
      w_norm <- w / w_sum
      m <- sum(w_norm * x)
      
      var_w_pop <- sum(w_norm * (x - m)^2)
      neff <- (w_sum^2) / sum(w^2)
      var_w_unb <- if (neff > 1) var_w_pop * (neff / (neff - 1)) else var_w_pop
      
      sd_w <- sqrt(var_w_unb)
      df_ <- max(neff - 1, 1)
      alpha <- 1 - ci
      crit <- if (dist == "t") stats::qt(1 - alpha/2, df = df_) else stats::qnorm(1 - alpha/2)
      se <- sqrt(var_w_unb / neff)
      
      return(tibble::tibble(
        mean = m, sd = sd_w, n = as.integer(N_raw), n_eff = as.numeric(neff),
        se = se, lower = m - crit * se, upper = m + crit * se, ci = ci
      ))
    }
  }
  
  dplyr::summarise(
    df,
    dplyr::across(
      {{ var }},
      ~ {
        if (wt_missing) {  # Use the flag instead
          compute_stats(.x, w = NULL, ci = ci, dist = dist)
        } else {
          wcol <- dplyr::pick(!!wt)[[1]]
          compute_stats(.x, w = wcol, ci = ci, dist = dist)
        }
      },
      .names = "{.col}"
    ),
    .groups = "drop"
  ) |>
    tidyr::unnest(cols = dplyr::all_of(rlang::as_name(var)))
}



gg1 <- gss_sec_idx %>%
  group_by(year) %>%
  mean_ci(idx, wt = wt) %>% filter(year != 2021) %>% 
  mutate(gen = "Entire Sample")





library(dplyr)
library(ggplot2)

# ---- 1) Build generations from birth year ----
# (GSS has age; derive birthyr = year - age)
gen_levels <- c("Gen Z","Millennial","Gen X","Boomer","Silent")

gss_gen_idx <- gss %>%
  filter(attend == 0, year >= 1988) %>%                    # spirituality starts in 1988
  mutate(
    birthyr = year - age,
    gen = case_when(
      birthyr >= 1997 ~ "Gen Z",
      birthyr >= 1981 & birthyr <= 1996 ~ "Millennial",
      birthyr >= 1965 & birthyr <= 1980 ~ "Gen X",
      birthyr >= 1946 & birthyr <= 1964 ~ "Boomer",
      birthyr <= 1945 ~ "Silent",
      TRUE ~ NA_character_
    ),
    gen = factor(gen, levels = gen_levels),
    
    # ---- 2) Secularity Index components (higher = more secular) ----
    none_affil   = as.integer(reltrad == 7),   # 1 = None, 0 = any religion
    god_score = case_when(
      god == 6 ~ 0.0,  # Believe, no doubts
      god == 5 ~ 0.2,
      god == 4 ~ 0.4,
      god == 3 ~ 0.6,
      god == 2 ~ 0.8,
      god == 1 ~ 1.0,  # Don't believe
      TRUE ~ NA_real_
    ),
    spirit_score = case_when(                # sprtprsn: 1 Very .. 4 Not at all
      sprtprsn == 4 ~ 1.0,
      sprtprsn == 3 ~ 2/3,
      sprtprsn == 2 ~ 1/3,
      sprtprsn == 1 ~ 0.0,
      TRUE ~ NA_real_
    ),
    idx = rowMeans(across(c(none_affil, god_score, spirit_score)), na.rm = TRUE),
    
    # Era-correct weights
    wt = if_else(year <= 2002, wtssps, wtssnrps)
  ) %>%
  filter(!is.na(gen), !is.na(idx), is.finite(wt))

# ---- 3) Year × Generation summary with 84% CI ----
gg2 <- gss_gen_idx %>%
  group_by(year, gen) %>%
  mean_ci(idx, wt = wt, ci = .84) %>%
  ungroup() %>%
  filter(year != 2021)   # optional but recommended

both <- bind_rows(gg1, gg2)

library(ggplot2)
library(dplyr)
library(scales)

# filter out years with too few respondents per generation
both_clean <- both %>%
  filter(n_eff >= 75, !is.na(mean), !is.na(gen), year != 2021) %>% filter(gen != "Silent") %>% 
  mutate(
    gen = factor(gen, levels = c("Entire Sample", "Boomer", "Gen X", "Millennial"))
  ) 

# prettier version
both_clean %>%
  ggplot(aes(x = year, y = mean, color = gen, group = gen)) +
  geom_point(size = 2, shape = 21, stroke = .5) +
  geom_smooth(se = FALSE, linewidth = 1, linetype = "solid", alpha = .3) +
  scale_x_continuous(
    breaks = seq(1988, 2024, 4),
    expand = expansion(mult = c(.01, .01))
  ) +
  scale_color_manual(
    values = c(
      "Gen Z"       = "#FF6F00",
      "Millennial"  = "#FF9E00",
      "Gen X"       = "#00A6ED",
      "Boomer"      = "#7FB800",
      "Silent"      = "#9E0059", 
      "Entire Sample" = "black"
    )
  ) +
  theme_rb(legend = TRUE) +
  theme(
    legend.position = "bottom",
    legend.text = element_text(size = 20),
    plot.title = element_text(size = 18, face = "bold"),
    panel.grid.minor = element_blank()
  ) +
  labs(
    title = "Secularity Index Among Never-Attenders, by Generation (1988–2024)",
    caption = "@ryanburge + @religiondata | Data: General Social Survey, 1988-2024",
    x = "", y = ""
  )

save("gss_secularity_index_by_gen_clean.png", wd = 10, ht = 6)