#----- Religious Nones - GSS data 1972 - 2018
# https://twitter.com/pjmccann3/status/1109127955907837952
# https://t.co/3tBqjPBXeb  dta download url
# D:\Downloads\gss_full.dta
#rm(list=ls())

library(tidyverse)

d = read_dta("D:/Downloads/gss_full.dta")

saveRDS(d, "C:/Data/R/religous_nones_gss_data.rds")  # CYA save
# FYI, the *.dta file was 445 MB while the *.rds file was a mere 33 MB.  R compression FTW!

#----- Get raw data of desired columns -----
datraw <- select(d, year, evangelical, mainline, blackprot, catholic, jewish, otherfaith, nofaith)

#----- counts by year -----
datcount <- d %>%
  group_by(year) %>%
  summarize(
    total = sum(evangelical + mainline + blackprot + catholic + jewish + otherfaith + nofaith),
    evangelical = sum(evangelical),
    mainline = sum(mainline) ,
    blackprot = sum(blackprot),
    catholic = sum(catholic) ,
    jewish = sum(jewish) ,
    otherfaith = sum(otherfaith),
    nofaith = sum(nofaith)
  )

#----- percents by year -----
datpercent <- d %>%
  group_by(year) %>%
  summarize(
    total = sum(evangelical + mainline + blackprot + catholic + jewish + otherfaith + nofaith),        
    evangelical = sum(evangelical) / total,
            mainline = sum(mainline) / total,
            blackprot = sum(blackprot) / total,
            catholic = sum(catholic) / total,
            jewish = sum(jewish) / total,
            otherfaith = sum(otherfaith) / total,
            nofaith = sum(nofaith) / total
         )

#----- save off data as TSV files -----
write.table(datraw, file = "C:/Data/R/religous_nones_gss_data_raw.tsv", sep = "\t", row.names = FALSE, quote = FALSE)
write.table(datcount, file = "C:/Data/R/religous_nones_gss_data_count.tsv", sep = "\t", row.names = FALSE, quote = FALSE)
write.table(datpercent, file = "C:/Data/R/religous_nones_gss_data_percent.tsv", sep = "\t", row.names = FALSE, quote = FALSE)

# tidy up the dataset for ggplot
tidydat <- gather(dpercent, key="denomination", value="count",  -year, -total)

ggplot(tidydat, aes(x = year, y = count)) +
  geom_line(aes(color=denomination))