# # Description of the development of religious denomination, religious attendance, prayer, and religiosity

library(magrittr)
library(memisc)



load("ESScumul.RData")


ESScumul %<>% within({
    reldenom[reldenom.any == 2] <- 9
    labels(reldenom) <- labels(reldenom) + c('No religion'=9)
    denom <- recode(reldenom,
                    Catholic              = 1 <- 1,
                    Protestant            = 2 <- 2,
                    'Other Christian'     = 3 <- 3:4,
                    Islamic               = 4 <- 6,
                    'Other Non-Christian' = 5 <- c(5,7:8),
                    'No religion'         = 6 <- 9,
                    copy=TRUE
                    )
    western.europe <- country %in% c(
                                   "Austria",
                                   "Belgium",
                                   "Denmark",
                                   "France",
                                   "Germany (West)",
                                   "Ireland",
                                   "Italy",
                                   "Netherlands",
                                   "Norway",
                                   "Portugal",
                                   "Spain",
                                   "Sweden",
                                   "Switzerland",
                                   "United Kingdom"
                                 )
    # There were obviously problems in measuring religious denomination
    # in these countries and years
    denom[country=="France" & year <= 2005] <- NA
    denom[country=="Denmark" & year >= 2019] <- NA
    denom[country=="Sweden" & year == 2015] <- NA
    denom[country=="United Kingdom" & year >= 2004 & year <= 2007 ] <- NA
    relatt <- recode(relig.attend,
        Weekly = 5 <- 1:3,
        Monthly = 4 <- 4,
        Holidays = 3 <- 5,
        Rarely = 2 <- 6,
        Never = 1 <- 7,
        copy = TRUE
    )
    # There were obviously problems in measuring religious attendance
    # in these countries and years
    relatt[(country == "Sweden" & year == 2015) |
           (country == "Norway" & year == 2017) |
           (country == "Ireland" & year == 2021) | 
            year == 2020] <- NA
    pray[(country == "Sweden" & year == 2015) |
         (country == "Norway" & year == 2017)] <- NA
    religious[(country == "Norway" & year == 2017)] <- NA
})


ESSwork <- subset(ESScumul,
                  select=c(
                      vote,
                      denom,
                      religious,
                      relig.attend,
                      pray,
                      gender,
                      age,
                      year,
                      ybirth,
                      w.design,
                      w.poststrat,
                      id,
                      CHES_id,
                      oesch9,
                      oesch6,
                      country,
                      western.europe
                  ))


revlevels <- function(x){
    factor(x,levels=rev(levels(x)))
}


perc.denom.country.year <- percentages(denom~country+year,
                                          data=ESSwork,
                                          #weights=w.design,
                                          subset=western.europe==1)


perc.denom.country.year.table <- perc.denom.country.year
perc.denom.country.year %<>% as.data.frame()
perc.denom.country.year %<>% subset(is.finite(Percentage))
perc.denom.country.year %<>% within({
    year <- as.numeric(as.character(year))
})


library(ggplot2)


str(perc.denom.country.year)


library(RColorBrewer)


#pal <- c(brewer.pal(8,"Accent") [c(1:5)],few_pal("Dark")(2)[2:1])
pal <- brewer.pal(8,"Accent") [c(1:2,7,3:5)]
#pal <- few_pal("Light")(7)


options(jupyter.plot.width=9,
        jupyter.plot.height=6)
ggplot(perc.denom.country.year,aes(x=year,y=Percentage)) + 
  geom_bar(aes(fill=revlevels(denom)),stat="identity",width=1.1) +
  scale_fill_manual(values=pal) +
  #scale_fill_grey(start=.9,end=.4) +
  #scale_fill_brewer(palette = "Accent") +
  #scale_fill_viridis_d(begin=1,end=0) +
  facet_wrap(~country,nrow=3) + theme_bw() + 
  labs(x="Year of fieldwork",
       fill="Regious denomination")


ggsave("figures/percentages-denomination-time-countries-bars.pdf")
ggsave("figures/percentages-denomination-time-countries-bars.png")


options(jupyter.plot.width=9,
        jupyter.plot.height=6)
ggplot(perc.denom.country.year,aes(x=year,y=Percentage)) + 
  geom_area(aes(fill=revlevels(denom))) +
  scale_fill_manual(values=pal) +
  #scale_fill_grey(start=.9,end=.4) +
  #scale_fill_brewer(palette = "Accent") +
  #scale_fill_viridis_d(begin=1,end=0) +
  facet_wrap(~country,nrow=3) + theme_bw() + 
  labs(x="Year of fieldwork",
       fill="Regious denomination")


ggsave("figures/percentages-denomination-time-countries-areas.pdf")
ggsave("figures/percentages-denomination-time-countries-areas.png")


options(jupyter.plot.width=9,
        jupyter.plot.height=6)
ggplot(perc.denom.country.year,aes(x=year,y=Percentage)) + 
  geom_area(aes(fill=revlevels(denom)),stat = "smooth",span=2) +
  scale_fill_manual(values=pal) +
  #scale_fill_grey(start=.9,end=.4) +
  #scale_fill_brewer(palette = "Accent") +
  # scale_fill_brewer(palette = "Dark2") +
  #scale_fill_viridis_d(option = "inferno") +
  #scale_fill_few() +
  #scale_fill_viridis_d(begin=1,end=0) +
  facet_wrap(~country,nrow=3) + theme_bw() + 
  labs(x="Year of fieldwork",
       fill="Regious denomination")


ggsave("figures/percentages-denomination-time-countries-areas-smooth.pdf")
ggsave("figures/percentages-denomination-time-countries-areas-smooth.png")


perc.denom.country <- percentages(denom~country,
                                          data=ESScumul,
                                          weights=w.design,
                                          subset=western.europe==1)


perc.denom.country.table <- perc.denom.country
perc.denom.country %<>% as.data.frame()
perc.denom.country %<>% subset(is.finite(Percentage))


perc.denom.country.wide <- to.data.frame(perc.denom.country.table)


perc.denom.country.wide %<>% subset(is.finite(Protestant))


perc.denom.country.wide %$$% {
    Cath_Prot <- 100*(Catholic/(Catholic+Protestant))
}
perc.denom.country.wide %<>% rename('No religion'='No_religion')


perc.denom.country.wide %$% Cor(cbind(Catholic,No_religion,Cath_Prot))


perc.denom.country.wide


save(perc.denom.country.wide, file="perc-denom-country-wide.RData")


perc.relatt.country.year <- percentages(relatt~country+year,
                                          data=ESScumul,
                                          weights=w.design,
                                          subset=western.europe==1 & country != "Luxembourg" & 
                                          !(country == "Sweden" & year == 2015) &
                                          !(country == "Norway" & year == 2017) &
                                          !(country == "Ireland" & year == 2021) &
                                          as.numeric(denom) %in% 1:3
                                       )
perc.relatt.country.year %<>% as.data.frame()
perc.relatt.country.year %<>% subset(is.finite(Percentage))
perc.relatt.country.year %<>% within({
    year <- as.numeric(as.character(year))
})


options(jupyter.plot.width=9,
        jupyter.plot.height=6)
ggplot(perc.relatt.country.year,aes(x=year,y=Percentage)) + 
  geom_area(aes(fill=relatt)) +
  scale_fill_brewer() +
  facet_wrap(~country,nrow=3) + theme_bw() + 
  labs(x="Year of fieldwork",
       fill="Attends religious\nservices")


ggsave("figures/percentages-religatt-time-countries-areas.pdf")
ggsave("figures/percentages-religatt-time-countries-areas.png")


options(jupyter.plot.width=9,
        jupyter.plot.height=6)
ggplot(perc.relatt.country.year,aes(x=year,y=Percentage)) + 
  geom_area(aes(fill=relatt),stat = "smooth",span=2) +
  scale_fill_brewer() +
  facet_wrap(~country,nrow=3) + theme_bw() + 
  labs(x="Year of fieldwork",
       fill="Attends religious\nservices")


ggsave("figures/percentages-religatt-time-countries-areas-smooth.pdf")
ggsave("figures/percentages-religatt-time-countries-areas-smooth.png")


perc.pray.country.year <- percentages(pray~country+year,
                                          data=ESScumul,
                                          weights=w.design,
                                          subset=western.europe==1 & country != "Luxembourg"  & 
                                          !(country == "Sweden" & year == 2015) &
                                          !(country == "Norway" & year == 2017) &
                                          as.numeric(denom) %in% 1:3)
perc.pray.country.year %<>% as.data.frame()
perc.pray.country.year %<>% subset(is.finite(Percentage))
perc.pray.country.year %<>% within({
    year <- as.numeric(as.character(year))
})


options(jupyter.plot.width=9,
        jupyter.plot.height=6)
ggplot(perc.pray.country.year,aes(x=year,y=Percentage)) + 
  geom_area(aes(fill=pray)) +
  scale_fill_brewer() +
  facet_wrap(~country,nrow=3) + theme_bw() + 
  labs(x="Year of fieldwork",
       fill="Prays how often")


ggsave("figures/percentages-prayer-time-countries-areas.pdf")
ggsave("figures/percentages-prayer-time-countries-areas.png")


options(jupyter.plot.width=9,
        jupyter.plot.height=6)
ggplot(perc.pray.country.year,aes(x=year,y=Percentage)) + 
  geom_area(aes(fill=pray),stat = "smooth",span=2) +
  scale_fill_brewer() +
  facet_wrap(~country,nrow=3) + theme_bw() + 
  labs(x="Year of fieldwork",
       fill="Prays how often")


ggsave("figures/percentages-prayer-time-countries-areas-smooth.pdf")
ggsave("figures/percentages-prayer-time-countries-areas-smooth.png")


perc.religious.country.year <- percentages(religious~country+year,
                                          data=ESScumul,
                                          weights=w.design,
                                          subset=western.europe==1 
                                           & country != "Luxembourg" 
                                           & !(country == "Norway" & year == 2017)
                                           & as.numeric(denom) %in% 1:3)
perc.religious.country.year %<>% as.data.frame()
perc.religious.country.year %<>% subset(is.finite(Percentage))
perc.religious.country.year %<>% within({
    year <- as.numeric(as.character(year))
})


startcol <- RColorBrewer::brewer.pal(9,"Blues")[1]
endcol <- RColorBrewer::brewer.pal(9,"Blues")[9]


Blues11 <- colorRampPalette(RColorBrewer::brewer.pal(9,"Blues")[c(1,9)])(11)


options(jupyter.plot.width=12,
        jupyter.plot.height=10)
ggplot(perc.religious.country.year,aes(x=year,y=Percentage)) + 
  geom_area(aes(fill=religious)) +
  scale_fill_manual(values=Blues11) +
  facet_wrap(~country) + theme_bw() + 
  labs(x="Year of fieldwork",
       fill="How religious")


ggsave("figures/percentages-religious-time-countries-areas.pdf")
ggsave("figures/percentages-religious-time-countries-areas.png")


perc.relatt.denom.year <- percentages(relatt~denom+year,
                                          data=ESScumul,
                                          weights=w.design,
                                          subset=western.europe==1 
                                          & country != "Luxembourg"
                                          & year != 2020)
perc.relatt.denom.year %<>% as.data.frame()
perc.relatt.denom.year %<>% subset(is.finite(Percentage))
perc.relatt.denom.year %<>% within({
    year <- as.numeric(as.character(year))
    #denom5 <- factor(denom5,levels=levels(denom5)[c(1,2,5,3,4)])
})


options(jupyter.plot.width=9,
        jupyter.plot.height=6)
ggplot(perc.relatt.denom.year,aes(x=year,y=Percentage)) + 
  geom_area(aes(fill=relatt),stat = "identity") +
  scale_fill_brewer() +
  facet_wrap(~denom,nrow=2) + theme_bw() + 
  labs(x="Year of fieldwork",
       fill="Attends religious\nservices")


ggsave("figures/percentages-religatt-denomination-time-areas.pdf")
ggsave("figures/percentages-religatt-denomination-time-areas.png")


options(jupyter.plot.width=9,
        jupyter.plot.height=6)
ggplot(subset(perc.relatt.denom.year,!(year == 2020)),aes(x=year,y=Percentage)) + 
  geom_area(aes(fill=relatt),stat = "smooth",span=2) +
  scale_fill_brewer() +
  facet_wrap(~denom,nrow=2) + theme_bw() + 
  labs(x="Year of fieldwork",
       fill="Attends religious\nservices")


ggsave("figures/percentages-religatt-denomination-time-areas-smooth.pdf")
ggsave("figures/percentages-religatt-denomination-time-areas-smooth.png")


perc.pray.denom.year <- percentages(pray~denom+year,
                                          data=ESScumul,
                                          weights=w.design,
                                          subset=western.europe==1 
                                          & country != "Luxembourg"
                                          & year != 2020)
perc.pray.denom.year %<>% as.data.frame()
perc.pray.denom.year %<>% subset(is.finite(Percentage))
perc.pray.denom.year %<>% within({
    year <- as.numeric(as.character(year))
    #denom5 <- factor(denom5,levels=levels(denom5)[c(1,2,5,3,4)])
})


options(jupyter.plot.width=9,
        jupyter.plot.height=6)
ggplot(perc.pray.denom.year,aes(x=year,y=Percentage)) + 
  geom_area(aes(fill=pray),stat = "identity") +
  scale_fill_brewer() +
  facet_wrap(~denom,nrow=2) + theme_bw() + 
  labs(x="Year of fieldwork",
       fill="Prays how often")


ggsave("figures/percentages-prayer-denomination-time-areas.pdf")
ggsave("figures/percentages-prayer-denomination-time-areas.png")


options(jupyter.plot.width=9,
        jupyter.plot.height=6)
ggplot(subset(perc.pray.denom.year,!(year == 2020)),
       aes(x=year,y=Percentage)) + 
  geom_area(aes(fill=pray),stat = "smooth",span=2) +
  scale_fill_brewer() +
  facet_wrap(~denom,nrow=2) + theme_bw() + 
  labs(x="Year of fieldwork",
       fill="Prays how often")


ggsave("figures/percentages-prayer-denomination-time-areas-smooth.pdf")
ggsave("figures/percentages-prayer-denomination-time-areas-smooth.png")


perc.religious.denom.year <- percentages(religious~denom+year,
                                          data=ESScumul,
                                          weights=w.design,
                                          subset=western.europe==1 
                                          & country != "Luxembourg"
                                          & year != 2020)
perc.religious.denom.year %<>% as.data.frame()
perc.religious.denom.year %<>% subset(is.finite(Percentage))
perc.religious.denom.year %<>% within({
    year <- as.numeric(as.character(year))
    #denom <- factor(denom,levels=levels(denom)[c(1,2,5,3,4)])
})


options(jupyter.plot.width=9,
        jupyter.plot.height=6)
ggplot(perc.religious.denom.year,aes(x=year,y=Percentage)) + 
  geom_area(aes(fill=religious)) +
  scale_fill_manual(values=Blues11) +
  facet_wrap(~denom,nrow=2) + theme_bw() + 
  labs(x="Year of fieldwork",
       fill="How religious")


ggsave("figures/percentages-religious-denomination-time-areas.pdf")
ggsave("figures/percentages-religious-denomination-time-areas.png")




