
suppressMessages(library(RKernel))

# ## 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")

# We restrict our analysis to the denominationally mixed countries in Western
# Europe, that is, Germany, the Netherlands, Switzerland, and the United Kingdom.

ESS_work_west <- subset(ESS_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"
    ),
    select=c(CHES_id,
             year,
             essround,
             denom,
             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
    relmemb <- recode(as.numeric(denom),
                      "Christian"     <- 1:4,
                      "Non-Christian" <- 5,
                      "No religion"   <- 6
                     )
})


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(relmemb) <- contr("sum", base = 1)
    contrasts(oesch8) <- contr("sum", base = 1)
    })

# ## Loading the models

library(mclogit)


load("model-relmemb-relprin-context.RData")


source("Wald.R")


(Wald_relprin_relmemb_class <- 
 Wald(model_relprin,
     list(~religious_principles:time,
          ~religious_principles:relmemb,
          ~religious_principles:relmemb:time,
          ~religious_principles:oesch8,
          ~religious_principles:oesch8:time
           )))

# ## Visualisation

source("vis-helpers.R")

# ### Comparing religious membership and class regarding religious principles positions

probgrid.relprin <- prob_grid_create(model=model_relprin,
                                      data=subset(ESS_work_west,is.finite(religious_principles)),
                                      relmemb=levels(relmemb),
                                      oesch8=levels(oesch8),
                                      lrecon=0,
                                      religious_principles=.25,
                                      immigration=0,
                                      Secularization=0,
                                      center_time=2014,
                                      reference=c(religious_principles=0))
probs.relprin.relmemb <- subset(probgrid.relprin, oesch8 == "Clerks" & year == 2014)
probs.relprin.oesch8 <- subset(probgrid.relprin, relmemb == "Christian" & year == 2014)



library(ggplot2)


library(gridExtra)


(g.probs.relprin.relmemb <- ggpoints(data=probs.relprin.relmemb,
         x=relmemb,
         y=prob,
         ymin=prob.lower,
         ymax=prob.upper,
         labels=c(title="Religious membership",y="Predicted probability",x=""),
         percentages=TRUE,
         ylim=c(.25,.75)
         ))


(g.probs.relprin.oesch8 <- ggpoints(data=probs.relprin.oesch8,
         x=oesch8,
         y=prob,
         ymin=prob.lower,
         ymax=prob.upper,
         labels=c(title="Occupational class",y="Predicted probability",x=""),
         percentages=TRUE,
         ylim=c(.25,.75)
         ))


#!cell_opt: jupyter.plot.width=8,jupyter.plot.height=4
(gg <- grid.arrange(g.probs.relprin.relmemb + labs(title="(Non-)religious group"),
             drop_annot_y(g.probs.relprin.oesch8) + labs(title="Occupational class"),
             nrow=1,widths=c(4,8)))
ggsave(gg,file="relmemb-context-comparison/compare-relmemb-class-relprin.pdf")
ggsave(gg,file="relmemb-context-comparison/compare-relmemb-class-relprin.png")

# ## Visualization: Any divergence between classes?

probs.relprin.relmemb.time <- subset(probgrid.relprin, oesch8 == "Clerks")
probs.relprin.oesch8.time <- subset(probgrid.relprin, relmemb == "Christian")


#!cell_opt: jupyter.plot.height=4.5, jupyter.plot.width=4.5
(g.relprin.relmemb.time <- ggcurves(data=probs.relprin.relmemb.time,
         x=year,
         y=prob,
         ymin=prob.lower,
         ymax=prob.upper,
         grp=relmemb,
         labels=c(x="Year of ESS fieldwork",y="Predicted probability",linetype=""),
         vline=2014,
         percentages=TRUE
         ))


#!cell_opt: jupyter.plot.height=4.5, jupyter.plot.width=4.5
(g.relprin.oesch8.time <- ggcurves(data=probs.relprin.oesch8.time,
         x=year,
         y=prob,
         ymin=prob.lower,
         ymax=prob.upper,
         grp=oesch8,
         labels=c(x="Year of ESS fieldwork",y="Predicted probability",linetype=""),
         vline=2014,
         percentages=TRUE
         ))

# ## Visualisation: Comparing religious membership and class regarding economic left-right positions

probgrid.lrecon <- prob_grid_create(model=model_relprin,
                                      data=subset(ESS_work_west,is.finite(religious_principles)),
                                      relmemb=levels(relmemb),
                                      oesch8=levels(oesch8),
                                      lrecon=.25,
                                      religious_principles=0,
                                      immigration=0,
                                      Secularization=0,
                                      year=2014,
                                      center_time=2014,
                                      reference=c(lrecon=0))
probs.lrecon.relmemb <- subset(probgrid.lrecon, oesch8 == "Clerks")
probs.lrecon.oesch8 <- subset(probgrid.lrecon, relmemb == "Christian")


library(ggplot2)


library(gridExtra)


(g.probs.lrecon.relmemb <- ggpoints(data=probs.lrecon.relmemb,
         x=relmemb,
         y=prob,
         ymin=prob.lower,
         ymax=prob.upper,
         labels=c(title="Religious membership",y="Predicted probability",x=""),
         percentages=TRUE,
         ylim=c(.25,.75)
         ))


(g.probs.lrecon.oesch8 <- ggpoints(data=probs.lrecon.oesch8,
         x=oesch8,
         y=prob,
         ymin=prob.lower,
         ymax=prob.upper,
         labels=c(title="Occupational class",y="Predicted probability",x=""),
         percentages=TRUE,
         ylim=c(.25,.75)
         ))


#!cell_opt: jupyter.plot.width=8,jupyter.plot.height=4
(gg <- grid.arrange(g.probs.lrecon.relmemb + labs(title="(Non-)religious group"),
             drop_annot_y(g.probs.lrecon.oesch8) + labs(title="Occupational class"),
             nrow=1,widths=c(4,8)))
ggsave(gg,file="relmemb-context-comparison/compare-relmemb-class-lrecon.pdf")
ggsave(gg,file="relmemb-context-comparison/compare-relmemb-class-lrecon.png")




# ## Visualisation: Comparing religious membership and class regarding immigration positions

probgrid.immigration <- prob_grid_create(model=model_relprin,
                                      data=subset(ESS_work_west,is.finite(religious_principles)),
                                      relmemb=levels(relmemb),
                                      oesch8=levels(oesch8),
                                      lrecon=0,
                                      religious_principles=0,
                                      immigration=.25,
                                      Secularization=0,
                                      year=2014,
                                      center_time=2014,
                                      reference=c(immigration=0))
probs.immigration.relmemb <- subset(probgrid.immigration, oesch8 == "Clerks")
probs.immigration.oesch8 <- subset(probgrid.immigration, relmemb == "Christian")


library(ggplot2)


library(gridExtra)


(g.probs.immigration.relmemb <- ggpoints(data=probs.immigration.relmemb,
         x=relmemb,
         y=prob,
         ymin=prob.lower,
         ymax=prob.upper,
         labels=c(title="Religious membership",y="Predicted probability",x=""),
         percentages=TRUE,
         ylim=c(.25,.75)
         ))


(g.probs.immigration.oesch8 <- ggpoints(data=probs.immigration.oesch8,
         x=oesch8,
         y=prob,
         ymin=prob.lower,
         ymax=prob.upper,
         labels=c(title="Occupational class",y="Predicted probability",x=""),
         percentages=TRUE,
         ylim=c(.25,.75)
         ))


#!cell_opt: jupyter.plot.width=8,jupyter.plot.height=4
(gg <- grid.arrange(g.probs.immigration.relmemb + labs(title="(Non-)religious group"),
             drop_annot_y(g.probs.immigration.oesch8) + labs(title="Occupational class"),
             nrow=1,widths=c(4,8)))
ggsave(gg,file="relmemb-context-comparison/compare-relmemb-class-immigration.pdf")
ggsave(gg,file="relmemb-context-comparison/compare-relmemb-class-immigration.png")










