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

# ## 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.

ls()


ESS_work_west <- subset(ESS_relatt_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,
             relatt,
             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.treatment","contr.treatment"))


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)))
    EarlyLate <- factor()
    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 %<>% subset(relmemb == "Christian")


options(width=120)
ESS_work_west %$% levels(Country)


options(width=120)
ESS_work_west %$% table(Country)


ESS_work_west %<>% within({
    Stage <- recode(Country,
                           "Early" <- c("Denmark","Belgium","France","Netherlands","Norway","Sweden","United Kingdom"),
                           "Intermediate" <- c("Austria","Germany","Luxembourg","Switzerland"),
                           "Late" <- c("Ireland","Italy","Portugal","Spain")
                            )
})


mid <- function(x) mean(range(x,na.rm=TRUE))
std <- function(x) (x-mid(x))/abs(diff(range(x,na.rm=TRUE)))
revert_std <- function(z,x){
    y <- abs(diff(range(x,na.rm=TRUE)))*z
    y + mid(x)
}


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

# ## Model estimation

library(mclogit)


model_relprin <- mclogit(wgt|id~religious_principles.sq+
                                immigration.sq+
                                lrecon.sq+
                                religious_principles/(oesch8+relatt*Stage)*time+
                                immigration/(oesch8+relatt)*time+
                                lrecon/(oesch8+relatt)*time
                                -time,
                                random=~1|CHES_id/essround,
                                data=ESS_work_west,
                                 trace=TRUE,
                                estimator="REML"
                                )



save(model_relprin,file="churchat-stage-relprin/churchat-relprin-model.RData")


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

# ## Wald tests

source("Wald.R")


(Wald_relprin <- 
 Wald(model_relprin,
     list(~religious_principles,
          ~religious_principles:time,
          ~religious_principles:relatt,
          ~religious_principles:relatt:time,
          ~religious_principles:Stage,
          ~religious_principles:relatt:Stage,
          ~religious_principles:relatt:Stage:time
           )))


Wald_relprin %>%
        dimrename(          
        "religious_principles" = "Religious principles",
        "religious_principles:relatt" = "$\\times$ church attendance",
        "religious_principles:time" = "$\\times$ time",
        "religious_principles:relatt:time" = "$\\times$ time $\\times$ church attendance",
        "religious_principles:Stage" = "$\\times$ stage",
        "religious_principles:relatt:Stage" = "$\\times$ stage $\\times$ church attendance",
        "religious_principles:relatt:Stage:time" = "$\\times$ stage $\\times$ church attendance $\\times$ time"
        ) %>% rename(
            "Chisq"="$W$",
            "Pr(>Chisq)"="$p$-value") %>%
        toLatex(digits=c(1,0,3)) %>% writeLines(con="churchat-stage-relprin/Wald-stage-churchat.tex")

# ## Visualisation

predigrid.relprin <- with(subset(ESS_work_west,is.finite(religious_principles)),
                    expand.grid(
                            relatt=levels(relatt),
                            oesch8=levels(oesch8),
                            Stage=levels(Stage),
                            lrecon=0,
                            religious_principles=seq(from=-.5,to=.5,length=21),
                            immigration=0,
                            year=unique(year),
                            party="a"
                        ))


predigrid.relprin %<>% within({
    id <- interaction(relatt,religious_principles)
    id <- match(id,unique(id))
})  

predigrid.relprin %<>% within({
    lrecon.sq <- lrecon^2
    religious_principles.sq <- religious_principles^2
    immigration.sq <- immigration^2
    CHES_id <- 1
    essround <- 3
    country <- "Germany"
    time <- (year - 2010)/10
})    

predigrid.relprin %<>% sort(by=~id)


predigrid.relprin.aver <- subset(predigrid.relprin,
                                oesch8=="Service worker" &
                                year==2014 &
                                relatt %in% c("Never", "Holidays","Weekly")
                               )


pred.relprin.aver <- predict(model_relprin,
                               newdata=predigrid.relprin.aver,conditional=FALSE,
                               se.fit=TRUE)


predigrid.relprin.aver %<>% within({
        pred <- pred.relprin.aver$fit
        pred.upper <- pred.relprin.aver$fit + 2*pred.relprin.aver$se.fit
        pred.lower <- pred.relprin.aver$fit - 2*pred.relprin.aver$se.fit
        religious_principles <- 10*religious_principles + 5
})


library(ggplot2)


#!cell_opt: jupyter.plot.width=10, jupyter.plot.height=4.5
(g1 <- ggplot(predigrid.relprin.aver,
       aes(x=religious_principles)) +
    geom_vline(xintercept=7.5,color="gray80") +
    geom_line(aes(y=pred,linetype=relatt)) +
    geom_ribbon(aes(ymax=pred.upper,ymin=pred.lower,group=relatt),alpha=.1) +
    labs(x="Party religious principles position", y="Party utility",linetype=NULL) +
    coord_cartesian(ylim=c(-1.5,1)) +
    facet_wrap(~Stage) +
    theme_bw() +
    theme(legend.position="bottom"))


#!cell_opt: jupyter.plot.width=10, jupyter.plot.height=5
(g1a <- ggplot(predigrid.relprin.aver,
       aes(x=religious_principles)) +
    geom_vline(xintercept=7.5,color="gray80") +
    geom_line(aes(y=pred,linetype=Stage)) +
    geom_ribbon(aes(ymax=pred.upper,ymin=pred.lower,group=Stage),alpha=.1) +
    labs(x="Party religious principles position", y="Party utility",linetype="Stage:") +
    coord_cartesian(ylim=c(-1.5,1)) +
    facet_wrap(~relatt) +
    theme_bw() +
    theme(legend.position="bottom"))
ggsave(g1a,file="churchat-stage-relprin/churchat-stage.pdf")
ggsave(g1a,file="churchat-stage-relprin/churchat-stage.png")


predigrid.relprin.change.time <- subset(predigrid.relprin,
                                oesch8=="Service worker" &
                                religious_principles==.25 &
                                relatt %in% c("Never", "Holidays","Weekly"))


pred.relprin.change.time <- predict(model_relprin,
                               newdata=predigrid.relprin.change.time,conditional=FALSE,
                               se.fit=TRUE)
predigrid.relprin.change.time %<>% within({
        pred <- pred.relprin.change.time$fit
        pred.upper <- pred.relprin.change.time$fit + 2*pred.relprin.change.time$se.fit
        pred.lower <- pred.relprin.change.time$fit - 2*pred.relprin.change.time$se.fit
})


#!cell_opt: jupyter.plot.width=10, jupyter.plot.height=5
(g2 <- ggplot(predigrid.relprin.change.time,
       aes(x=year)) +
    geom_vline(xintercept=2014,color="gray80") +
    geom_line(aes(y=pred,linetype=relatt)) +
    geom_ribbon(aes(ymax=pred.upper,ymin=pred.lower,group=relatt),alpha=.1) +
    labs(x="Year of ESS wave", y="Party utility at pos = 7.5", linetype="Group:") +
    coord_cartesian(ylim=c(-1.5,1)) +
    facet_wrap(~Stage) +
    theme_bw() +
    theme(legend.position="bottom"))


#!cell_opt: jupyter.plot.width=10, jupyter.plot.height=4.5
(g2a <- ggplot(predigrid.relprin.change.time,
       aes(x=year)) +
    geom_vline(xintercept=2014,color="gray80") +
    geom_line(aes(y=pred,linetype=Stage)) +
    geom_ribbon(aes(ymax=pred.upper,ymin=pred.lower,group=Stage),alpha=.1) +
    labs(x="Year of ESS wave", y="Party utility at pos = 7.5", linetype="Stage:") +
    coord_cartesian(ylim=c(-1.5,1)) +
    facet_wrap(~relatt) +
    theme_bw() +
    theme(legend.position="bottom"))
ggsave(g2a,file="churchat-stage-relprin/churchat-stage-change.pdf")
ggsave(g2a,file="churchat-stage-relprin/churchat-stage-change.png")




