###############################################################################
# Energie-Enquete
#
# What makes incentive-based instruments unpopular
# Isabelle Stadelmann-Steffen & Clau Dermont
# University of Bern
#
# clau.dermont@protonmail.ch
# as of December 2017
###############################################################################

# based on Conjoint-Dataset cjIncentives.RData
# estimate AMCEs, interaction effects and produce figures
# using cjoint 2.0.6

###############################################################################

if (!require("pacman", quietly=T)) install.packages("pacman")
#### Setup
pacman::p_load('tidyverse','cjoint','lme4','likert',
               'rtf')

setwd("~/Git/paper_incentives")
cjPolicy <- read.table('cjPolicy.csv', fileEncoding = 'UTF-8', header=T, 
                       sep=';', stringsAsFactors = F)
taxreform <- read.table('taxreform.csv', fileEncoding = 'UTF-8', header=T, 
                     sep=';')

# custom function
source("https://raw.githubusercontent.com/cdermont/r_cjoint_cd/master/plotamce.R")


###############################################################################
# Conjoint attributes labeling

cjPolicy$'Source.of.Funding' <- recode_factor(
  cjPolicy$'Source.of.Funding',
  "3" = "tax on electricity",
  "4" = "tax on non-renewable electricity",
  "1" = "income and revenue tax",
  "2" = "value added tax"
)

cjPolicy$'Policy.Measure' <- recode_factor(
  cjPolicy$'Policy.Measure',
  "1" = "investment grants",
  "2" = "feed-in tariff",
  "3" = "tax release",
  "4" = "redistribution"
)

cjPolicy$'Energy.Source.Priority' <-
  recode_factor(
    cjPolicy$'Energy.Source.Priority',
    "1" = "renewable energy in general",
    "2" = "solar power",
    "3" = "wind power",
    "4" = "small-scale hydro power",
    "5" = "geothermal power",
    "6" = "no priority"
  )

cjPolicy$'Costs' <- 
  recode_factor(
    cjPolicy$'Costs',
    "1" = "none",
    "2" = "+8 CHF",
    "3" = "+15 CHF",
    "4" = "+23 CHF",
    "5" = "+30 CHF"
  )

cjPolicy$'Exceptions' <- 
  recode_factor(
    cjPolicy$'Exceptions',
    "1" = "no exceptions",
    "2" = "industry"
  )

cjPolicy$'Running.Time' <- 
  recode_factor(
    cjPolicy$'Running.Time',
    "1" = "10 years",
    "2" = "20 years",
    "3" = "35 years"
  )

cjPolicy$'Nuclear.Power.Plants' <-
  recode_factor(
    cjPolicy$'Nuclear.Power.Plants',
    "1" = "switch off",
    "2" = "60 years run-time limit",
    "3" = "no run-time limit"
  )

# Conjoint Design

cjPolAttr <- list()
cjPolAttr$'Source.of.Funding' <- levels(cjPolicy$'Source.of.Funding')
cjPolAttr$'Policy.Measure' <- levels(cjPolicy$'Policy.Measure')
cjPolAttr$'Costs' <- levels(cjPolicy$'Costs')
cjPolAttr$'Exceptions' <- levels(cjPolicy$'Exceptions')
cjPolAttr$'Running.Time' <- levels(cjPolicy$'Running.Time')
cjPolAttr$'Energy.Source.Priority' <- levels(cjPolicy$'Energy.Source.Priority')
cjPolAttr$'Nuclear.Power.Plants' <- levels(cjPolicy$'Nuclear.Power.Plants')

cjPolConst <- list()
cjPolConst[[1]] <- list()
cjPolConst[[1]][["Policy.Measure"]] <- c("redistribution")
cjPolConst[[1]][["Source.of.Funding"]] <- c("income and revenue tax", 
                                            "value added tax")
cjPolConst[[2]] <- list()
cjPolConst[[2]][["Exceptions"]] <- c("industry")
cjPolConst[[2]][["Source.of.Funding"]] <- c("income and revenue tax", 
                                            "value added tax")

cjPolMargW <- list()
cjPolMargW$"Source.of.Funding" <- c(1/4, 1/4, 1/4, 1/4)
cjPolMargW$"Policy.Measure" <- c(27/100, 27/100, 26/100, 20/100)
cjPolMargW$"Costs" <- c(1/5, 1/5, 1/5, 1/5, 1/5)
cjPolMargW$"Exceptions" <- c(1/4, 3/4)
cjPolMargW$"Running.Time" <- c(1/3, 1/3, 1/3)
cjPolMargW$"Energy.Source.Priority" <- c(26/100, 16/100, 16/100, 16/100, 
                                         16/100, 10/100)
cjPolMargW$"Nuclear.Power.Plants" <- c(1/3, 1/3, 1/3)

cjPolDesign <- makeDesign(type='constraints', attribute.levels=cjPolAttr, 
                          constraints=cjPolConst, level.probs=cjPolMargW)
rm(cjPolAttr, cjPolConst, cjPolMargW)

cjPolBase <- list()
cjPolBase$"Source.of.Funding" <- "tax on electricity"
cjPolBase$"Policy.Measure" <- "feed-in tariff"
cjPolBase$"Costs" <- "none"
cjPolBase$"Exceptions" <- "no exceptions"
cjPolBase$"Running.Time" <- "10 years"
cjPolBase$"Energy.Source.Priority" <- "no priority"
cjPolBase$"Nuclear.Power.Plants" <- "no run-time limit"

# Taxreform labeling

for (i in 3:7) {
  taxreform[,i] <- recode_factor(
    taxreform[,i], 
    "1"="Disagree",
    "2"="Rather disagree",
    "5"="Don't know",
    "3"="Rather agree",
    "4"="Agree"
  )
}

###############################################################################
# Figures and Tables

###### Support for concepts
# Output: Fig. 1

cjPolicy$concat <- paste(cjPolicy$Source.of.Funding,cjPolicy$Policy.Measure,
                         cjPolicy$Costs,cjPolicy$Exceptions,
                         cjPolicy$Running.Time,cjPolicy$Nuclear.Power.Plants,
                         cjPolicy$Energy.Source.Priority,sep="_")

policySupport <- cjPolicy %>% 
  group_by(concat) %>% 
  summarise(mean=mean(support-1)/10) # scale is 1-11

policySupport <- policySupport %>% 
  separate(concat,
           c("Source.of.Funding","Policy.Measure","Costs","Exceptions",
             "Running.Time","Nuclear.Power.Plants","Energy.Source.Priority"),
           sep="_", remove=T) %>%
  bind_cols(policySupport)

policySupport$colorize <- "no"
policySupport$colorize[
  policySupport$Source.of.Funding=="tax on electricity" &
    policySupport$Policy.Measure=="redistribution"] <- "yes"

ggplot(
  policySupport,aes(mean)) + 
  geom_density(fill="grey80") + 
  geom_vline(aes(xintercept=0.5), color="grey30", linetype=5) +
  scale_fill_manual(values=c("#2b8cbe","#ece7f2")) +
  scale_linetype_manual(values=c(1,2)) +
  scale_x_continuous(labels=scales::percent_format()) +
  labs(title="Support for Electricity Policy Proposals",
       x="Mean support per concept",
       y="Density") + 
  theme_minimal() + 
  theme(legend.position="none")

ggsave("img/fig1.eps",
       dpi=300,
       units=c("in"),
       width=9, height=4)
ggsave("img/fig1.png",
       dpi=300,
       units=c("in"),
       width=9, height=4)


###### Model 1: general Conjoint / AMCE
# Output: Fig. 2

m1 <- amce(support ~ Source.of.Funding + Policy.Measure + 
             Costs + Exceptions + Running.Time + 
             Energy.Source.Priority +  Nuclear.Power.Plants,  
           data=cjPolicy, 
           cluster=TRUE, respondent.id="id", 
           baselines=cjPolBase, design=cjPolDesign)

p.m1 <- plot.amce2(m1)

# extract data from figure for further refining
p.m1.d <- p.m1$plot$data
p.m1.d$x <- seq(1,nrow(p.m1.d), 1)

p.m1.d$x2 <- c(11,12,13,14,15,16,24,25,26,27,28,
               29,30,17,18,19,31,32,33,34,6,7,8,
               9,10,20,21,22,23,1,2,3,4,5)
p.m1.d <- filter(p.m1.d, x2<17)

ggplot(p.m1.d) +
  geom_hline(yintercept = 0, size = 0.5, colour = "black", 
             linetype = "dotted") +
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper,
                      x = reorder(printvar, -x2)),
                  size=0.4) +
  ylim(-0.76,0.25) + coord_flip() + 
  labs(title="Policy Support", 
       y="Change in Support", 
       x="") +
  theme_minimal() +
  theme(axis.text.y = element_text(hjust = 0), 
        axis.title.x=element_text(size=12), 
        legend.position="bottom")

ggsave("img/fig2.eps",
       dpi=300,
       units=c("in"),
       width=7, height=4)
ggsave("img/fig2.png",
       dpi=300,
       units=c("in"),
       width=7, height=4)
rm(p.m1, p.m1.d)


###### Understanding tax reforms
# Output: Fig. 3

str(taxreform)

steer <- dplyr::select(taxreform, starts_with("es"))
colnames(steer) <- c("Tax reduces consumption",
                     "Redistribution doesn't reduce consumption",
                     "Low consumption means paying less",
                     "Redistribution puts a strain on public finances",
                     "There is double dividend")

steer <- likert(steer)

likert.bar.plot(steer, 
                          low.color="#ca0020", high.color="#0571b0",
                          legend="Response") + 
  labs(title="Opinion on steering tax mechanism",
       x="",
       y="Share of responses in %") + 
  theme_minimal() + 
  theme(axis.text.y = element_text(hjust = 0), 
        legend.position="bottom")

ggsave("img/fig3.eps",
       dpi=300,
       units=c("in"),
       width=9, height=4)
ggsave("img/fig3.png",
       dpi=300,
       units=c("in"),
       width=9, height=4)
rm(steer)


###### Model 2: include political predisposition
# Output: Fig. 4
cjPolicy$lirecat <- factor(cjPolicy$lirecat, 
                              levels=c("left","middle","right"))
cjPolicy.na <- filter(cjPolicy, !is.na(lirecat))

m2 <- amce(support ~ Source.of.Funding:lirecat + Policy.Measure:lirecat +
             Costs:lirecat + Exceptions:lirecat + Running.Time:lirecat +
             Energy.Source.Priority:lirecat + Nuclear.Power.Plants:lirecat, 
           data=cjPolicy.na, respondent.varying="lirecat",
           cluster=TRUE, respondent.id="id", 
           baselines=cjPolBase, design=cjPolDesign)

p.m2 <- plot.amce2(m2, plot.display="interaction")

p.m2.d <- p.m2$plot$data
p.m2.d$x <- seq(1,nrow(p.m2.d), 1)
p.m2.d$x2 <- seq(1,nrow(p.m2.d)/3, 1)
p.m2.d$facet <- recode_factor(p.m2.d$facet, 
                            'Conditional on\nlirecat = left' = "left",
                            'Conditional on\nlirecat = middle' = "middle", 
                            'Conditional on\nlirecat = right' = "right"
)
p.m2.d <- filter(p.m2.d, x2<17)

ggplot(p.m2.d) +
  geom_hline(yintercept = 0, size = 0.5, colour = "black", 
             linetype = "dotted") +
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper,
                      x = reorder(printvar, -x), color=facet, 
                      shape=facet),
                  size=0.4,
                  position=position_dodge(width=0.95)) +
  ylim(-1,0.75) + coord_flip() + 
  facet_grid(~facet) +
  labs(title="Policy Support by political ideology", 
       y="Change in Support", 
       x="", 
       color="Political predisposition",
       shape="Political predisposition") +
  scale_color_brewer(palette="Set1") + 
  theme_minimal() +
  theme(axis.text.y = element_text(hjust = 0), 
        axis.title.x=element_text(size=12), 
        legend.position="bottom")

ggsave("img/fig4.eps",
       dpi=300,
       units=c("in"),
       width=9, height=4)
ggsave("img/fig4.png",
       dpi=300,
       units=c("in"),
       width=9, height=4)
rm(p.m2, p.m2.d)


###### Support for concepts by political predisposition
# Output: Fig. 5

policySupport <- filter(cjPolicy, !is.na(lirecat)) %>% 
  group_by(lirecat) %>% 
  summarise(groupmean=mean(support-1)/10)

policySupport <- filter(cjPolicy, !is.na(lirecat)) %>% 
  group_by(lirecat,concat) %>% 
  summarise(mean=mean(support-1)/10) %>%
  full_join(policySupport, by="lirecat")

ggplot(
  policySupport,aes(x=mean, fill=as.factor(lirecat),
                    linetype=as.factor(lirecat))) + 
  geom_density() + 
  geom_vline(aes(xintercept=groupmean), color="grey20", linetype=5) +
  facet_grid(lirecat~.) +
  scale_fill_brewer(palette="Pastel1") + 
  scale_x_continuous(labels=scales::percent_format()) +
  labs(title="Support for Electricity Policy Proposal",
       x="Mean support per concept (line: group mean)",
       y="Density",
       fill="Political predisposition",
       linetype="Political predisposition") + 
  theme_minimal() + 
  theme(legend.position="none")

ggsave("img/fig5.eps",
       dpi=300,
       units=c("in"),
       width=9, height=5)
ggsave("img/fig5.png",
       dpi=300,
       units=c("in"),
       width=9, height=5)
rm(policySupport)


###### Model results from M1
# Output: Tab. A3

res.m1.amce <- summary(m1)$amce[,c(1:4,7)]

res.m1 <- select(res.m1.amce, Attribute, Level)
res.m1$nr <- seq(1,nrow(res.m1),1)
res.m1.amce <- full_join(res.m1, res.m1.amce, by=c("Attribute","Level"))
res.m1.amce$m <- 1

res.m1 <- rbind(res.m1.amce)
colnames(res.m1) <- c("Attribute","Level","nr","Estimate","SE","star","m")
res.m1 <- gather(res.m1, "var", "val", 4:5)
res.m1 <- arrange(res.m1, m, nr)

res.m1$nrx <- seq(1,2,1)
res.m1$p[res.m1$nrx==1] <- paste0(round(res.m1$val[res.m1$nrx==1], 2), 
                                  res.m1$star[res.m1$nrx==1])
res.m1$p[res.m1$nrx==2] <- paste0("(", round(res.m1$val[res.m1$nrx==2], 2), ")")

res.m1 <- select(res.m1, Attribute, Level, p, m, nr, nrx)
res.m1 <- spread(res.m1, m, p)
res.m1 <- arrange(res.m1, nr, nrx)
res.m1$Level <- paste0("    ",res.m1$Level)
res.m1$Level[res.m1$nrx==2] <- ""

res.m1.attr <- as.data.frame(cbind(
  c(0.5, 4.5, 9.5, 10.5, 12.5, 15.5, 17.5),
  c(
    "Costs (Baseline = none)",
    "Energy Source Priority (Baseline = no priority)",
    "Exceptions (Baseline = no exceptions)",
    "Nuclear Power Plants (Baseline = switch off)",
    "Policy Measure (Baseline = feed-in tariff)",
    "Running Time (Baseline = 10 years)",
    "Source of Funding (Baseline = tax on electricity)"
  )
)
)
colnames(res.m1.attr) <- c("nr","Level")
res.m1.attr$nr <- as.numeric(as.character(res.m1.attr$nr))
res.m1 <- full_join(res.m1, res.m1.attr, by=c("Level","nr"))

res.m1 <- arrange(res.m1, nr, nrx)
res.m1 <- select(res.m1, -Attribute, -nr, -nrx)
colnames(res.m1) <- c("","AMCE")

rtffile <- RTF("m1.doc")
addParagraph(rtffile, "Table A.3: Full results of model 1\n")
addTable(rtffile, as.data.frame(res.m1))
done(rtffile)
rm(res.m1, res.m1.attr, res.m1.amce)

###### Model results from M2
# Output: Tab. A5

res.m2.amce <- summary(m2)$amce[,c(1:4,7)]
res.m2.cond0 <- summary(m2)$lirecat1amce[c(1:4,7)]
res.m2.cond1 <- summary(m2)$lirecat2amce[c(1:4,7)]
res.m2.cond2 <- summary(m2)$lirecat3amce[c(1:4,7)]

res.m2 <- select(res.m2.amce, Attribute, Level)
res.m2$nr <- seq(1,nrow(res.m2),1)
res.m2.amce <- full_join(res.m2, res.m2.amce, by=c("Attribute","Level"))
res.m2.cond0 <- full_join(res.m2, res.m2.cond0, by=c("Attribute","Level"))
res.m2.cond1 <- full_join(res.m2, res.m2.cond1, by=c("Attribute","Level"))
res.m2.cond2 <- full_join(res.m2, res.m2.cond2, by=c("Attribute","Level"))
res.m2.amce$m <- 1
res.m2.cond0$m <- 2
res.m2.cond1$m <- 3
res.m2.cond2$m <- 4

res.m2 <- rbind(res.m2.amce, res.m2.cond0, res.m2.cond1, res.m2.cond2)
colnames(res.m2) <- c("Attribute","Level","nr","Estimate","SE","star","m")
res.m2 <- gather(res.m2, "var", "val", 4:5)
res.m2 <- arrange(res.m2, m, nr)

res.m2$nrx <- seq(1,2,1)
res.m2$p[res.m2$nrx==1] <- paste0(round(res.m2$val[res.m2$nrx==1], 2), 
                                    res.m2$star[res.m2$nrx==1])
res.m2$p[res.m2$nrx==2] <- paste0("(", round(res.m2$val[res.m2$nrx==2], 2), ")")

res.m2 <- select(res.m2, Attribute, Level, p, m, nr, nrx)
res.m2 <- spread(res.m2, m, p)
res.m2 <- arrange(res.m2, nr, nrx)
res.m2$Level <- paste0("    ",res.m2$Level)
res.m2$Level[res.m2$nrx==2] <- ""

res.m2.attr <- as.data.frame(cbind(
  c(0.5, 4.5, 9.5, 10.5, 12.5, 15.5, 17.5),
  c(
    "Costs (Baseline = none)",
    "Energy Source Priority (Baseline = no priority)",
    "Exceptions (Baseline = no exceptions)",
    "Nuclear Power Plants (Baseline = switch off)",
    "Policy Measure (Baseline = feed-in tariff)",
    "Running Time (Baseline = 10 years)",
    "Source of Funding (Baseline = tax on electricity)"
  )
)
)
colnames(res.m2.attr) <- c("nr","Level")
res.m2.attr$nr <- as.numeric(as.character(res.m2.attr$nr))
res.m2 <- full_join(res.m2, res.m2.attr, by=c("Level","nr"))

res.m2 <- arrange(res.m2, nr, nrx)
res.m2 <- select(res.m2, -Attribute, -nr, -nrx)
colnames(res.m2) <- c("","AMCE","left","middle","right")

rtffile <- RTF("m2.doc")
addParagraph(rtffile, "Table A.5: Full results of model 2\n")
addTable(rtffile, as.data.frame(res.m2))
done(rtffile)
rm(res.m2, res.m2.attr, res.m2.amce)


###### Understanding tax reforms by political predisposition
# Output: Fig. A1

taxreform$lirecat <- NA
taxreform$lirecat[taxreform$lire>=0] <- "left"
taxreform$lirecat[taxreform$lire>=4] <- "middle"
taxreform$lirecat[taxreform$lire>=7] <- "right"
taxreform$lirecat <- factor(taxreform$lirecat, 
                            levels=c("left","middle","right"))

steer.na <- na.omit(taxreform)
steer <- dplyr::select(steer.na, starts_with("es"))
colnames(steer) <- c("Tax reduces consumption",
                     "Redistribution doesn't reduce consumption",
                     "Low consumption means paying less",
                     "Redistribution puts a strain on public finances",
                     "There is double dividend")

steer <- likert(steer, grouping=steer.na$lirecat)

likert.bar.plot(steer, low.color="#ca0020", high.color="#0571b0",
                legend="Response") + 
  labs(title="Opinion on steering tax mechanism",
       x="",
       y="Share of responses in %") + 
  theme_minimal() + 
  theme(axis.text.y = element_text(hjust = 0), 
        legend.position="bottom")

ggsave("img/figA2.eps",
       dpi=300,
       units=c("in"),
       width=9, height=6)
ggsave("img/figA2.png",
       dpi=300,
       units=c("in"),
       width=9, height=6)
rm(steer, steer.na)


###############################################################################

###### Further analysis of support

# cutoff-model
cjPolicy$supportx <- as.numeric(cjPolicy$support>=8)

m3 <- amce(supportx ~ Source.of.Funding + Policy.Measure + 
             Costs + Exceptions + Running.Time + 
             Energy.Source.Priority +  Nuclear.Power.Plants,  
           data = cjPolicy, 
           cluster = TRUE, respondent.id = "id", 
           baselines = cjPolBase, design = cjPolDesign)

p.m3 <- plot.amce2(m3)

# extract data from figure for further refining
p.m3.d <- p.m3$plot$data
p.m3.d$x <- seq(1,nrow(p.m3.d), 1)

p.m3.d$x2 <- c(11,12,13,14,15,16,24,25,26,27,28,
               29,30,17,18,19,31,32,33,34,6,7,8,
               9,10,20,21,22,23,1,2,3,4,5)
p.m3.d <- filter(p.m3.d, x2<17)

ggplot(p.m3.d) +
  geom_hline(yintercept = 0, size = 0.5, colour = "black", 
             linetype = "dotted") +
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper,
                      x = reorder(printvar, -x2)),
                  size=0.4) +
  ylim(-0.1,0.05) + coord_flip() + 
  labs(title="Policy Support (>=8 as 'yes')", 
       y="Change in Support", 
       x="") +
  theme_minimal() +
  theme(axis.text.y = element_text(hjust = 0), 
        axis.title.x=element_text(size=12), 
        legend.position="bottom")

ggsave("img/figA3.eps",
       dpi=300,
       units=c("in"),
       width=7, height=4)
ggsave("img/figA3.png",
       dpi=300,
       units=c("in"),
       width=7, height=4)


###### Marginal effects for subgroups
cjPolicy$female <- recode_factor(cjPolicy$female, 
                          '0' = 'male',
                          '1' = 'female')
cjPolicy$agecat <- factor(cjPolicy$agecat, 
                          levels=c("18-35 years","36-50 years","51-65 years",
                                   "66 years and older"))
cjPolicy$educat <- factor(cjPolicy$educat, 
                      levels=c("low","middle","high"))
cjPolicy$incomecat <- factor(cjPolicy$incomecat, 
                         levels=c("low","middle","high"))

## Model with gender
cjPolicy$xfactor <- cjPolicy$female
cjPolicy.na <- filter(cjPolicy, !is.na(xfactor))

m3.sex <- amce(support ~ Source.of.Funding:xfactor + Policy.Measure:xfactor + 
                 Costs:xfactor + Exceptions:xfactor + Running.Time:xfactor +
                 Energy.Source.Priority:xfactor + Nuclear.Power.Plants:xfactor,
               data = cjPolicy.na,
               respondent.id = "id",
               respondent.varying = "xfactor",
               baselines = cjPolBase, design = cjPolDesign)

summary(m3.sex)

p.m3 <- plot.amce2(m3.sex, plot.display="interaction")

p.m3.d <- p.m3$plot$data
p.m3.d$x <- seq(1,nrow(p.m3.d), 1)
p.m3.d$x2 <- seq(1,nrow(p.m3.d)/2, 1)
p.m3.d$facet <- recode_factor(p.m3.d$facet, 
                              'Conditional on\nxfactor = 0' = "male",
                              'Conditional on\nxfactor = 1' = "female"
)

ggplot(filter(p.m3.d, x2 < 17)) +
  geom_hline(yintercept = 0, size = 0.5, colour = "black", 
             linetype = "dotted") +
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper,
                      x = reorder(printvar, -x), color=facet, 
                      shape=facet),
                  size=0.4,
                  position=position_dodge(width=0.95)) +
  ylim(-1,0.5) + coord_flip() + 
  facet_grid(~facet) +
  labs(title="Policy Support by gender", 
       y="Change in Support", 
       x="") +
  scale_color_brewer(palette="Set1") + 
  theme_minimal() +
  theme(axis.text.y = element_text(hjust = 0), 
        axis.title.x=element_text(size=12), 
        legend.position="none")


## Model with age
cjPolicy$xfactor <- cjPolicy$agecat
cjPolicy.na <- filter(cjPolicy, !is.na(xfactor))

m3.age <- amce(support ~ Source.of.Funding:xfactor + Policy.Measure:xfactor + 
                 Costs:xfactor + Exceptions:xfactor + Running.Time:xfactor +
                 Energy.Source.Priority:xfactor + Nuclear.Power.Plants:xfactor,
               data = cjPolicy.na,
               respondent.id = "id",
               respondent.varying = "xfactor",
               baselines = cjPolBase, design = cjPolDesign)

summary(m3.age)

p.m3 <- plot.amce2(m3.age, plot.display="interaction")

p.m3.d <- p.m3$plot$data
p.m3.d$x <- seq(1,nrow(p.m3.d), 1)
p.m3.d$x2 <- seq(1,nrow(p.m3.d)/4, 1)
p.m3.d$facet <- recode_factor(p.m3.d$facet, 
                              'Conditional on\nxfactor = 18-35 years' = "18-35 years",
                              'Conditional on\nxfactor = 36-50 years' = "36-50 years",
                              'Conditional on\nxfactor = 51-65 years' = "51-65 years",
                              'Conditional on\nxfactor = 66 years and older' = "66 years and older"
)

ggplot(filter(p.m3.d, x2 < 17)) +
  geom_hline(yintercept = 0, size = 0.5, colour = "black", 
             linetype = "dotted") +
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper,
                      x = reorder(printvar, -x), color=facet, 
                      shape=facet),
                  size=0.4,
                  position=position_dodge(width=0.95)) +
  ylim(-1,0.8) + coord_flip() + 
  facet_grid(~facet) +
  labs(title="Policy Support by age category", 
       y="Change in Support", 
       x="") +
  scale_color_brewer(palette="Set1") + 
  theme_minimal() +
  theme(axis.text.y = element_text(hjust = 0), 
        axis.title.x=element_text(size=12), 
        legend.position="none")



## Model with education
cjPolicy$xfactor <- cjPolicy$educat
cjPolicy.na <- filter(cjPolicy, !is.na(xfactor))

m3.edu <- amce(support ~ Source.of.Funding:xfactor + Policy.Measure:xfactor + 
                 Costs:xfactor + Exceptions:xfactor + Running.Time:xfactor +
                 Energy.Source.Priority:xfactor + Nuclear.Power.Plants:xfactor,
               data = cjPolicy.na,
               respondent.id = "id",
               respondent.varying = "xfactor",
               baselines = cjPolBase, design = cjPolDesign)

summary(m3.edu)

p.m3 <- plot.amce2(m3.edu, plot.display="interaction")

p.m3.d <- p.m3$plot$data
p.m3.d$x <- seq(1,nrow(p.m3.d), 1)
p.m3.d$x2 <- seq(1,nrow(p.m3.d)/3, 1)
p.m3.d$facet <- recode_factor(p.m3.d$facet, 
                              'Conditional on\nxfactor = low' = "low",
                              'Conditional on\nxfactor = middle' = "middle",
                              'Conditional on\nxfactor = high' = "high"
)

ggplot(filter(p.m3.d, x2 < 17)) +
  geom_hline(yintercept = 0, size = 0.5, colour = "black", 
             linetype = "dotted") +
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper,
                      x = reorder(printvar, -x), color=facet, 
                      shape=facet),
                  size=0.4,
                  position=position_dodge(width=0.95)) +
  ylim(-1,0.5) + coord_flip() + 
  facet_grid(~facet) +
  labs(title="Policy Support by education level", 
       y="Change in Support", 
       x="") +
  scale_color_brewer(palette="Set1") + 
  theme_minimal() +
  theme(axis.text.y = element_text(hjust = 0), 
        axis.title.x=element_text(size=12), 
        legend.position="none")



## Model with income
cjPolicy$xfactor <- cjPolicy$incomecat
cjPolicy.na <- filter(cjPolicy, !is.na(xfactor))

m3.inc <- amce(support ~ Source.of.Funding:xfactor + Policy.Measure:xfactor + 
                 Costs:xfactor + Exceptions:xfactor + Running.Time:xfactor +
                 Energy.Source.Priority:xfactor + Nuclear.Power.Plants:xfactor,
               data = cjPolicy.na,
               respondent.id = "id",
               respondent.varying = "xfactor",
               baselines = cjPolBase, design = cjPolDesign)

summary(m3.inc)

p.m3 <- plot.amce2(m3.inc, plot.display="interaction")

p.m3.d <- p.m3$plot$data
p.m3.d$x <- seq(1,nrow(p.m3.d), 1)
p.m3.d$x2 <- seq(1,nrow(p.m3.d)/3, 1)
p.m3.d$facet <- recode_factor(p.m3.d$facet, 
                              'Conditional on\nxfactor = low' = "low",
                              'Conditional on\nxfactor = middle' = "middle",
                              'Conditional on\nxfactor = high' = "high"
)

ggplot(filter(p.m3.d, x2 < 17)) +
  geom_hline(yintercept = 0, size = 0.5, colour = "black", 
             linetype = "dotted") +
  geom_pointrange(aes(y = pe, ymin = lower, ymax = upper,
                      x = reorder(printvar, -x), color=facet, 
                      shape=facet),
                  size=0.4,
                  position=position_dodge(width=0.95)) +
  ylim(-1,0.5) + coord_flip() + 
  facet_grid(~facet) +
  labs(title="Policy Support by income level", 
       y="Change in Support", 
       x="") +
  scale_color_brewer(palette="Set1") + 
  theme_minimal() +
  theme(axis.text.y = element_text(hjust = 0), 
        axis.title.x=element_text(size=12), 
        legend.position="none")
rm(p.m3, p.m3.d)