Code and Artifacts for Mormon Typology

mormon_cluster.R

library(tidyverse)
library(haven)
library(car)
library(cluster)
library(pewmethods)
library(openxlsx)
library(fpc)

# Download 3 waves of the PRLS more info here: https://www.pewresearch.org/religious-landscape-study/
data07 <- haven::read_sav("./data/rls/rls2007.sav")
data14 <- haven::read_sav("./data/rls/rls2014.sav")
data24 <- haven::read_sav("./data/rls/rls2024.sav")

# I have gone through the codebooks of all 3 surveys and picked out 23 features
# to use for clustering analysis. The problem is the features are coded and
# named differently across the surveys. So, I will do the recoding and naming
# here.

# I created a google sheet where I map each question I use across the 3 surveys to a common set of 23 variables. 
# You can find it here: https://docs.google.com/spreadsheets/d/1E-kRTii0wMrwOlQlX9-if2I19k0O8l7Va4A35Lh8DAs/edit?usp=sharing

#######################
### Data Cleaning
#######################

# Step 1.1 Rename and recode features for consistency across waves

features <- c(
	'Happiness',
	'Health',
	'Family',
	'Homosexuality',
	'Right_and_Wrong',
	'Environment',
	'Bigger_Government',
    'Abortion',
    'Guidance',
    'Evolution',
    'Church_Attendance',
    'Religious_Importance',
    'Belief_in_God',
    'Certainty',
    'Heaven',
    'Hell',
    'Meditate',
    'Prayer_Frequency',
    'Church_Education',
    'Scriptures',
    'Proselyte',
    'Peace',
    'Wonder'
)

demos <- c(
    "age_rc",
    "educ_rc",
    "non_white",
    "female",
    "party_rc",
    "parent",
	"marital_rc"
)

other <- c(
	'weight',
	'respondent_id'
)

data14$weight <- data14$WEIGHT
data24$weight <- data24$WEIGHT

# Filtering to Mormons
data07 <- data07[data07$reltrad == 20000,]
data14 <- data14[data14$RELTRAD == 20000,]
data24 <- data24[data24$RELTRAD == 20000,]

# normalize weights
data07$weight <- data07$weight * (nrow(data07) / sum(data07$weight))
data14$weight <- data14$weight * (nrow(data14) / sum(data14$weight))
data24$weight <- data24$weight * (nrow(data24) / sum(data24$weight))

# respondent id
data07$respondent_id <- paste0("2007X", data07$psraid)
data14$respondent_id <- paste0("2014X", data14$resp)
data24$respondent_id <- paste0("2024X", data24$P_SUID)

###############################
# Happiness / Life Satisfaction
data07$Happiness <- dplyr::case_when(
	data07$q2 == 1 & data07$q2a == 1 ~ 1,
	data07$q2 == 1 & data07$q2a == 2 ~ 2,
	data07$q2 == 2 ~ 3,
	TRUE ~ NA_real_
)

data14$Happiness <- ifelse(data14$qa1 == 9, NA, data14$qa1)

data24$Happiness <- ifelse(data24$HAPPY == 99, NA, data24$HAPPY)

###############################
# Health Satisfaction
data07$Health <- ifelse(data07$q3a == 9, NA, data07$q3a)

data14$Health <- ifelse(data14$qa2b == 9, NA, data14$qa2b)

data24$Health <- dplyr::case_when(
	data24$SATIS_A == 1 ~ 1,
	data24$SATIS_A == 2 ~ 2,
	data24$SATIS_A == 3 ~ 2,
	data24$SATIS_A == 4 ~ 3,
	data24$SATIS_A == 5 ~ 4,
	data24$SATIS_A == 99 ~ NA_real_
)

###############################
# Family Satisfaction
data07$Family <- ifelse(data07$q3b == 9, NA, data07$q3b)

data14$Family <- ifelse(data14$qa2a == 9, NA, data14$qa2a)

data24$Family <- dplyr::case_when(
	data24$SATIS_B == 1 ~ 1,
	data24$SATIS_B == 2 ~ 2,
	data24$SATIS_B == 3 ~ 2,
	data24$SATIS_B == 4 ~ 3,
	data24$SATIS_B == 5 ~ 4,
	data24$SATIS_B == 99 ~ NA_real_
)

###############################
# Homosexuality
data07$Homosexuality <- dplyr::case_when(
	data07$q5b == 1 ~ 1,
	data07$q5b == 2 ~ 2,
	data07$q5b == 3 ~ NA_real_,
	TRUE ~ NA_real_
)

data14$Homosexuality <- dplyr::case_when(
	data14$qb2a == 1 ~ 1,
	data14$qb2a == 2 ~ 2,
	data14$qb2a == 3 ~ NA_real_,
	TRUE ~ NA_real_
)

data24$Homosexuality <- dplyr::case_when(
	data24$QB2A == 1 ~ 1,
	data24$QB2A == 2 ~ 2,
	TRUE ~ NA_real_
)

###############################
# Right And Wrong

data07$Right_and_Wrong <- dplyr::case_when(
	data07$q10b == 1 ~ 1,
	data07$q10b == 2 ~ 1,
	data07$q10b == 3 ~ 2,
	data07$q10b == 4 ~ 2,
	TRUE ~ NA_real_
)

data14$Right_and_Wrong <- dplyr::case_when(
	data14$qb2d == 1 ~ 1,
	data14$qb2d == 2 ~ 2,
	TRUE ~ NA_real_
)

data24$Right_and_Wrong <- dplyr::case_when(
	data24$QB2D == 1 ~ 1,
	data24$QB2D == 2 ~ 2,
	TRUE ~ NA_real_
)

###############################
# Environment
data07$Environment <- dplyr::case_when(
	data07$q5e == 1 ~ 1,
	data07$q5e == 2 ~ 2,
	TRUE ~ NA_real_
)

data14$Environment <- dplyr::case_when(
	data14$qb2c == 1 ~ 1,
	data14$qb2c == 2 ~ 2,
	TRUE ~ NA_real_
)

data24$Environment <- dplyr::case_when(
	data24$QB2C == 1 ~ 1,
	data24$QB2C == 2 ~ 2,
	TRUE ~ NA_real_
)

###############################
# Bigger Government
data07$Bigger_Government <- dplyr::case_when(
	data07$q6 == 1 ~ 1,
	data07$q6 == 2 ~ 2,
	TRUE ~ NA_real_
)

data14$Bigger_Government <- dplyr::case_when(
	data14$qb20 == 1 ~ 1,
	data14$qb20 == 2 ~ 2,
	TRUE ~ NA_real_
)

data24$Bigger_Government <- dplyr::case_when(
	data24$GOVSIZE1 == 1 ~ 1,
	data24$GOVSIZE1 == 2 ~ 2,
	TRUE ~ NA_real_
)

###############################
# Abortion

data07$Abortion <- ifelse(data07$q7 == 9, NA, data07$q7)

data14$Abortion <- ifelse(data14$qb21 == 9, NA, data14$qb21)

data24$Abortion <- ifelse(data24$ABRTLGL == 99, NA, data24$ABRTLGL)

###############################
# Guidance

data07$Guidance <- ifelse(data07$q10d == 1, 1, ifelse(
	data07$q10d == 9, NA,
	0
))

data14$Guidance <- ifelse(data14$qb31 == 1, 1, ifelse(
	data14$qb31 == 9, NA,
	0
))

data24$Guidance <- ifelse(data24$GUIDE_A == 1, 1, ifelse(
	data24$GUIDE_A == 99, NA,
	0
))

###############################
# Evolution

data07$Evolution <- dplyr::case_when(
	data07$q10c == 1 ~ 1,
	data07$q10c == 2 ~ 1,
	data07$q10c == 3 ~ 2,
	data07$q10c == 4 ~ 2,
	data07$q10c == 9 ~ NA_real_
)

data14$Evolution <- ifelse(data14$qb30 == 9, NA, data14$qb30)

data24$Evolution <- dplyr::case_when(
	data24$EVOL == 1 ~ 1,
	data24$EVOL == 2 ~ 1,
	data24$EVOL == 3 ~ 2,
	data24$EVOL == 99 ~ NA_real_
)

###############################
# Church Attendance

data07$Church_Attendance <- ifelse(data07$q20 == 9, NA, data07$q20)

data14$Church_Attendance <- ifelse(data14$attend == 9, NA, data14$attend)

data24$Church_Attendance <- ifelse(data24$ATTNDPERRLS == 99, NA, data24$ATTNDPERRLS)

###############################
# Religious Importance

data07$Religious_Importance <- ifelse(data07$q21 == 9, NA, data07$q21)

data14$Religious_Importance <- ifelse(data14$qf2 == 9, NA, data14$qf2)

data24$Religious_Importance <- ifelse(data24$RELIMP == 99, NA, data24$RELIMP)

###############################
# Belief in God

data07$Belief_in_God <- ifelse(data07$q30 == 9, NA, ifelse(data07$q30 == 3, NA, data07$q30))

data14$Belief_in_God <- ifelse(data14$qg1 == 9, NA, ifelse(data14$qg1 == 3, NA, data14$qg1))

data24$Belief_in_God <- ifelse(data24$GOD == 99, NA, data24$GOD)

###############################
# Certainty

data07$Certainty <- ifelse(data07$q31 == 9, NA, data07$q31)

data14$Certainty <- ifelse(data14$qg1b == 9, NA, data14$qg1b)

data24$Certainty <- ifelse(data24$GOD2 == 99, NA, data24$GOD2)

###############################
# Heaven
data07$Heaven <- ifelse(data07$q35 == 9, NA, ifelse(data07$q35 == 3, NA, data07$q35))

data14$Heaven <- ifelse(data14$qg5 == 9, NA, ifelse(data14$qg5 == 3, NA, data14$qg5))

data24$Heaven <- ifelse(data24$HVN == 99, NA, data24$HVN)

###############################
# Hell

data07$Hell <- ifelse(data07$q36 == 9, NA, ifelse(data07$q36 == 3, NA, data07$q36))

data14$Hell <- ifelse(data14$qg6 == 9, NA, ifelse(data14$qg6 == 3, NA, data14$qg6))

data24$Hell <- ifelse(data24$HLL == 99, NA, data24$HLL)

###############################
# Meditate

data07$Meditate <- dplyr::case_when(
	data07$q42c == 1 ~ 1,
	data07$q42c == 2 ~ 2,
	data07$q42c == 3 ~ 3,
	data07$q42c == 4 ~ 4,
	data07$q42c == 5 ~ 4,
	data07$q42c == 9 ~ NA_real_
)

data14$Meditate <- dplyr::case_when(
	data14$qi2c == 1 ~ 1,
	data14$qi2c == 2 ~ 2,
	data14$qi2c == 3 ~ 3,
	data14$qi2c == 4 ~ 4,
	data14$qi2c == 5 ~ 4,
	data14$qi2c == 9 ~ NA_real_
)

data24$Meditate <- dplyr::case_when(
	data24$SPIRACT_F == 1 ~ 1,
	data24$SPIRACT_F == 2 ~ 1,
	data24$SPIRACT_F == 3 ~ 2,
	data24$SPIRACT_F == 4 ~ 3,
	data24$SPIRACT_F == 5 ~ 4,
	data24$SPIRACT_F == 99 ~ NA_real_
)

# Adding quick function for NA cleaning of variables
list_to_it <- list(
	"Prayer_Frequency" = list('2007' = "q41", '2014' = "qi1", '2024' = "PRAY"),
	"Church_Education" = list('2007' = "q42a", '2014' = "qi2a", '2024' = "PRAC_A"),
	"Scriptures" = list('2007' = "q42b", '2014' = "qi2b", '2024' = "PRAC_B"),
	"Proselyte" = list('2007' = "q42d", '2014' = "qi2g", '2024' = "PRAC_D"),
	"Peace" = list('2007' = "q43a", '2014' = "qi4a", '2024' = "EXP_A"),
	"Wonder" = list('2007' = "q43b", '2014' = "qi4b", '2024' = "EXP_B")
)

data07$marital_rc <- as.numeric(car::recode(data07$marital, "4:5=2;6=4;9=NA"))
data07$age_rc <- as.numeric(car::recode(
    as.numeric(data07$age), 
    "1:35=1;36:45=2;46:64=3;65:97=4;99=NA", as.numeric = FALSE))
data07$educ_rc <- as.numeric(car::recode(data07$educ, "1:3=1;4:5=2;6=3;7=4;9=NA"))
data07$non_white <- as.numeric(car::recode(data07$race, "1=0;2:5=1;9=NA"))
data07$female <- as.numeric(car::recode(data07$sex, "1=0;2=1;9=NA"))
data07$party_rc <- dplyr::case_when(
    data07$party %in% c(1) | (data07$partyln  %in% c(1)) ~ 1,
    data07$party %in% c(2) | (data07$partyln  %in% c(2)) ~ 3,
    data07$party %in% c(4,5) ~ 2,
    TRUE ~ NA
)
data07$parent <- as.numeric(car::recode(data07$children, "0=0;1:50=1;99=NA"))


data14$marital_rc <- as.numeric(car::recode(data14$marital, "4:5=2;6=4;9=NA"))
data14$age_rc <- as.numeric(car::recode(data14$agerec, "1:3=1;4:5=2;6:9=3;10:15=4;99=NA"))
data14$educ_rc <- as.numeric(car::recode(data14$educ, "1:3=1;4:5=2;6=3;7:8=4;9=NA"))
data14$non_white <- as.numeric(car::recode(data14$racethn, "1=0;2:4=1;9=NA"))
data14$female <- as.numeric(car::recode(data14$SEX, "1=0;2=1"))
data14$fertrec_rc <- as.numeric(car::recode(data14$fertrec, "1:6=1;99=NA"))
data14$party_rc <- as.numeric(dplyr::case_when(
    data14$party %in% c(1) | (data14$party  %in% c(3) & data14$partyln  %in% c(1)) ~ 1,
    data14$party %in% c(2) | (data14$party  %in% c(3) & data14$partyln  %in% c(2)) ~ 3,
    data14$party %in% c(4,5) ~ 2,
    TRUE ~ NA
))
data14$parent <- as.numeric(car::recode(data14$childrenrec, "0=0;1:50=1;99=NA"))


data24$marital_rc <- as.numeric(car::recode(data24$MARITAL, "4:5=2;6=4;99=NA"))
data24$age_rc <- as.numeric(car::recode(
    as.numeric(data24$BIRTHDECADE), 
    "6:7=1;5=2;3:4=3;1:2=4;99=NA"
))
data24$educ_rc <- as.numeric(car::recode(data24$EDUCREC, "99=NA"))
data24$non_white <- as.numeric(car::recode(data24$RACECMB, "1=0;2:5=1;9=NA"))
data24$GENDER <- as.numeric(car::recode(data24$GENDER, "99=NA"))
data24$female <- as.numeric(car::recode(data24$GENDER, "1=0;2=1;3=NA;99=NA"))
data24$inc_rc <- as.numeric(car::recode(data24$INC_SDT1, "1:3=1;4:5=2;6:8=3;99=NA"))
data24$parent <- as.numeric(car::recode(data24$CHILDREN, "1=1;2=0;99=NA"))

data24$party_rc <- as.numeric(dplyr::case_when(
    data24$PARTY %in% c(1) | (data24$PARTY  %in% c(3) & data24$PARTYLN  %in% c(1)) ~ 1,
    data24$PARTY %in% c(2) | (data24$PARTY  %in% c(3) & data24$PARTYLN  %in% c(2)) ~ 3,
    data24$PARTY %in% c(4) ~ 2,
    TRUE ~ NA
))

replace_99_9_w_na <- function(data07, data14, data24, instr){
	for (var in names(instr)){
        print(var)
        data07[[var]] <- ifelse(data07[[instr[[var]][['2007']]]] == 9, NA, data07[[instr[[var]][['2007']]]])
        data14[[var]] <- ifelse(data14[[instr[[var]][['2014']]]] == 9, NA, data14[[instr[[var]][['2014']]]])
        data24[[var]] <- ifelse(data24[[instr[[var]][['2024']]]] == 99, NA, data24[[instr[[var]][['2024']]]])
	}
	return(list(data07 = data07, data14 = data14, data24 = data24))
}

res <- replace_99_9_w_na(data07, data14, data24, list_to_it)

# Step 1.2 Combine Data

data07 <- res$data07[c(features,demos,other)] |> mutate(year = 2007) |> mutate(across(-respondent_id, as.numeric))
data14 <- res$data14[c(features,demos,other)] |> mutate(year = 2014) |> mutate(across(-respondent_id, as.numeric))
data24 <- res$data24[c(features,demos,other)] |> mutate(year = 2024) |> mutate(across(-respondent_id, as.numeric))

data <- bind_rows(data07, data14, data24)

# normalize weights
data$weight <- data$weight * (nrow(data) / sum(data$weight))

#######################
### Clustering Set Up & Implementation
#######################

# Step 2.1 Imputing and Scaling

X <- data |> 
	select(features) |>
	#impute at the mean
	mutate(across(everything(), function(x) ifelse(is.na(x), mean(x, na.rm = TRUE), x))) |> 
	#scale the features
	scale(center = TRUE, scale = TRUE)

# Step 2.2 Dimensionality Reduction

res_pca <- prcomp(X, center = TRUE, scale. = TRUE)
summary(res_pca)
# plot(res_pca, type = "l")

pca_data <- res_pca$x[, 1:10] # first 10 explain 70% of variance; 1st explains 25%

# Step 2.3 Clustering

set.seed(123)  # for reproducibility
wss <- sapply(1:8, function(k){
	kmeans(pca_data, centers = k, nstart = 25)$tot.withinss
})
# did not converge at 1 obviously

# plot(
# 	1:8, 
# 	wss, 
# 	type = "b", 
# 	xlab = "Number of clusters", 
# 	ylab = "Within groups sum of squares"
# )
# There is a clear elbow at k=2, but 3 and 4 still explain some variance before mostly flattening after that.


# Step 2.4 Second Diagnostic Check
sil_width <- c()
for(k in 2:8){
	km <- kmeans(pca_data, centers = k, nstart = 25)
	sil <- cluster::silhouette(km$cluster, dist(pca_data))
	sil_width[k] <- mean(sil[, 3])  # average silhouette width
}

# plot(
# 	2:8, 
# 	sil_width[2:8], 
# 	type = "b",
# 	xlab = "Number of clusters (k)",
# 	ylab = "Average silhouette width"
# )

# Again highest silhouette width at k=2 followed by 4 then 3.
# I'll run tabs at 2,3,4 and see what breaks out.
# After running tabs and analysing adding an extra step....

kmeans_res <- kmeans(pca_data, centers = 2, nstart = 25)
data$cluster_2 <- factor(kmeans_res$cluster)

kmeans_res <- kmeans(pca_data, centers = 3, nstart = 25)
data$cluster_3 <- factor(kmeans_res$cluster)

kmeans_res <- kmeans(pca_data, centers = 4, nstart = 25)
data$cluster_4 <- factor(kmeans_res$cluster)

#######################
### Packaging Results
#######################

# # Step 3.1 Run Crosstabulations

vars_to_it <- c(demos, features)
cluster_vars <- c("cluster_2", "cluster_3", "cluster_4")
xtabs <- purrr::map(cluster_vars, function(cluster_var){
	purrr::map(vars_to_it, ~get_totals(.x, data, wt = "weight",
    by = cluster_var, by_total = TRUE, digits = 1))
})

vars_to_it <- c(demos, features)
cluster_vars <- c("cluster_2", "cluster_3", "cluster_4")
xtabs_2024 <- purrr::map(cluster_vars, function(cluster_var){
	purrr::map(vars_to_it, ~get_totals(.x, data |> filter(year == 2024), wt = "weight",
    by = cluster_var, by_total = TRUE, digits = 1))
})

# # Step 3.2 Save to Excel

pewmethods::df_list_to_xlsx(
	df_list = xtabs,
	sheet_name = cluster_vars,
	title = cluster_vars,
	outfile = "./tabs/15_lds_cluster_tabs.xlsx"
)

pewmethods::df_list_to_xlsx(
	df_list = xtabs_2024,
	sheet_name = cluster_vars,
	title = cluster_vars,
	outfile = "./tabs/15_lds_cluster_tabs_2024.xlsx"
)

# # Step 3.3 Export Clusters

data <- data |> 
	mutate(mormon_typology = case_when(
		cluster_4 == 1 ~ "Cultural Mormons",
		cluster_4 == 2 ~ "Devout Traditionalists",
		cluster_4 == 3 ~ "In-Betweeners",
		cluster_4 == 4 ~ "Adaptive Believers"
	))

mormon_typology_data <- data |> 
	select(respondent_id, mormon_typology)

write.csv(mormon_typology_data, "./data/rls/mormon_typology.csv", row.names = FALSE)

#######################
### Additional Robustness Checks
#######################

# Step 4.1 I have a group that is small between 2-6% testing robustness...

# for (i in 1:40) {
# 	idx <- sample(1:1810, 1000)
# 	data_copy <- data[idx,]
# 	boot_res <- kmeans(pca_data[idx,], centers = 4, nstart = 25)
# 	data_copy$cluster_4 <- factor(boot_res$cluster)
# 	print(get_totals('cluster_4', data_copy, by = 'year', wt = "weight"))
# 	print(get_totals('Belief_in_God', data_copy, by = 'cluster_4', wt = "weight"))
# 	print(get_totals('non_white', data_copy, by = 'cluster_4', wt = "weight"))
# }

# # This small group consistently captures disbelief in god and often more nonwhite

# bootstab <- fpc::clusterboot(
# 	pca_data, 
# 	B=500, 
# 	clustermethod=kmeansCBI,
#     krange=4, 
# 	seed=123, 
# 	count=FALSE, 
# 	nstart=50
# )
# bootstab$bootmean  # stability measures per cluster

# orig_clusters <- factor(bootstab$result$partition, levels = 1:4)
# levels(orig_clusters) <- c(
# 	"Devout Traditionalists",
# 	"Cultural but Distant",
# 	"Thoughtful Believers",
# 	"Emerging & Questioning"
# )
# prop.table(table(orig_clusters)) * 100

# # all clusters are stable and strong according to bootstrapped Jaccard 0.78–0.91
# # including the small one in question

# set.seed(123)
# sizes <- numeric(500)
# for (b in 1:500) {
# 	idx <- sample(seq_len(nrow(pca_data)), replace = TRUE)
# 	km  <- kmeans(pca_data[idx,], centers = 4, nstart = 50, iter.max = 100)
# 	# size of the smallest cluster in this resample:
# 	sizes[b] <- min(table(km$cluster))/length(idx)
# }
# summary(sizes)
# mean(sizes <= 0.10)

# # 95% of the time, the smallest cluster is smaller than 10% of the sample size.

mormon_typology.csv