
suppressMessages(library(RKernel))
options(width=120)

# ## Preparing the data
# The data we are using are in the file "ESS-relig-CHES.RData".

suppressMessages(library(memisc))
library(magrittr)


load(file="ESS-relig-CHES.RData")


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


ESS_work_west <- subset(ESS_pray_denom_CHES,
    country2 %in% c(
        "Austria",
        "Belgium/Wallonia","Belgium/Flanders","Belgium/Brussels",
        "Denmark",
        "France",
        "Germany (West)","Germany/Bavaria",
        "Ireland",
        "Italy",
        "Luxembourg",
        "Netherlands",
        "Norway",
        "Portugal",
        "Spain/Castilia","Spain/Catalunia","Spain/Basque",
        "Sweden",
        "Switzerland",
        "United Kingdom/England",
        "United Kingdom/Scotland",
        "United Kingdom/Wales"
    ) & !is.na(galtan),
    select=c(CHES_id,
             year,
             essround,
             denom,
             pray,
             oesch8,
             wgt,
             id,
             country2,
             lrecon,
             galtan,
             spendvtax,
             deregulation,
             redistribution,
             civlib_laworder,
             sociallifestyle,
             religious_principles,
             immigration))


ESS_work_west %<>% as.data.frame()

# The information about religious denomination appears to be incorrect or misleading 
# in France until 2005 and in the United Kingdom from 2004 to 2007 and after 2020.

options(contrasts=c("contr.sum","contr.sum"))


ESS_work_west %<>% within({
    Country <- recode(as.character(country2),
                      c("Belgium/Wallonia","Belgium/Flanders","Belgium/Brussels") -> "Belgium",
                      c("Germany (West)","Germany/Bavaria") -> "Germany",
                      c("Spain/Castilia","Spain/Catalunia","Spain/Basque") -> "Spain",
                      c("United Kingdom/England","United Kingdom/Scotland","United Kingdom/Wales") 
                          -> "United Kingdom",
                      copy=TRUE)
    Country <- factor(Country,levels=sort(levels(Country)))
    denom[Country=="France" & year <= 2005 | 
         Country=="United Kingdom" & year >= 2004 & year <= 2007] <- NA
    cathprot <- recode(as.numeric(denom),
                      "Catholic"     <- 1,
                      "Protestant"   <- 2
                     )
})


ESS_work_west <- merge(ESS_work_west, perc.denom.country.wide,
                       by.x = 'Country',
                       by.y = 'country',
                       all.x = TRUE,
                       all.y = FALSE)


ESS_work_west %<>% within({
    foreach(x = c(lrecon,
                  galtan,
                  spendvtax,
                  deregulation,
                  redistribution,
                  civlib_laworder,
                  sociallifestyle,
                  religious_principles,
                  immigration),
            y = c(lrecon.sq,
                  galtan.sq,
                  spendvtax.sq,
                  deregulation.sq,
                  redistribution.sq,
                  civlib_laworder.sq,
                  sociallifestyle.sq,
                  religious_principles.sq,
                  immigration.sq),
            {
                x <- (x-5)/10
                y <- x^2
            })
    oesch8 <- factor(oesch8,levels=levels(oesch8)[c(3,1:2,4:8)]) # Production workers as baseline
    country <- factor(country2) # Drop countries not included in this subset
    time <- (year - 2010)/10
    contrasts(cathprot) <- contr("sum", base = 1)
    contrasts(oesch8) <- contr("sum", base = 1)
    contrasts(pray) <- contr("sum", base = 1)
    No_religion <- (No_religion - 50)/100
    Denom_Het <- 2*(1-(Cath_Prot/100)^2-(1-Cath_Prot/100)^2)
    Cath_Prot <- (Cath_Prot - 50)/100
    })


library(mclogit)


load("model-cathprot-prayer-galtan-context.RData")

# ## Visualisation

source("vis-helpers.R")


gridpoints <- sort(unique(c(-.25,.25,seq(from=-.5,to=.5,length=21))))


predigrid.galtan <- prob_grid_create(model=model_galtan,
                                      data=ESS_work_west,
                                      cathprot=levels(cathprot),
                                      pray=levels(pray),
                                      oesch8=levels(oesch8),
                                      lrecon=0,
                                      galtan=gridpoints,
                                      immigration=0,
                                      No_religion=0,
                                      Cath_Prot = c(95,50,5),
                                      reference=c(galtan=0),
                                      transformations=list(
                                          Denom_Het = 2*(1-(Cath_Prot/100)^2-(1-Cath_Prot/100)^2),
                                          Cath_Prot = (Cath_Prot-50)/100
                                      ))


predigrid.galtan %$$% {
    Cath_Prot <- Cath_Prot*100 + 50
    Cath_Prot <- recode(Cath_Prot,
                         "Mostly Catholic" <- 95,
                         "Mixed" <- 50,
                         "Mostly Protestant" <- 5) 
}


predigrid.galtan.churchat.aver <- subset(predigrid.galtan, time == 0 & oesch8 == "Clerks" & 
                                          pray %in% c("Never","Weekly","Daily") & Cath_Prot == "Mixed")
predigrid.galtan.churchat.time <- subset(predigrid.galtan, galtan == 7.5 & oesch8 == "Clerks" & 
                                          pray %in% c("Never","Weekly","Daily") & Cath_Prot == "Mixed")


predigrid.galtan.context.aver <- subset(predigrid.galtan, time == 0 & oesch8 == "Clerks" & 
                                          pray == "Daily")
predigrid.galtan.context.time <- subset(predigrid.galtan, galtan == 7.5 & oesch8 == "Clerks" & 
                                          pray == "Daily")


library(ggplot2)


add_pray <- function(string) paste("Prays:",string)
add_country <- function(string) paste("Country:",string)


#!cell_opt: jupyter.plot.width=10, jupyter.plot.height=4.5
(g.galtan.cathprot.aver <- ggcurves(data=predigrid.galtan.churchat.aver,
         x=galtan,
         y=prob,
         ymin=prob.lower,
         ymax=prob.upper,
         grp=cathprot,
         labels=c(x="GAL-TAN position",y="Predicted probability",linetype=""),
         vline=7.5,
         percentages=TRUE                   
         ) + facet_grid(~pray, labeller=labeller(.default=add_pray)))
ggsave(g.galtan.cathprot.aver,file="cathprot-context-galtan/cathprot-context-galtan-prayer.pdf")
ggsave(g.galtan.cathprot.aver,file="cathprot-context-galtan/cathprot-context-galtan-prayer.png")


#!cell_opt: jupyter.plot.width=4.5, jupyter.plot.height=4.5
(g.galtan.cathprot.aver <- ggcurves(data=subset(predigrid.galtan.churchat.aver, pray=="Daily"),
         x=galtan,
         y=prob,
         ymin=prob.lower,
         ymax=prob.upper,
         grp=cathprot,
         labels=c(x="GAL-TAN position",y="Predicted probability",linetype=""),
         vline=7.5,
         percentages=TRUE                   
         ))
ggsave(g.galtan.cathprot.aver,file="cathprot-context-galtan/cathprot-context-galtan-praydaily-aver.pdf")
ggsave(g.galtan.cathprot.aver,file="cathprot-context-galtan/cathprot-context-galtan-praydaily-aver.png")


#!cell_opt: jupyter.plot.width=10, jupyter.plot.height=4.5
(g.galtan.cathprot.time <- ggcurves(data=predigrid.galtan.churchat.time,
         x=year,
         y=prob,
         ymin=prob.lower,
         ymax=prob.upper,
         grp=cathprot,
         labels=c(x="GAL-TAN position",y="Predicted probability",linetype=""),
         vline=2014,
         percentages=TRUE                   
         ) + facet_grid(~pray, labeller=labeller(.default=add_pray)))
ggsave(g.galtan.cathprot.time,file="cathprot-context-galtan/cathprot-context-galtan-prayer-time.pdf")
ggsave(g.galtan.cathprot.time,file="cathprot-context-galtan/cathprot-context-galtan-prayer-time.png")


#!cell_opt: jupyter.plot.width=4.5, jupyter.plot.height=4.5
(g.galtan.cathprot.time <- ggcurves(data=subset(predigrid.galtan.churchat.time,pray=="Daily"),
         x=year,
         y=prob,
         ymin=prob.lower,
         ymax=prob.upper,
         grp=cathprot,
         labels=c(x="GAL-TAN position",y="Predicted probability",linetype=""),
         vline=2014,
         percentages=TRUE                   
         ))

ggsave(g.galtan.cathprot.time,file="cathprot-context-galtan/cathprot-context-galtan-praydaily-time.pdf")
ggsave(g.galtan.cathprot.time,file="cathprot-context-galtan/cathprot-context-galtan-praydaily-time.png")


#!cell_opt: jupyter.plot.width=8, jupyter.plot.height=5
(g.galtan.cathprot.aver <- ggcurves(data=predigrid.galtan.churchat.aver,
         x=galtan,
         y=prob,
         ymin=prob.lower,
         ymax=prob.upper,
         grp=pray,
         labels=c(x="GAL-TAN position",y="Predicted probability",linetype="Prays how often"),
         vline=7.5,
         percentages=TRUE
         ) + facet_grid(~cathprot))


#!cell_opt: jupyter.plot.width=10, jupyter.plot.height=4.5
(g.galtan.cathprot.context.aver <- ggcurves(data=predigrid.galtan.context.aver,
         x=galtan,
         y=prob,
         ymin=prob.lower,
         ymax=prob.upper,
         grp=cathprot,
         labels=c(x="GAL-TAN position",y="Predicted probability",linetype=""),
         vline=7.5,
         percentages=TRUE
         ) + facet_grid(~Cath_Prot,labeller=labeller(.default=add_country)))
ggsave(g.galtan.cathprot.context.aver,file="cathprot-context-galtan/cathprot-context-galtan-context-prayer.pdf")
ggsave(g.galtan.cathprot.context.aver,file="cathprot-context-galtan/cathprot-context-galtan-context-prayer.png")


#!cell_opt: jupyter.plot.width=10, jupyter.plot.height=4.5
(g.galtan.cathprot.context.time <- ggcurves(data=predigrid.galtan.context.time,
         x=year,
         y=prob,
         ymin=prob.lower,
         ymax=prob.upper,
         grp=cathprot,
         labels=c(x="GAL-TAN position",y="Predicted probability",linetype=""),
         vline=2014,
         percentages=TRUE
         ) + facet_grid(~Cath_Prot,labeller=labeller(.default=add_country)))
ggsave(g.galtan.cathprot.time,file="cathprot-context-galtan/cathprot-context-galtan-context-prayer-time.pdf")
ggsave(g.galtan.cathprot.time,file="cathprot-context-galtan/cathprot-context-galtan-context-prayer-time.png")


#!cell_opt: jupyter.plot.width=8, jupyter.plot.height=5
(g.galtan.cathprot.context.aver <- ggcurves(data=predigrid.galtan.context.aver,
         x=galtan,
         y=prob,
         ymin=prob.lower,
         ymax=prob.upper,
         grp=Cath_Prot,
         labels=c(x="GAL-TAN position",y="Predicted probability",linetype="Country:"),
         vline=7.5,
         percentages=TRUE
         ) + facet_grid(~cathprot))




