ct_selection_untitled-38D84864.R

library(rio)
library(janitor)

many <- import("E://data/many.sav") %>% clean_names()

# Recode country variable to proper names
many <- many %>%
  mutate(country = case_when(
    country == 1 ~ "Australia",
    country == 2 ~ "Belgium",
    country == 3 ~ "Brazil",
    country == 4 ~ "Canada",
    country == 5 ~ "Chile",
    country == 6 ~ "China",
    country == 7 ~ "Croatia",
    country == 8 ~ "Denmark",
    country == 9 ~ "France",
    country == 10 ~ "Germany",
    country == 11 ~ "India",
    country == 12 ~ "Ireland",
    country == 13 ~ "Israel",
    country == 14 ~ "Italy",
    country == 15 ~ "Japan",
    country == 16 ~ "Lithuania",
    country == 17 ~ "Morocco",
    country == 18 ~ "Netherlands",
    country == 19 ~ "Romania",
    country == 20 ~ "Singapore",
    country == 21 ~ "Spain",
    country == 22 ~ "Turkey",
    country == 23 ~ "United Kingdom",
    country == 24 ~ "United States",
    TRUE ~ as.character(country)
  ),
  region = case_when(
    country %in% c("Australia", "Japan", "Singapore") ~ "Asia-Pacific",
    country %in% c("Belgium", "Croatia", "Denmark", "France", "Germany", "Ireland", "Italy", "Lithuania", "Netherlands", "Romania", "Spain", "United Kingdom") ~ "Europe",
    country %in% c("Brazil", "Chile") ~ "Latin America",
    country %in% c("China", "India") ~ "Asia",
    country %in% c("Morocco", "Turkey") ~ "Middle East & North Africa",
    country %in% c("Canada", "United States") ~ "North America",
    country == "Israel" ~ "Middle East",
    TRUE ~ "Other"
  ))



gg <- many %>% 
  mutate(rel = frcode(rel_3 == 0 ~ "Atheist", 
                      rel_3 == .5 ~ "Not Religious", 
                      rel_3 == 1 ~ "Religious")) %>% 
  group_by(country) %>% 
  ct(rel, show_na = FALSE)


lvl <- gg %>% 
  filter(rel == "Atheist") %>% 
  group_by(country) %>% 
  summarise(sum = sum(pct)) %>% 
  select(country, sort = sum)

gg <- left_join(gg, lvl)

gg$country <- fct_reorder(gg$country, gg$sort, .desc = TRUE)

distinguishable_colors <- diverging_hcl(3, palette = "Purple-Blue")


gg %>% 
  mutate(lab = round(pct, 2)) %>% 
  ggplot(., aes(x = 1, y = pct, fill = fct_rev(rel))) +
  geom_col(color = "black") + 
  coord_flip() +
  facet_wrap(~ country, ncol =1, strip.position = "left") +
  theme_rb() +
  scale_fill_manual(values = c("#60267F", "#28676F", "#B4B681")) +
  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 >.05 & rel == "Religious", paste0(lab*100, '%'), '')), position = position_stack(vjust = 0.5), size = 4, family = "font", color = "white") +
  geom_text(aes(label = ifelse(pct >.05 & rel == "Not Religious", paste0(lab*100, '%'), '')), position = position_stack(vjust = 0.5), size = 4, family = "font", color = "white") +
  geom_text(aes(label = ifelse(pct >.05 & rel == "Atheist", paste0(lab*100, '%'), '')), position = position_stack(vjust = 0.5), size = 4, family = "font", color = "black") +
  theme(plot.title = element_text(size = 16)) +
  theme(strip.text.y.left = element_text(angle = 0, hjust = 1)) +
  theme(legend.text = element_text(size = 20)) +
  labs(x = "", y = "", title = "Independently of whether you attend religious services or not, would you say you are:", 
       caption = "@ryanburge + @religiondata\nData: Many-Analysts Religion Project, 2019") 
save("many_rel_person.png", wd = 9, ht = 8)


gg1 <- many %>% 
  group_by(country) %>% 
  mean_ci(rel_1) %>% 
  select(country, att = mean)

gg2 <- many %>% 
  group_by(country) %>% 
  mean_ci(rel_2) %>% 
  select(country, pray = mean, n)

scatter <- left_join(gg1, gg2)

reg <- many %>% 
  select(country, region) %>% 
  distinct()

scatter <- left_join(scatter, reg)

pop <- data.frame(
  country = c("Australia", "Belgium", "Brazil", "Canada", "Chile", "China", "Croatia", "Denmark", "France", "Germany", "India", "Ireland", "Israel", "Italy", "Japan", "Lithuania", "Morocco", "Netherlands", "Romania", "Singapore", "Spain", "Turkey", "United Kingdom", "United States"),
  total_population = c(25687041, 11589623, 213993437, 38005238, 19116201, 1411778724, 3882226, 5831404, 67391582, 83166711, 1393409038, 5022202, 9551000, 59066225, 125584839, 2794700, 36910560, 17475435, 19237691, 5703569, 47351567, 85319200, 67215293, 332915073)
)

scatter <- left_join(scatter, pop)

library(ggrepel)

color_mapping_region <- c(
  "Asia" = "#f08080",              # Light Coral
  "Asia-Pacific" = "#ffa500",       # Orange
  "Europe" = "#5f2680",             # Same as Protestant/All Others
  "Latin America" = "#800080",      # Sienna
  "Middle East" = "#4b0082",        # Indigo
  "Middle East & North Africa" = "#008000", # Green
  "North America" = "#add8e6"        # Light Blue
)




scatter %>% 
  ggplot(., aes(x = att, y = pray)) +
  geom_point(aes(color = region, size = total_population), stroke = 1, shape = 21) +
  geom_smooth(se = FALSE, method = lm, linetype = "twodash", color = "azure4") +
  theme_rb() + 
  scale_color_manual(values = color_mapping_region) +
  geom_text_repel(aes(x = att, y = pray, label = country), family = "font", size = 5) +
  labs(x = "Religious Attendance Index", y = "Prayer Frequency Index", title = "Relationship Between Religious Attendance and Prayer Frequency in 25 Countries", 
       caption = "@ryanburge + @religiondata\nData: Many-Analysts Religion Project, 2019")
save("many_att_pray_scatter.png")



country_colors <- c(
  "Australia" = "#ffa500",        # Asia-Pacific
  "Belgium" = "#5f2680",          # Europe
  "Brazil" = "#800080",           # Latin America
  "Canada" = "#add8e6",           # North America
  "Chile" = "#800080",            # Latin America
  "China" = "#f08080",            # Asia
  "Croatia" = "#5f2680",          # Europe
  "Denmark" = "#5f2680",          # Europe
  "France" = "#5f2680",           # Europe
  "Germany" = "#5f2680",          # Europe
  "India" = "#f08080",            # Asia
  "Ireland" = "#5f2680",          # Europe
  "Israel" = "#4b0082",           # Middle East
  "Italy" = "#5f2680",            # Europe
  "Japan" = "#ffa500",            # Asia-Pacific
  "Lithuania" = "#5f2680",        # Europe
  "Morocco" = "#008000",          # Middle East & North Africa
  "Netherlands" = "#5f2680",      # Europe
  "Romania" = "#5f2680",          # Europe
  "Singapore" = "#ffa500",        # Asia-Pacific
  "Spain" = "#5f2680",            # Europe
  "Turkey" = "#008000",           # Middle East & North Africa
  "United Kingdom" = "#5f2680",   # Europe
  "United States" = "#add8e6"     # North America
)


gg <- many %>%
  mutate(rel_7 = round(rel_5 * 6)) %>% 
  group_by(country) %>% 
  ct(rel_7)

mean_scores <- many %>%
  group_by(country) %>%
  mutate(rel_7 = round(rel_5 * 6)) %>% 
  summarise(mean_rel_7 = round(mean(rel_7, na.rm = TRUE), 1))


# Merge mean scores with gg to ensure alignment
gg <- gg %>%
  left_join(mean_scores, by = "country")

# Reorder the country factor in gg based on mean_rel_7
gg$country <- fct_reorder(gg$country, gg$mean_rel_7, .desc = TRUE)

# Plot with reordered facets and mean score annotations
gg %>% 
  ggplot(aes(x = rel_7, y = pct, fill = country)) +
  geom_col(color = "black") +
  facet_wrap(~ country) +  # Uses reordered country variable
  theme_rb() + 
  scale_fill_manual(values = country_colors) +
  scale_y_continuous(labels = scales::percent) +
  geom_text(aes(label = paste("Mean:", mean_rel_7), x = Inf, y = Inf), hjust = 2.25, vjust = 1.35, size = 4, family = "font") +
  labs(x = "", y = "", title = "To what extent do you believe in God? 0 = Not At All, 6 = Very Much", 
       caption = "@ryanburge + @religiondata\nData: Many-Analysts Religion Project, 2019")
save("many_believe_god.png", ht = 10)



gg <- many %>%
  mutate(rel_7 = round(rel_6 * 6)) %>% 
  mutate(rel = frcode(rel_7 == 0 ~ "Not\nAt All", 
                      rel_7 >= 1 & rel_7 <= 5 ~ "In\nBetween", 
                      rel_7 == 6 ~ "Very\nMuch")) %>% 
  group_by(country) %>% 
  ct(rel) 

gg %>% 
  mutate(lab = round(pct, 2)) %>% 
  ggplot(., aes(x = rel, y = pct, fill = country)) +
  geom_col(color = "black") +
  facet_wrap(~ country) +
  theme_rb() + 
  geom_text(aes(y = pct + .085, label = paste0(lab*100, '%')), position = position_dodge(width = .9), size = 6, family = "font") +
  scale_fill_manual(values = country_colors) +
  scale_y_continuous(labels = scales::percent) +
  labs(x = "", y = "", title = "To what extent do you believe in an afterlife?", 
       caption = "@ryanburge + @religiondata\nData: Many-Analysts Religion Project, 2019")
save("many_believe_afterlife.png", ht = 10)


gg <- many %>% 
  mutate(scale = (rel_1 + rel_2 + rel_5 + rel_6 + rel_7 + rel_8 + rel_9) / max(rel_1 + rel_2 + rel_5 + rel_6 + rel_7 + rel_8 + rel_9, na.rm = TRUE) * 100) %>% 
  select(country, scale) 


gg_meaned <- gg %>% 
  group_by(country) %>% 
  mutate(mean = mean(scale)) %>% 
  ungroup() 

library(ggridges)


gg_meaned %>% 
  ggplot() +
  geom_density_ridges(
    aes(x = scale, y = reorder(country, mean), fill = country), 
    quantile_lines = T, quantile_fun = mean, alpha = .5)  +   
  coord_cartesian(clip = "off") +
  theme_rb() +
  scale_fill_manual(values = country_colors) +
  theme(plot.title = element_text(size = 14)) +
  labs(x = "", y = "", title = "Distribution of Religiosity Index Scores", 
       caption = "@ryanburge + @religiondata\nData: Many-Analysts Religion Project, 2019")
save("many_rel_ridges.png", ht = 10)




gg <- many %>% 
  mutate(scale = (rel_1 + rel_2 + rel_5 + rel_6 + rel_7 + rel_8 + rel_9) / max(rel_1 + rel_2 + rel_5 + rel_6 + rel_7 + rel_8 + rel_9, na.rm = TRUE) * 100) %>% 
  mutate(psych = wb_psych_1 + vwb_psych_2 + wb_psych_3 + wb_psych_4 + wb_psych_5 + wb_psych_6/ max(wb_psych_1 + vwb_psych_2 + wb_psych_3 + wb_psych_4 + wb_psych_5 + wb_psych_6, na.rm = TRUE) * 100) %>% 
  select(country, scale, psych) %>% as_tibble()


gg %>% 
  ggplot(., aes(x = scale, y = psych, color = country)) +
  geom_point(stroke = .25, shape = 21, alpha = .5) + 
  geom_smooth(se = FALSE, method = lm, color = "black", linetype = "twodash") +
  scale_color_manual(values = country_colors) +
  facet_wrap(~ country) +
  theme_rb() + 
  labs(x = "Religiosity Scale", y = "Psychological Well-Being Scale", 
       title = "The Relationship Between Religiosity and Psychological Well-Being Across 25 Countries", 
       caption = "@ryanburge + @religiondata\nData: Many-Analysts Religion Project, 2019")
save("many_psych_rel_scales_scatter.png", ht = 10)


regression_results <- gg %>%
  group_by(country) %>%
  do(tidy(lm(psych ~ scale, data = .), conf.int = TRUE)) %>%
  filter(term == "scale") %>%
  select(country, estimate, conf.low, conf.high)

# Plot the coefficients with confidence intervals
regression_results %>%
  ggplot(aes(x = reorder(country, estimate), y = estimate, color = country)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) +
  geom_point(stroke = 1, shape = 21, fill = "white") +
  coord_flip() +
  theme_rb() +
  geom_hline(yintercept = 0) +
  scale_color_manual(values = country_colors) +
  labs(x = "Country", y = "Coefficient of Scale", title = "The Relationship Between the Religiosity Scale and the Psychological Well-Being Scale with 95% CI", 
       caption = "@ryanburge + @religiondata\nData: Many-Analysts Religion Project, 2019")
save("many_coefficients_two_scale.png")