#author: Vanessa Schwaiger / vanessa.schwaiger@sowi.uni-stuttgart.de
# Install Packages --------------------------------------------------------

#install.packages("pacman")
library(pacman)
pacman::p_load(ggplot2, tidyr, qdapTools, ggalluvial, ggpubr,tidycomm,
               nnet, lsr, sjmisc, magrittr, MASS, factoextra, FactoMineR, 
               ggstance, jtools, stargazer, stargazer, survey, readxl, broom,
               sjPlot, texreg, tidyverse, magrittr, haven, sjmisc, ggeffects,
               texreg, car, psych, knitr, labelled, margins, lme4,sjstats,
               DescTools, corrplot, cregg,openxlsx, janitor, openxlsx, dplyr,cregg,
               GDAtools,expss)





# Read Data Set -----------------------------------------------------------
#con1 <- read.csv(".../dataset_sticktothestatusquo.csv")

# Rename Feature Labels for Analysis --------------------------------------
con1 <- con1 %>% mutate_at(c('feature1', 'feature2', 'feature3','feature4', 'feature5', 'feature6' ), as.factor)
is.numeric(con1$ID)
table(con1$ID)
table(con1$feature1)
is.factor(con1$feature1)
levels(con1$feature1)<-c("Citizen Forum", "Assertive Leader", "Parliament")
levels(con1$feature1)
table(con1$feature1)

table(con1$feature2)
is.factor(con1$feature2)
levels(con1$feature2)<-c("Consultation with Experts", "Consideration of Public Opinion", "On their own")
levels(con1$feature2)
table(con1$feature2)


table(con1$feature3)
is.factor(con1$feature3)
levels(con1$feature3)<-c("Considered and Slow", "Efficient and Fast")
levels(con1$feature3)
table(con1$feature3)

table(con1$feature4)
is.factor(con1$feature4)
levels(con1$feature4)<-c("Binding Decision", "Inclusion Parliament", "Final Referendum")
levels(con1$feature4)
table(con1$feature4)


table(con1$feature5)
is.factor(con1$feature5)
levels(con1$feature5)<-c("In Favor of Measurement", "Against Measurement")
levels(con1$feature5)
table(con1$feature5)


table(con1$feature6)
is.factor(con1$feature6)
levels(con1$feature6)<-c("Preference Mismatch", "Preference Match")
levels(con1$feature6)
table(con1$feature6)

summary(con1$ID)

#feature names

feature_names <- c(
  'feature1' = "Main institution/actor",
  'feature2' = "Additional Involvement",
  'feature3' = "Process",
  'feature4' = "Decision",
  'feature5' = "Result",
  'feature6' = "Outcome Favorability"
)





# AMCE Baseline Model Choice-----------------------------------------------------



# Relevel Features --------------------------------------------------------

con1$feature1 <- factor(con1$feature1, levels = c("Parliament", "Assertive Leader", "Citizen Forum"))
con1$feature2 <- factor(con1$feature2, levels = c("On their own", "Consultation with Experts", "Consideration of Public Opinion"))




#### baseline models ####




amces <- chosen ~ feature1 * feature4 + feature2 + feature3   +feature6
plot(cj_freqs(con1, amces, id = ~ID))

amce_choice <- cj(con1, amces, id = ~ID,  estimate="amce", 
                  feature_labels = feature_names, weights = con1$scale_base_weights_f,
                  feature_order = c("feature1", "feature2", "feature3", "feature4", "feature6"))


amce_choice

#plot model
amce_plot<-plot(amce_choice, xlim = c(-0.2,0.2))+  
  theme(axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        axis.ticks = element_blank(), 
        legend.position = "none",
        panel.spacing.y = unit(1.5, "lines"),
        panel.border = element_rect(linetype = "blank"), 
        strip.background= element_rect(size = 2, colour="grey89", fill = "white"),
        strip.text.y.left = element_text(angle = 0),
        strip.placement = "outside") + 
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_brewer(palette="Dark2")

amce_plot


#save plot
#ggsave("amce_basemodel.png", width = 20, height=17, units = "cm", dpi=800, scale=1, bg="white") 
#ggsave("amce_basemodel_highres.png", width = 20, height=17, units = "cm", dpi=1200, scale=1, bg="white") 






# Subgroup Analysis Choice  ------------------------------------------------------
# make subsets (recommended by Leeper) ------------------------------------

#make subset without "inclusion parliament"
table(con1$feature4)
con1_minu_inclusion_parl<-con1 %>% filter(feature4 == "Binding Decision" |
                                            feature4 == "Final Referendum") %>% droplevels()

table(con1_minu_inclusion_parl$feature4)



#make subset without "parliament"
table(con1$feature1)
con1_minu_parl<-con1 %>% filter(feature1 == "Citizen Forum" | feature1 == "Assertive Leader") %>%
  droplevels()
table(con1_minu_parl$feature1)




#compare level distribution in base dataset, without inclusion parliament and without parliament
f1 <- chosen ~ feature1 * feature4 + feature2 + feature3   +feature6
f2 <- chosen ~ feature1  + feature2 + feature3 +feature4+  +feature6

base_distr<-plot(cj_freqs(con1, f2, id = ~ID))
base_distr + theme(legend.position="none")
#base data --> inclusion parliament underrepr. because of restriction
#ggsave("basemodel_distribution.png", width = 12, height=10, units = "cm", dpi=800, scale=1.2, bg="white") 



plot(cj_freqs(con1_minu_inclusion_parl, f2, id = ~ID))
#without inclusion parliament --> parliament overrepresented 

subsetB<-plot(cj_freqs(con1_minu_parl, f2, id = ~ID))
subsetB + theme(legend.position="none")

#without parliament --> equally distribution across features
#ggsave("subsetB_distribution.png", width = 12, height=10, units = "cm", dpi=800, scale=1.2, bg="white") 


#without inclusion parliament, parliament is oversampled. Select randomly 2/3 of the cases so it is distributed equally

# Define the fraction of cases to select
fraction <- 2/3


# Set a seed for reproducibility
set.seed(456)

# Randomly select fraction of cases from each category

con1_minu_inclusion_parl_rand <- do.call(rbind, lapply(unique(con1_minu_inclusion_parl$feature1), 
                                                       function(x) {
                                                         if (x == "Parliament") {
                                                           n <- sum(con1_minu_inclusion_parl$feature1 == "Parliament")
                                                           n_select <- ceiling(fraction * n)
                                                           con1_minu_inclusion_parl[sample(which(con1_minu_inclusion_parl$feature1 == "Parliament"), 
                                                                                           n_select), ]
                                                         } else {
                                                           con1_minu_inclusion_parl[con1_minu_inclusion_parl$feature1 == x, ]
                                                         }
                                                       }))

#save dataset
write.csv(con1_minu_inclusion_parl_rand, "random_selection_parliament.csv", row.names = FALSE)

#check if equally distributed
subsetA<-plot(cj_freqs(con1_minu_inclusion_parl_rand, f2, id = ~ID))
subsetA + theme(legend.position="none")
#now equally distrbuted
#ggsave("subsetA_distribution.png", width = 12, height=10, units = "cm", dpi=800, scale=1.2, bg="white") 



#data set for subgroup analysis
#con1 = all data
#con1_minu_inclusion_parl_rand = without inclusion parliament (randomized 2/3 of parliament cases) 
#con1_minu_parl


# Differences in Conditional AMCEs Choice ----------------------------------------




# Satisfaction (AMCEs) ------------------------------------------------------------


#compute binary variable satisfaction for each dataset / subset

#base
table(con1$demo_satis_dicho)
con1$demo_satis_dicho_nu<-as.numeric(con1$demo_satis_dicho)
table(con1$demo_satis_dicho_nu)
con1$demo_satis_dicho_rec<-(con1$demo_satis_dicho_nu-2)*(-1)
table(con1$demo_satis_dicho_rec)

con1$demo_satis_dicho_rec<-as.factor(con1$demo_satis_dicho_rec)
table(con1$demo_satis_dicho_rec)
is.factor(con1$demo_satis_dicho_rec)
levels(con1$demo_satis_dicho_rec)<-c("High Satis.", "Low Satis.")
table(con1$demo_satis_dicho_rec)

#without inclusion parliament

table(con1_minu_inclusion_parl_rand$demo_satis_dicho)
con1_minu_inclusion_parl_rand$demo_satis_dicho_nu<-as.numeric(con1_minu_inclusion_parl_rand$demo_satis_dicho)
table(con1_minu_inclusion_parl_rand$demo_satis_dicho_nu)
con1_minu_inclusion_parl_rand$demo_satis_dicho_rec<-(con1_minu_inclusion_parl_rand$demo_satis_dicho_nu-2)*(-1)
table(con1_minu_inclusion_parl_rand$demo_satis_dicho_rec)

con1_minu_inclusion_parl_rand$demo_satis_dicho_rec<-as.factor(con1_minu_inclusion_parl_rand$demo_satis_dicho_rec)
table(con1_minu_inclusion_parl_rand$demo_satis_dicho_rec)
is.factor(con1_minu_inclusion_parl_rand$demo_satis_dicho_rec)
levels(con1_minu_inclusion_parl_rand$demo_satis_dicho_rec)<-c("High Satis.", "Low Satis.")
table(con1_minu_inclusion_parl_rand$demo_satis_dicho_rec)



#without parliament
table(con1_minu_parl$demo_satis_dicho)
con1_minu_parl$demo_satis_dicho_nu<-as.numeric(con1_minu_parl$demo_satis_dicho)
table(con1_minu_parl$demo_satis_dicho_nu)
con1_minu_parl$demo_satis_dicho_rec<-(con1_minu_parl$demo_satis_dicho_nu-2)*(-1)
table(con1_minu_parl$demo_satis_dicho_rec)

con1_minu_parl$demo_satis_dicho_rec<-as.factor(con1_minu_parl$demo_satis_dicho_rec)
table(con1_minu_parl$demo_satis_dicho_rec)
is.factor(con1_minu_parl$demo_satis_dicho_rec)
levels(con1_minu_parl$demo_satis_dicho_rec)<-c("High Satis.", "Low Satis.")
table(con1_minu_parl$demo_satis_dicho_rec)


    # all - no subsets by constraint ------------------------------------------


#amces satisfaction
amces_satis<- cj((con1), chosen ~ feature1+feature2+ feature3+feature4 +feature6, id = ~ID, estimate = "amce", 
                 weights = ~scale_base_weights_f, 
                 by = ~demo_satis_dicho_rec, feature_labels = feature_names,
                 feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))
diff_amces_satis <- cj((con1), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff",
                       weights = ~scale_base_weights_f,
                       by = ~demo_satis_dicho_rec, feature_labels=feature_names,
                       feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1")) 


#plot conditional amces satisfaction
plot(rbind(amces_satis, diff_amces_satis), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#save conditional amces satisfaction
#ggsave("subgroup_satisfaction_amces.png", width = 22, height=19, units = "cm", dpi=800, scale=1.2, bg="white") 



    # subsets by constraint ---------------------------------------------------
#test if mediansplit still correct within subsets
#full
weighted.mean(con1$demo_satis, con1$scale_base_weights_f, na.rm=TRUE)
weighted_median(con1$demo_satis, weights=con1$scale_base_weights_f)

#subset A
weighted.mean(con1_minu_inclusion_parl_rand$demo_satis, con1_minu_inclusion_parl_rand$scale_base_weights_f, na.rm=TRUE)
weighted_median(con1_minu_inclusion_parl_rand$demo_satis, weights=con1_minu_inclusion_parl_rand$scale_base_weights_f)
#same implications for median split as full set

#subset B
weighted.mean(con1_minu_parl$demo_satis, con1_minu_parl$scale_base_weights_f, na.rm=TRUE)
weighted_median(con1_minu_parl$demo_satis, weights=con1_minu_parl$scale_base_weights_f)
#same implications for median split as full set



#without inclusion parliament

amces_satis_minu_inclusion_parl_rand<- cj((con1_minu_inclusion_parl_rand), chosen ~ feature1+feature2+ feature3+feature4 +feature6, id = ~ID, estimate = "amce", 
                                          weights = ~scale_base_weights_f, 
                                          by = ~demo_satis_dicho_rec, feature_labels = feature_names,
                                          feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))
diff_amces_satis_minu_inclusion_parl_rand <- cj((con1_minu_inclusion_parl_rand), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff",
                                                weights = ~scale_base_weights_f,
                                                by = ~demo_satis_dicho_rec, feature_labels=feature_names,
                                                feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1")) 


#plot conditional amces satisfaction
plot(rbind(amces_satis_minu_inclusion_parl_rand, diff_amces_satis_minu_inclusion_parl_rand), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#save conditional amces satisfaction without inclusion
#ggsave("subgroup_satisfaction_amces_without_inclusion_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1.2, bg="white") 


#without parliament

amces_satis_minu_parl<- cj((con1_minu_parl), chosen ~ feature1+feature2+ feature3+feature4 +feature6, id = ~ID, estimate = "amce", 
                           weights = ~scale_base_weights_f, 
                           by = ~demo_satis_dicho_rec, feature_labels = feature_names,
                           feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))
diff_amces_satis_minu_parl <- cj((con1_minu_parl), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff",
                                 weights = ~scale_base_weights_f,
                                 by = ~demo_satis_dicho_rec, feature_labels=feature_names,
                                 feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1")) 


#plot conditional amces satisfaction
plot(rbind(amces_satis_minu_parl, diff_amces_satis_minu_parl), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#save conditional amces satisfaction without parliament
#ggsave("subgroup_satisfaction_amces_without_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1.2, bg="white") 
















# Political Interest(AMCEs) ------------------------------------------------------

#compute binary variable pol interest


#base
con1$pol_int_dicho<-as.factor(con1$pol_int_dicho)
table(con1$pol_int_dicho)
is.factor(con1$pol_int_dicho)
levels(con1$pol_int_dicho)<-c("Low pol. Int.", "High pol. Int.")
table(con1$pol_int_dicho)
length(con1$pol_int_dicho)

#without inclusion parliament
con1_minu_inclusion_parl_rand$pol_int_dicho<-as.factor(con1_minu_inclusion_parl_rand$pol_int_dicho)
table(con1_minu_inclusion_parl_rand$pol_int_dicho)
is.factor(con1_minu_inclusion_parl_rand$pol_int_dicho)
levels(con1_minu_inclusion_parl_rand$pol_int_dicho)<-c("Low pol. Int.", "High pol. Int.")
table(con1_minu_inclusion_parl_rand$pol_int_dicho)
length(con1_minu_inclusion_parl_rand$pol_int_dicho)


#without parliament
con1_minu_parl$pol_int_dicho<-as.factor(con1_minu_parl$pol_int_dicho)
table(con1_minu_parl$pol_int_dicho)
is.factor(con1_minu_parl$pol_int_dicho)
levels(con1_minu_parl$pol_int_dicho)<-c("Low pol. Int.", "High pol. Int.")
table(con1_minu_parl$pol_int_dicho)
length(con1_minu_parl$pol_int_dicho)




    # all - no subsets by constraint ------------------------------------------

amces_polint<- cj((con1), chosen ~ feature1+feature2+ feature3+feature4 +feature6, id = ~ID, estimate = "amce",
                  weights = ~scale_base_weights_f,
                  by = ~pol_int_dicho, feature_labels = feature_names,
                  feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))
diff_amces_polint <- cj((con1), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff",
                        weights = ~scale_base_weights_f,
                        by = ~pol_int_dicho, feature_labels=feature_names,
                        feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces pol interest
plot(rbind(amces_polint, diff_amces_polint), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#save conditional amces pol interest
#ggsave("subgroup_pol_int_amces.png", width = 22, height=19, units = "cm", dpi=800, scale=1.2, bg="white") 


    # subsets by constraint ---------------------------------------------------
#test if mediansplit still correct within subsets
#full
weighted.mean(con1$pol_int, con1$scale_base_weights_f, na.rm=TRUE)
weighted_median(con1$pol_int, weights=con1$scale_base_weights_f)

#subset A
weighted.mean(con1_minu_inclusion_parl_rand$pol_int, con1_minu_inclusion_parl_rand$scale_base_weights_f, na.rm=TRUE)
weighted_median(con1_minu_inclusion_parl_rand$pol_int, weights=con1_minu_inclusion_parl_rand$scale_base_weights_f)
#same implications for median split as full set

#subset B
weighted.mean(con1_minu_parl$pol_int, con1_minu_parl$scale_base_weights_f, na.rm=TRUE)
weighted_median(con1_minu_parl$pol_int, weights=con1_minu_parl$scale_base_weights_f)
#same implications for median split as full set


#AMCEs subset without inclusion parliament
amces_polint_minu_inclusion_parl_rand<- cj((con1_minu_inclusion_parl_rand), chosen ~ feature1+feature2+ feature3+feature4 +feature6, id = ~ID, estimate = "amce",
                                           weights = ~scale_base_weights_f,
                                           by = ~pol_int_dicho, feature_labels = feature_names,
                                           feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))
diff_amces_polint_minu_inclusion_parl_rand <- cj((con1_minu_inclusion_parl_rand), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff",
                                                 weights = ~scale_base_weights_f,
                                                 by = ~pol_int_dicho, feature_labels=feature_names,
                                                 feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces pol interest
plot(rbind(amces_polint_minu_inclusion_parl_rand, diff_amces_polint_minu_inclusion_parl_rand), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#save conditional amces pol interest
#ggsave("subgroup_pol_int_amces_without_inclusion_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1.2, bg="white") 



#AMCEs subset without parliament
amces_polint_minu_parl<- cj((con1_minu_parl), chosen ~ feature1+feature2+ feature3+feature4 +feature6, id = ~ID, estimate = "amce",
                            weights = ~scale_base_weights_f,
                            by = ~pol_int_dicho, feature_labels = feature_names,
                            feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))
diff_amces_polint_minu_parl <- cj((con1_minu_parl), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff",
                                  weights = ~scale_base_weights_f,
                                  by = ~pol_int_dicho, feature_labels=feature_names,
                                  feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces pol interest
plot(rbind(amces_polint_minu_parl, diff_amces_polint_minu_parl), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#save conditional amces pol interest
#ggsave("subgroup_pol_int_amces_without_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1.2, bg="white") 











# Education (AMCEs) -------------------------------------------------------

#binary variable education


#base
table(con1$gym)

con1$gym<-as.factor(con1$gym)
levels(con1$gym)<-c("Educ. low", "Educ. high")
table(con1$gym)


#without inclusion parliament
table(con1_minu_inclusion_parl_rand$gym)

con1_minu_inclusion_parl_rand$gym<-as.factor(con1_minu_inclusion_parl_rand$gym)
levels(con1_minu_inclusion_parl_rand$gym)<-c("Educ. low", "Educ. high")
table(con1_minu_inclusion_parl_rand$gym)


#without parliament
table(con1_minu_parl$gym)

con1_minu_parl$gym<-as.factor(con1_minu_parl$gym)
levels(con1_minu_parl$gym)<-c("Educ. low", "Educ. high")
table(con1_minu_parl$gym)
is.factor(con1_minu_parl$gym)


    # all - no subsets by constraint ------------------------------------------

amces_gym<- cj((con1), chosen ~ feature1+feature2+ feature3+feature4 +feature6, id = ~ID, estimate = "amce", 
               weights = ~scale_base_weights_f,
               by = ~gym, feature_labels = feature_names,
               feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))
diff_amces_gym <- cj((con1), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff", 
                     weights = ~scale_base_weights_f,
                     by = ~gym, feature_labels=feature_names,
                     feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces education
plot(rbind(amces_gym, diff_amces_gym), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))


#save conditional amces education
#ggsave("subgroup_education_amces.png", width = 22, height=19, units = "cm", dpi=800, scale=1.2, bg="white") 



    # subsets by constraint ---------------------------------------------------


#amces education subset without inclusion parliament
amces_gym_minu_inclusion_parl<- cj((con1_minu_inclusion_parl_rand), chosen ~ feature1+feature2+ feature3+feature4 +feature6, id = ~ID, estimate = "amce", 
                                   weights = ~scale_base_weights_f,
                                   by = ~gym, feature_labels = feature_names,
                                   feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))
diff_amces_gym_minu_inclusion_parl <- cj((con1_minu_inclusion_parl_rand), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff", 
                                         weights = ~scale_base_weights_f,
                                         by = ~gym, feature_labels=feature_names,
                                         feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces education
plot(rbind(amces_gym_minu_inclusion_parl, diff_amces_gym_minu_inclusion_parl), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))


#save conditional amces education subset without inclusion parliament
#ggsave("subgroup_education_amces_without_inclusion_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1.2, bg="white") 



#amces education without parliament
amces_gym_minu_parl<- cj((con1_minu_parl), chosen ~ feature1+feature2+ feature3+feature4 +feature6, id = ~ID, estimate = "amce", 
                         weights = ~scale_base_weights_f,
                         by = ~gym, feature_labels = feature_names,
                         feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))
diff_amces_gym_minu_parl <- cj((con1_minu_parl), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff", 
                               weights = ~scale_base_weights_f,
                               by = ~gym, feature_labels=feature_names,
                               feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces education
plot(rbind(amces_gym_minu_parl, diff_amces_gym_minu_parl), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))


#save conditional amces education subset without inclusion parliament
#ggsave("subgroup_education_amces_without_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1.2, bg="white") 










# overview subgroups differences (AMCEs) -----------------------------------------------

    # all - no subsets by constraint ------------------------------------------


#compute AMCEs / conditional AMCEs
amces_subgroup_choice <- chosen ~ feature1 + feature2 + feature3  +feature4 +feature6


#base
diff_amce_gym_choice <- cj(con1, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~gym)
diff_amce_interest_choice <- cj(con1, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~pol_int_dicho)
diff_amce_dissatisfied_choice <- cj(con1, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~demo_satis_dicho_rec)

#prepare to plot
diff_estimates_choice <-bind_rows(
  diff_amce_dissatisfied_choice %>% mutate(model = "Dissatisfied"),
  diff_amce_interest_choice %>% mutate(model = "Politically interested"),
  diff_amce_gym_choice %>% mutate(model = "Education high"))

diff_choice <- diff_estimates_choice %>% dplyr::select(BY:level, estimate:model)
diff_choice$model <- factor(diff_choice$model, levels = c("Dissatisfied",  "Politically interested","Education high"))


#plot overview
plot(diff_choice, feature_headers = FALSE, xlab = "\nEstimated Difference (AMCE) for different groups of adolescents") + 
  ggplot2::facet_wrap(~model, ncol = 3) +
  theme(
    panel.border = element_rect(linetype = "dotted", colour = "black"),
    panel.background = element_blank(), 
    axis.text.y = element_text(colour = "black", size = 12),
    axis.text.x = element_text(colour = "black", size = 12),
    axis.title = element_text(colour = "black", size = 12),
    panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
    strip.background = element_blank(),
    strip.text.x = element_text(colour = "black", size = 15),
    legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))


#ggsave("subgroup_overview_differences_amces.png", width = 22, height=19, units = "cm", dpi=800, scale=1, bg="white") 
#ggsave("subgroup_overview_differences_amces_highres.png", width = 22, height=19, units = "cm", dpi=1200, scale=1, bg="white") 


    # subsets by constraint ---------------------------------------------------


#wihtout inclusion parliament
diff_amce_gym_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~gym)
diff_amce_interest_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~pol_int_dicho)
diff_amce_dissatisfied_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~demo_satis_dicho_rec)

#prepare to plot
diff_estimates_choice_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_dissatisfied_choice_minu_inclusion_parl_rand %>% mutate(model = "Dissatisfied"),
  diff_amce_interest_choice_minu_inclusion_parl_rand %>% mutate(model = "Politically interested"),
  diff_amce_gym_choice_minu_inclusion_parl_rand %>% mutate(model = "Education high"))

diff_choice_minu_inclusion_parl_rand <- diff_estimates_choice_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_choice_minu_inclusion_parl_rand$model <- factor(diff_choice_minu_inclusion_parl_rand$model, levels = c("Dissatisfied",  "Politically interested","Education high"))


#plot overview
plot(diff_choice_minu_inclusion_parl_rand, feature_headers = FALSE, xlab = "\nEstimated Difference (AMCE) for different groups of adolescents") + 
  ggplot2::facet_wrap(~model, ncol = 3) +
  theme(
    panel.border = element_rect(linetype = "dotted", colour = "black"),
    panel.background = element_blank(), 
    axis.text.y = element_text(colour = "black", size = 12),
    axis.text.x = element_text(colour = "black", size = 12),
    axis.title = element_text(colour = "black", size = 12),
    panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
    strip.background = element_blank(),
    strip.text.x = element_text(colour = "black", size = 15),
    legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#ggsave("subgroup_overview_differences_amces_without_inclusion_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1, bg="white") 


#without parliament
diff_amce_gym_choice_minu_parl <- cj(con1_minu_parl, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~gym)
diff_amce_interest_choice_minu_parl <- cj(con1_minu_parl, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~pol_int_dicho)
diff_amce_dissatisfied_choice_minu_parl <- cj(con1_minu_parl, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~demo_satis_dicho_rec)

#prepare to plot
diff_estimates_choice_minu_parl <-bind_rows(
  diff_amce_dissatisfied_choice_minu_parl %>% mutate(model = "Dissatisfied"),
  diff_amce_interest_choice_minu_parl %>% mutate(model = "Politically interested"),
  diff_amce_gym_choice_minu_parl %>% mutate(model = "Education high"))

diff_choice_minu_parl <- diff_estimates_choice_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_choice_minu_parl$model <- factor(diff_choice_minu_parl$model, levels = c("Dissatisfied",  "Politically interested","Education high"))


#plot overview
plot(diff_choice_minu_parl, feature_headers = FALSE, xlab = "\nEstimated Difference (AMCE) for different groups of adolescents") + 
  ggplot2::facet_wrap(~model, ncol = 3) +
  theme(
    panel.border = element_rect(linetype = "dotted", colour = "black"),
    panel.background = element_blank(), 
    axis.text.y = element_text(colour = "black", size = 12),
    axis.text.x = element_text(colour = "black", size = 12),
    axis.title = element_text(colour = "black", size = 12),
    panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
    strip.background = element_blank(),
    strip.text.x = element_text(colour = "black", size = 15),
    legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))


#ggsave("subgroup_overview_differences_amces_without_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1, bg="white") 












# ########## APPENDIX ########## ----------------------------------------------------------------





# Appendix 4 Variables / Mean ---------------------------------------------


# Rating / Decision Acceptance --------------------------------------------
summary(con1$rating)
weighted.mean(con1$rating, con1$scale_base_weights_f, na.rm=TRUE)
weighted_sd(con1$rating, w=con1$scale_base_weights_f)
sd(con1$rating, na.rm = TRUE)
weighted_median(con1$rating, weights=con1$scale_base_weights_f)


summary(con1$info_text_dauer_dicho)
summary(con1$timing_conjoint_num_test.x)

# time --------------------------------------------------------------------
#conjoint time
weighted.mean(con1$timing_conjoint_num_test.x, con1$scale_base_weights_f, na.rm=TRUE)
weighted_sd(con1$timing_conjoint_num_test.x, w=con1$scale_base_weights_f)
sd(con1$timing_conjoint_num_test.x, na.rm = TRUE)
weighted_median(con1$timing_conjoint_num_test.x, weights=con1$scale_base_weights_f)

#info text time
weighted.mean(con1$info_text_dauer, con1$scale_base_weights_f, na.rm=TRUE)
weighted_sd(con1$info_text_dauer, w=con1$scale_base_weights_f)
sd(con1$info_text_dauer, na.rm = TRUE)
weighted_median(con1$info_text_dauer, weights=con1$scale_base_weights_f)


# Issue Salience ----------------------------------------------------------
#clim_ch_1
con1$clim_ch_1<-as.numeric(con1$clim_ch_1)
weighted.mean(con1$clim_ch_1, con1$scale_base_weights_f, na.rm=TRUE)
summary(con1$clim_ch_1)
weighted_median(con1$clim_ch_1, weights=con1$scale_base_weights_f)
weighted_sd(con1$clim_ch_1, w=con1$scale_base_weights_f)
sd(con1$clim_ch_1, na.rm=TRUE)



# Worrys Climate Change ---------------------------------------------------
#worry_topics_1
con1$worry_topics_1<-as.numeric(con1$worry_topics_1)
con1$clim_ch_worry
weighted.mean(con1$worry_topics_1, con1$scale_base_weights_f, na.rm=TRUE)
summary(con1$worry_topics_1)
weighted_median(con1$worry_topics_1, weights=con1$scale_base_weights_f)
weighted_sd(con1$worry_topics_1, w=con1$scale_base_weights_f)
sd(con1$worry_topics_1, na.rm=TRUE)

con1$worry_topics_1<-as.numeric(con1$worry_topics_1)



# preference --------------------------------------------------------------
#measure_rating
weighted.mean(con1$measure_rating, con1$scale_base_weights_f, na.rm=TRUE)
summary(con1$measure_rating)
weighted_median(con1$measure_rating, weights=con1$scale_base_weights_f)
weighted_sd(con1$measure_rating, w=con1$scale_base_weights_f)
sd(con1$measure_rating, na.rm=TRUE)



#subgroups

#democratic satisfaction
table(con1$demo_satis)
length(con1$demo_satis)
table(con1$demo_satis_dicho_rec) #turned around
weighted.mean(con1$demo_satis, con1$scale_base_weights_f, na.rm=TRUE)
summary(con1$demo_satis)
weighted_median(con1$demo_satis, weights=con1$scale_base_weights_f)
weighted_sd(con1$measure_rating, w=con1$scale_base_weights_f)
sd(con1$demo_satis, na.rm=TRUE)

#pol interest
weighted.mean(con1$pol_int, con1$scale_base_weights_f, na.rm=TRUE)
summary(con1$pol_int)
weighted_median(con1$pol_int, weights=con1$scale_base_weights_f)
weighted_sd(con1$measure_rating, w=con1$scale_base_weights_f)
sd(con1$pol_int, na.rm=TRUE)

table(con1$info_text_dauer_dicho)
table(con1$timing_conjoint_num.x)


# t-test rating dissatisfaction pol. interest and gym ---------------------

#dissatisfaction
# Perform weighted t-test
weighted_t_test <- t.test(rating ~ demo_satis_dicho, 
                          data = con1, 
                          weights = scaled_base_weights_f)

# Print the results
print(weighted_t_test)


#pol interest
# Perform weighted t-test
weighted_t_test <- t.test(rating ~ pol_int_dicho, 
                          data = con1, 
                          weights = scaled_base_weights_f)

# Print the results
print(weighted_t_test)


#gym
# Perform weighted t-test
weighted_t_test <- t.test(rating ~ gym, 
                          data = con1, 
                          weights = scaled_base_weights_f)

# Print the results
print(weighted_t_test)





# Appendix 6 Robustness checks --------------------------------------------------------------


# learning and consistency ------------------------------------------------


#how many full data
table(con1$task)

con_task_1<- con1 %>% filter(task== 1)
con_task_2<- con1 %>% filter(task== 2)
con_task_3<- con1 %>% filter(task== 3)
con_task_4<- con1 %>% filter(task== 4)
con_task_5<- con1 %>% filter(task== 5)
con_task_6<- con1 %>% filter(task== 6)

con_task_1_2<- con1 %>% filter(task== 1 | task==2)


#how many subset 1
con_task_1_minu_inclusion_parl_rand<- con1_minu_inclusion_parl_rand %>% filter(task== 1)
con_task_2_minu_inclusion_parl_rand<- con1_minu_inclusion_parl_rand %>% filter(task== 2)
con_task_3_minu_inclusion_parl_rand<- con1_minu_inclusion_parl_rand %>% filter(task== 3)
con_task_4_minu_inclusion_parl_rand<- con1_minu_inclusion_parl_rand %>% filter(task== 4)
con_task_5_minu_inclusion_parl_rand<- con1_minu_inclusion_parl_rand %>% filter(task== 5)
con_task_6_minu_inclusion_parl_rand<- con1_minu_inclusion_parl_rand %>% filter(task== 6)


#how many subset 2
con_task_1_minu_parl<- con1_minu_inclusion_parl %>% filter(task== 1)
con_task_2_minu_parl<- con1_minu_inclusion_parl %>% filter(task== 2)
con_task_3_minu_parl<- con1_minu_inclusion_parl %>% filter(task== 3)
con_task_4_minu_parl<- con1_minu_inclusion_parl %>% filter(task== 4)
con_task_5_minu_parl<- con1_minu_inclusion_parl %>% filter(task== 5)
con_task_6_minu_parl<- con1_minu_inclusion_parl %>% filter(task== 6)


    #Option A: each task individual ------------------------------------------


#### BASE
con1 <- con1 %>% mutate_at(c('task'), as.factor)
levels(con1$task)

#prepare plot
amces_each_task<- cj((con1), chosen ~ feature1+feature2+ feature3+feature4 +feature6, id = ~ID, estimate = "amce", 
                     weights = ~scale_base_weights_f, 
                     by = ~task, feature_labels = feature_names)

#plot 

amces_individual_task_plot<- plot(amces_each_task, group = "task")+theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
                                                                         panel.background = element_blank(),
                                                                         panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
                                                                         strip.background = element_blank(),
                                                                         legend.position = "bottom")
amces_individual_task_plot

#ggsave("A6_individual_tasks.png", width = 22, height=19, units = "cm", dpi=800, scale=1, bg="white") 



    # Option B: add 3 task ----------------------------------------------------------------


#### BASE

#new dataset to test
con1_contest_data<-con1

#contest round 1

con1_contest_rounds_1<- con1%>% filter(task == 1)
table(con1_contest_rounds_1$task)
con1_contest_rounds_1$task_opt <- NA
con1_contest_rounds_1$task_opt[con1_contest_rounds_1$task == "1"]<-1
table(con1_contest_rounds_1$task_opt)


#contest round 1 - 2
con1_contest_rounds_1_2<- con1%>% filter(task == 1 | task ==2)
table(con1_contest_rounds_1_2$task)
con1_contest_rounds_1_2$task_opt <- NA
con1_contest_rounds_1_2$task_opt[con1_contest_rounds_1_2$task == "1"  | con1_contest_rounds_1_2$task == "2"]<-2
table(con1_contest_rounds_1_2$task_opt)



#contest round 1 -3 
con1_contest_rounds_1_3<- con1%>% filter(task == 1 | task ==2 | task==3)
table(con1_contest_rounds_1_3$task)
con1_contest_rounds_1_3$task_opt <- NA
con1_contest_rounds_1_3$task_opt[con1_contest_rounds_1_3$task == "1"  | con1_contest_rounds_1_3$task == "2" | con1_contest_rounds_1_3$task == "3"]<-3
table(con1_contest_rounds_1_3$task_opt)

#contest round 1-4
con1_contest_rounds_1_4<- con1%>% filter(task == 1 | task ==2 | task==3 | task==4)
table(con1_contest_rounds_1_4$task)
con1_contest_rounds_1_4$task_opt <- NA
con1_contest_rounds_1_4$task_opt[con1_contest_rounds_1_4$task == "1"  | con1_contest_rounds_1_4$task == "2" | con1_contest_rounds_1_4$task == "3"|
                                   con1_contest_rounds_1_4$task == "4"]<-4
table(con1_contest_rounds_1_4$task_opt)


#contest round 1-5
con1_contest_rounds_1_5<- con1%>% filter(task == 1 | task ==2 | task==3 | task==4 | task==5)
table(con1_contest_rounds_1_5$task)
con1_contest_rounds_1_5$task_opt <- NA
con1_contest_rounds_1_5$task_opt[con1_contest_rounds_1_5$task == "1"  | con1_contest_rounds_1_5$task == "2" | con1_contest_rounds_1_5$task == "3"|
                                   con1_contest_rounds_1_5$task == "4" | con1_contest_rounds_1_5$task == "5"]<-5
table(con1_contest_rounds_1_5$task_opt)

#contest round 1-6
con1_contest_rounds_1_6<- con1%>% filter(task == 1 | task ==2 | task==3 | task==4 | task==5 |task ==6)
table(con1_contest_rounds_1_6$task)
con1_contest_rounds_1_6$task_opt <- NA
con1_contest_rounds_1_6$task_opt[con1_contest_rounds_1_6$task == "1"  | con1_contest_rounds_1_6$task == "2" | con1_contest_rounds_1_6$task == "3"|
                                   con1_contest_rounds_1_6$task == "4" | con1_contest_rounds_1_6$task == "5" | con1_contest_rounds_1_6$task == "6"]<-6
table(con1_contest_rounds_1_6$task_opt)


#make dataset
con1_contest_added_task<-bind_rows(con1_contest_data,con1_contest_rounds_1, con1_contest_rounds_1_2, con1_contest_rounds_1_3, con1_contest_rounds_1_4,
                                   con1_contest_rounds_1_5, con1_contest_rounds_1_6)
table(con1_contest_added_task$task_opt)

con1_contest_added_task$task<-con1_contest_added_task$task_opt

#filter
con1_contest_rounds_add_one<- con1_contest_added_task%>% filter(task == 1 |task == 2 |task == 3 |task == 4 |task == 5 | task ==6)


#plot
con1_contest_rounds_add_one <- con1_contest_rounds_add_one %>% mutate_at(c('task'), as.factor)
con1_contest_rounds_add_one$task <- factor(con1_contest_rounds_add_one$task, levels = rev(levels(con1_contest_rounds_add_one$task)))
con1_contest_rounds_add_one$task<-plyr::revalue(con1_contest_rounds_add_one$task, c("1"="1", "2"="1-2", "3"= "1-3", "4" = "1-4", "5" = "1-5", "6" ="1-6"))
table(con1_contest_rounds_add_one$task)


amces_added_task<- cj((con1_contest_rounds_add_one), chosen ~ feature1+feature2+ feature3+feature4 +feature6, id = ~ID, estimate = "amce", 
                      weights = ~scale_base_weights_f, 
                      by = ~task, feature_labels = feature_names)



amces_added_task_plot<- plot(amces_added_task, group = "task")+theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
                                                                     panel.background = element_blank(),
                                                                     panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
                                                                     strip.background = element_blank(),
                                                                     legend.position = "bottom")
amces_added_task_plot

#ggsave("A6_added_task.png", width = 22, height=19, units = "cm", dpi=800, scale=1, bg="white") 


    # Option C: task three steps ----------------------------------------------------------------

#### BASE

#filter
con1_contest_rounds_three<- con1_contest_added_task%>% filter(task == 1 |task == 3 | task ==6)

#prepare plot
con1_contest_rounds_three <- con1_contest_rounds_three %>% mutate_at(c('task'), as.factor)
con1_contest_rounds_three$task <- factor(con1_contest_rounds_three$task, levels = rev(levels(con1_contest_rounds_three$task)))
con1_contest_rounds_three$task<-plyr::revalue(con1_contest_rounds_three$task, c("1"="1", "3"= "1-3", "6" ="1-6"))
table(con1_contest_rounds_three$task)


amces_contest_rounds_three<- cj((con1_contest_rounds_three), chosen ~ feature1+feature2+ feature3+feature4 +feature6, id = ~ID, estimate = "amce", 
                                weights = ~scale_base_weights_f, 
                                by = ~task, feature_labels = feature_names)

#plot
amces_three_steps_task_plot<- plot(amces_contest_rounds_three, group = "task")+theme(panel.border = element_rect(linetype = "dotted", colour = "gray17"),
                                                                                     panel.background = element_blank(),
                                                                                     panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
                                                                                     strip.background = element_blank(),
                                                                                     legend.position = "bottom")
amces_three_steps_task_plot

#save
#ggsave("A6_three_added_task.png", width = 22, height=19, units = "cm", dpi=800, scale=1, bg="white") 



# Winner Loser (AMCEs) ------------------------------------------------------------


#preparing

#base
table(con1$feature6)
con1$feature6_rec<-con1$feature6
is.factor(con1$feature6_rec)
levels(con1$feature6_rec)<-c("Mismatch", "Match")
levels(con1$feature6_rec)
table(con1$feature6_rec)
table(con1$feature6)


#without inclusion parl
table(con1_minu_inclusion_parl_rand$feature6)
con1_minu_inclusion_parl_rand$feature6_rec<-con1_minu_inclusion_parl_rand$feature6
is.factor(con1_minu_inclusion_parl_rand$feature6_rec)
levels(con1_minu_inclusion_parl_rand$feature6_rec)<-c("Mismatch", "Match")
levels(con1_minu_inclusion_parl_rand$feature6_rec)
table(con1_minu_inclusion_parl_rand$feature6_rec)
table(con1_minu_inclusion_parl_rand$feature6)

#without parl
table(con1_minu_parl$feature6)
con1_minu_parl$feature6_rec<-con1_minu_parl$feature6
is.factor(con1_minu_parl$feature6_rec)
levels(con1_minu_parl$feature6_rec)<-c("Mismatch", "Match")
levels(con1_minu_parl$feature6_rec)
table(con1_minu_parl$feature6_rec)
table(con1_minu_parl$feature6)


    # all - no subsets by constraint ------------------------------------------

#compute amces winner loser base
amces_winner_loser<- cj((con1), chosen ~ feature1+feature2+ feature3+feature4, id = ~ID, estimate = "amce", 
                        weights = ~scale_base_weights_f,
                        by = ~feature6_rec, feature_labels = feature_names,
                        feature_order = c("feature4", "feature3", "feature2", "feature1"))
diff_amces_winner_loser<- cj((con1), chosen ~ feature1 +feature2+ feature3 +feature4, id = ~ID, estimate = "amce_diff",
                             weights = ~scale_base_weights_f,
                             by = ~feature6_rec, feature_labels=feature_names,
                             feature_order = c("feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces winner loser
plot(rbind(amces_winner_loser, diff_amces_winner_loser), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#save conditional amces winnerloser
#ggsave("subgroup_winner_loser_amces.png", width = 22, height=19, units = "cm", dpi=800, scale=1.2, bg="white") 


    # subsets by constraint ----------------------------------------------------------



#compute amces winner loser without inclusion parl
amces_winner_loser_minu_inclusion_parl_rand<- cj((con1_minu_inclusion_parl_rand), chosen ~ feature1+feature2+ feature3+feature4, id = ~ID, estimate = "amce", 
                                                 weights = ~scale_base_weights_f,
                                                 by = ~feature6_rec, feature_labels = feature_names,
                                                 feature_order = c("feature4", "feature3", "feature2", "feature1"))
diff_amces_winner_loser_minu_inclusion_parl_rand<- cj((con1_minu_inclusion_parl_rand), chosen ~ feature1 +feature2+ feature3 +feature4, id = ~ID, estimate = "amce_diff",
                                                      weights = ~scale_base_weights_f,
                                                      by = ~feature6_rec, feature_labels=feature_names,
                                                      feature_order = c("feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces winner loser
plot(rbind(amces_winner_loser_minu_inclusion_parl_rand, diff_amces_winner_loser_minu_inclusion_parl_rand), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#save conditional amces winnerloser without inclusion parliament
#ggsave("subgroup_winner_loser_amces_without_inclusion_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1.2, bg="white") 




#compute amces winner loser base without parl
amces_winner_loser_minu_parl<- cj((con1_minu_parl), chosen ~ feature1+feature2+ feature3+feature4, id = ~ID, estimate = "amce", 
                                  weights = ~scale_base_weights_f,
                                  by = ~feature6_rec, feature_labels = feature_names,
                                  feature_order = c("feature4", "feature3", "feature2", "feature1"))
diff_amces_winner_loser_minu_parl<- cj((con1_minu_parl), chosen ~ feature1 +feature2+ feature3 +feature4, id = ~ID, estimate = "amce_diff",
                                       weights = ~scale_base_weights_f,
                                       by = ~feature6_rec, feature_labels=feature_names,
                                       feature_order = c("feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces winner loser
plot(rbind(amces_winner_loser_minu_parl, diff_amces_winner_loser_minu_parl), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#save conditional amces winnerloser without parliament
#ggsave("subgroup_winner_loser_amces_without_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1.2, bg="white") 






# Attention Check (AMCEs) ---------------------------------------------------------

#preparing

#base
table(con1$attention)

con1$attention<-as.factor(con1$attention)
levels(con1$attention)<-c("Low Att.", "High Att.")
table(con1$attention)

#without inclusion parl
table(con1_minu_inclusion_parl_rand$attention)

con1_minu_inclusion_parl_rand$attention<-as.factor(con1_minu_inclusion_parl_rand$attention)
levels(con1_minu_inclusion_parl_rand$attention)<-c("Low Att.", "High Att.")
table(con1_minu_inclusion_parl_rand$attention)

#without parl
table(con1_minu_parl$attention)

con1_minu_parl$attention<-as.factor(con1_minu_parl$attention)
levels(con1_minu_parl$attention)<-c("Low Att.", "High Att.")
table(con1_minu_parl$attention)




    # all - no subsets by constraint ------------------------------------------


#computing amces attention check base
amces_attention<- cj((con1), chosen ~ feature1+feature2+ feature3+feature4+feature6, id = ~ID, estimate = "amce",
                     weights = ~scale_base_weights_f,
                     by = ~attention, feature_labels = feature_names,
                     feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))
diff_amces_attention<- cj((con1), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff",
                          weights = ~scale_base_weights_f,
                          by = ~attention, feature_labels=feature_names,
                          feature_order = c("feature6","feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces attention
plot(rbind(amces_attention, diff_amces_attention), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))



#save conditional amces attention
#ggsave("subgroup_attention_amces.png", width = 22, height=19, units = "cm", dpi=800, scale=1.5, bg="white") 



    # subsets by constraint ---------------------------------------------------

#computing amces attention check without inclusion parl
amces_attention_minu_inclusion_parl_rand<- cj((con1_minu_inclusion_parl_rand), chosen ~ feature1+feature2+ feature3+feature4+feature6, id = ~ID, estimate = "amce",
                                              weights = ~scale_base_weights_f,
                                              by = ~attention, feature_labels = feature_names,
                                              feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))
diff_amces_attention_minu_inclusion_parl_rand<- cj((con1_minu_inclusion_parl_rand), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff",
                                                   weights = ~scale_base_weights_f,
                                                   by = ~attention, feature_labels=feature_names,
                                                   feature_order = c("feature6","feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces attention
plot(rbind(amces_attention_minu_inclusion_parl_rand, diff_amces_attention_minu_inclusion_parl_rand), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))



#save conditional amces attention
#ggsave("subgroup_attention_amces_without_inclusion_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1.5, bg="white") 


#computing amces attention check without parl
amces_attention_minu_parl<- cj((con1_minu_parl), chosen ~ feature1+feature2+ feature3+feature4+feature6, id = ~ID, estimate = "amce",
                               weights = ~scale_base_weights_f,
                               by = ~attention, feature_labels = feature_names,
                               feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))
diff_amces_attention_minu_parl<- cj((con1_minu_parl), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff",
                                    weights = ~scale_base_weights_f,
                                    by = ~attention, feature_labels=feature_names,
                                    feature_order = c("feature6","feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces attention
plot(rbind(amces_attention_minu_parl, diff_amces_attention_minu_parl), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))



#save conditional amces attention
#ggsave("subgroup_attention_amces_without_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1.5, bg="white") 








# Conjoint Time (AMCEs) ------------------------------------------------------

#preparing

#base
table(con1$timing_conjoint_num.x)
con1$timing_conjoint_num<-as.factor(con1$timing_conjoint_num.x)
levels(con1$timing_conjoint_num)<-c("less time Conj.", "more time Conj.")
table(con1$timing_conjoint_num)


#without inclusion parl
table(con1_minu_inclusion_parl_rand$timing_conjoint_num.x)
con1_minu_inclusion_parl_rand$timing_conjoint_num<-as.factor(con1_minu_inclusion_parl_rand$timing_conjoint_num.x)
levels(con1_minu_inclusion_parl_rand$timing_conjoint_num)<-c("less time Conj.", "more time Conj.")
table(con1_minu_inclusion_parl_rand$timing_conjoint_num)

#without parl
table(con1_minu_parl$timing_conjoint_num.x)
con1_minu_parl$timing_conjoint_num<-as.factor(con1_minu_parl$timing_conjoint_num.x)
levels(con1_minu_parl$timing_conjoint_num)<-c("less time Conj.", "more time Conj.")
table(con1_minu_parl$timing_conjoint_num)


    # all - no subsets by constraint ------------------------------------------


#computing conditional amces time conj. base
amces_timing_conjoint_num<- cj((con1), chosen ~ feature1+feature2+ feature3+feature4+feature6, id = ~ID, estimate = "amce", 
                               weights = ~scale_base_weights_f,
                               by = ~timing_conjoint_num, feature_labels = feature_names,
                               feature_order = c("feature6","feature4", "feature3", "feature2", "feature1"))
diff_amces_timing_conjoint_num<- cj((con1), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff",
                                    weights = ~scale_base_weights_f,
                                    by = ~timing_conjoint_num, feature_labels=feature_names,
                                    feature_order = c("feature6","feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces time conj.
plot(rbind(amces_timing_conjoint_num, diff_amces_timing_conjoint_num), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#save conditional amces time conj.
#ggsave("subgroup_time_conj_amces.png", width = 22, height=19, units = "cm", dpi=800, scale=1.5, bg="white") 


    # subsets by constraint ---------------------------------------------------
#computing conditional amces time conj. base without inclusion parl
amces_timing_conjoint_num_minu_inclusion_parl_rand<- cj((con1_minu_inclusion_parl_rand), chosen ~ feature1+feature2+ feature3+feature4+feature6, id = ~ID, estimate = "amce", 
                                                        weights = ~scale_base_weights_f,
                                                        by = ~timing_conjoint_num, feature_labels = feature_names,
                                                        feature_order = c("feature6","feature4", "feature3", "feature2", "feature1"))
diff_amces_timing_conjoint_num_minu_inclusion_parl_rand<- cj((con1_minu_inclusion_parl_rand), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff",
                                                             weights = ~scale_base_weights_f,
                                                             by = ~timing_conjoint_num, feature_labels=feature_names,
                                                             feature_order = c("feature6","feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces time conj.
plot(rbind(amces_timing_conjoint_num_minu_inclusion_parl_rand, diff_amces_timing_conjoint_num_minu_inclusion_parl_rand), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#save conditional amces time conj.
#ggsave("subgroup_time_conj_amces_without_inclusion_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1.5, bg="white") 


#computing conditional amces time conj. base without parl
amces_timing_conjoint_num_minu_parl<- cj((con1_minu_parl), chosen ~ feature1+feature2+ feature3+feature4+feature6, id = ~ID, estimate = "amce", 
                                         weights = ~scale_base_weights_f,
                                         by = ~timing_conjoint_num, feature_labels = feature_names,
                                         feature_order = c("feature6","feature4", "feature3", "feature2", "feature1"))
diff_amces_timing_conjoint_num_minu_parl<- cj((con1_minu_parl), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff",
                                              weights = ~scale_base_weights_f,
                                              by = ~timing_conjoint_num, feature_labels=feature_names,
                                              feature_order = c("feature6","feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces time conj.
plot(rbind(amces_timing_conjoint_num_minu_parl, diff_amces_timing_conjoint_num_minu_parl), feature_headers = FALSE) + ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#save conditional amces time conj.
#ggsave("subgroup_time_conj_amces_without_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1.5, bg="white") 





# Information Time (AMCEs) --------------------------------------------------------

#preparing


#base
table(con1$info_text_dauer_dicho)
con1$info_text_dauer_dicho<-as.factor(con1$info_text_dauer_dicho)
levels(con1$info_text_dauer_dicho)<-c("less time Info", "more time Info")
table(con1$info_text_dauer_dicho)

#without inclusion parl
table(con1_minu_inclusion_parl_rand$info_text_dauer_dicho)
con1_minu_inclusion_parl_rand$info_text_dauer_dicho<-as.factor(con1_minu_inclusion_parl_rand$info_text_dauer_dicho)
levels(con1_minu_inclusion_parl_rand$info_text_dauer_dicho)<-c("less time Info", "more time Info")
table(con1_minu_inclusion_parl_rand$info_text_dauer_dicho)

#without parl
table(con1_minu_parl$info_text_dauer_dicho)
con1_minu_parl$info_text_dauer_dicho<-as.factor(con1_minu_parl$info_text_dauer_dicho)
levels(con1_minu_parl$info_text_dauer_dicho)<-c("less time Info", "more time Info")
table(con1_minu_parl$info_text_dauer_dicho)


    # all - no subsets by constraint ------------------------------------------


#computing conditional amces time information base
amces_info_text_dauer_dicho<- cj((con1), chosen ~ feature1+feature2+ feature3+feature4+feature6, id = ~ID, estimate = "amce", 
                                 weights = ~scale_base_weights_f,
                                 by = ~info_text_dauer_dicho, feature_labels = feature_names,
                                 feature_order = c("feature6","feature4", "feature3", "feature2", "feature1"))
diff_amces_info_text_dauer_dicho<- cj((con1), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff",
                                      weights = ~scale_base_weights_f,
                                      by = ~info_text_dauer_dicho, feature_labels=feature_names,
                                      feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces time information 
plot(rbind(amces_info_text_dauer_dicho, diff_amces_info_text_dauer_dicho), feature_headers = FALSE) +ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#save conditional amces time information
#ggsave("subgroup_time_info_amces.png", width = 22, height=19, units = "cm", dpi=800, scale=1.5, bg="white") 



    # subgroup by constraints -------------------------------------------------


#computing conditional amces time information base without inclusion parl
amces_info_text_dauer_dicho_minu_inclusion_parl_rand<- cj((con1_minu_inclusion_parl_rand), chosen ~ feature1+feature2+ feature3+feature4+feature6, id = ~ID, estimate = "amce", 
                                                          weights = ~scale_base_weights_f,
                                                          by = ~info_text_dauer_dicho, feature_labels = feature_names,
                                                          feature_order = c("feature6","feature4", "feature3", "feature2", "feature1"))
diff_amces_info_text_dauer_dicho_minu_inclusion_parl_rand<- cj((con1_minu_inclusion_parl_rand), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff",
                                                               weights = ~scale_base_weights_f,
                                                               by = ~info_text_dauer_dicho, feature_labels=feature_names,
                                                               feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces time information 
plot(rbind(amces_info_text_dauer_dicho_minu_inclusion_parl_rand, diff_amces_info_text_dauer_dicho_minu_inclusion_parl_rand), feature_headers = FALSE) +ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#save conditional amces time information
#ggsave("subgroup_time_info_amces_without_inclusion_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1.5, bg="white") 


#computing conditional amces time information base without parl
amces_info_text_dauer_dicho_minu_parl<- cj((con1_minu_parl), chosen ~ feature1+feature2+ feature3+feature4+feature6, id = ~ID, estimate = "amce", 
                                           weights = ~scale_base_weights_f,
                                           by = ~info_text_dauer_dicho, feature_labels = feature_names,
                                           feature_order = c("feature6","feature4", "feature3", "feature2", "feature1"))
diff_amces_info_text_dauer_dicho_minu_parl<- cj((con1_minu_parl), chosen ~ feature1 +feature2+ feature3 +feature4+feature6, id = ~ID, estimate = "amce_diff",
                                                weights = ~scale_base_weights_f,
                                                by = ~info_text_dauer_dicho, feature_labels=feature_names,
                                                feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1")) 

#plot conditional amces time information 
plot(rbind(amces_info_text_dauer_dicho_minu_parl, diff_amces_info_text_dauer_dicho_minu_parl), feature_headers = FALSE) +ggplot2::facet_wrap(~BY, ncol = 3L)+
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 13),
        axis.text.x = element_text(colour = "black", size = 13),
        axis.title = element_text(colour = "black", size = 13),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 15),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#save conditional amces time information
#ggsave("subgroup_time_info_amces_without_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1.5, bg="white") 








# Subgroup Overview Robustness Checks (AMCEs) -----------------------------

#### BASE 
#prepare 
diff_amce_attention_choice <- cj(con1, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~attention)
diff_amce_timeinfo_choice <- cj(con1, amces_subgroup_choice, id = ~ID, estimate="amce_diff", weights = ~scale_base_weights_f,feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_amce_timeconjoint_choice <- cj(con1, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)

diff_estimates_choice_robust <-bind_rows(
  diff_amce_attention_choice %>% mutate(model = "Attentive"),
  diff_amce_timeinfo_choice %>% mutate(model = "Time Info."),
  diff_amce_timeconjoint_choice %>% mutate(model = "Time Conj."))


diff_choice_robust <- diff_estimates_choice_robust %>% dplyr::select(BY:level, estimate:model)
diff_choice_robust$model <- factor(diff_choice_robust$model, levels = c("Attentive",  "Time Info.",  "Time Conj." ))


#plot overview robustness checks
robust_diff<-plot(diff_choice_robust, feature_headers = FALSE, xlab = "Estimated Difference (AMCE) for different groups of adolescents\nRobustness Checks") + 
  ggplot2::facet_wrap(~model, ncol = 3) +
  theme(
    panel.border = element_rect(linetype = "dotted", colour = "black"),
    panel.background = element_blank(), 
    axis.text.y = element_text(colour = "black", size = 13),
    axis.text.x = element_text(colour = "black", size = 13),
    axis.title = element_text(colour = "black", size = 13),
    panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
    strip.background = element_blank(),
    strip.text.x = element_text(colour = "black", size = 15),
    legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

robust_diff


#save overview robustness checks
#ggsave("robustness_overview_differences_amces.png", width = 22, height=19, units = "cm", dpi=800, scale=1, bg="white") 


#### WITHOUT INCLUSION PARLIAMENT
#prepare 
diff_amce_attention_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~attention)
diff_amce_timeinfo_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_choice, id = ~ID, estimate="amce_diff", weights = ~scale_base_weights_f,feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_amce_timeconjoint_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)

diff_estimates_choice_robust_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_attention_choice_minu_inclusion_parl_rand %>% mutate(model = "Attentive"),
  diff_amce_timeinfo_choice_minu_inclusion_parl_rand %>% mutate(model = "Time Info."),
  diff_amce_timeconjoint_choice_minu_inclusion_parl_rand %>% mutate(model = "Time Conj."))


diff_choice_robust_minu_inclusion_parl_rand <- diff_estimates_choice_robust_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_choice_robust_minu_inclusion_parl_rand$model <- factor(diff_choice_robust_minu_inclusion_parl_rand$model, levels = c("Attentive",  "Time Info.",  "Time Conj." ))


#plot overview robustness checks
robust_diff_minu_inclusion_parl_rand<-plot(diff_choice_robust_minu_inclusion_parl_rand, feature_headers = FALSE, xlab = "Estimated Difference (AMCE) for different groups of adolescents\nRobustness Checks") + 
  ggplot2::facet_wrap(~model, ncol = 3) +
  theme(
    panel.border = element_rect(linetype = "dotted", colour = "black"),
    panel.background = element_blank(), 
    axis.text.y = element_text(colour = "black", size = 13),
    axis.text.x = element_text(colour = "black", size = 13),
    axis.title = element_text(colour = "black", size = 13),
    panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
    strip.background = element_blank(),
    strip.text.x = element_text(colour = "black", size = 15),
    legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

robust_diff_minu_inclusion_parl_rand


#save overview robustness checks
#ggsave("robustness_overview_differences_amces_without_inclusion_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1, bg="white") 


#### WITHOUT PARLIAMENT
#prepare 
diff_amce_attention_choice_minu_parl <- cj(con1_minu_parl, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~attention)
diff_amce_timeinfo_choice_minu_parl <- cj(con1_minu_parl, amces_subgroup_choice, id = ~ID, estimate="amce_diff", weights = ~scale_base_weights_f,feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_amce_timeconjoint_choice_minu_parl <- cj(con1_minu_parl, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)

diff_estimates_choice_robust_minu_parl <-bind_rows(
  diff_amce_attention_choice_minu_parl %>% mutate(model = "Attentive"),
  diff_amce_timeinfo_choice_minu_parl %>% mutate(model = "Time Info."),
  diff_amce_timeconjoint_choice_minu_parl %>% mutate(model = "Time Conj."))


diff_choice_robust_minu_parl <- diff_estimates_choice_robust_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_choice_robust_minu_parl$model <- factor(diff_choice_robust_minu_parl$model, levels = c("Attentive",  "Time Info.",  "Time Conj." ))


#plot overview robustness checks
robust_diff_minu_parl<-plot(diff_choice_robust_minu_parl, feature_headers = FALSE, xlab = "Estimated Difference (AMCE) for different groups of adolescents\nRobustness Checks") + 
  ggplot2::facet_wrap(~model, ncol = 3) +
  theme(
    panel.border = element_rect(linetype = "dotted", colour = "black"),
    panel.background = element_blank(), 
    axis.text.y = element_text(colour = "black", size = 13),
    axis.text.x = element_text(colour = "black", size = 13),
    axis.title = element_text(colour = "black", size = 13),
    panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
    strip.background = element_blank(),
    strip.text.x = element_text(colour = "black", size = 15),
    legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

robust_diff_minu_parl


#save overview robustness checks
#ggsave("robustness_overview_differences_amces_without_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1, bg="white") 







# MARGINAL MEANS ----------------------------------------------------------



# Baseline Marginal Means -------------------------------------------------
#baseline model
mm <- chosen ~ feature1 * feature4 + feature2 + feature3   +feature6
mm_choice <- cj(con1, mm, id = ~ID,  estimate="mm",weights = ~scale_base_weights_f, feature_labels = feature_names)
mm_choice

#plot baseline model
mm_plot<-plot(mm_choice,xlim = c(0.4,0.6), vline=0.5)+  
  theme(axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        axis.ticks = element_blank(), 
        legend.position = "none",
        panel.spacing.y = unit(1.5, "lines"),
        panel.border = element_rect(linetype = "blank"), 
        strip.background= element_rect(size = 2, colour="grey89", fill = "white"),
        strip.text.y.left = element_text(angle = 0),
        strip.placement = "outside") + 
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_brewer(palette="Dark2")

mm_plot

#save baseline model
#ggsave("mm_basemodel.png", width = 20, height=17, units = "cm", dpi=800, scale=1, bg="white") 






# Subgroup Marginal Means Characteristics -------------------------------------------------

f1 <- chosen ~ feature1 * feature4+feature2+ feature3 +feature6
f2<-chosen ~ feature1 + feature2 + feature3 + feature4 + feature6


# Dissatisfaction (MMs) ---------------------------------------------------
    # all - no subsets by constraint ------------------------------------------


#### BASE 

#prepare
con1$Dissatisfaction <- con1$demo_satis_dicho_rec
dissatisfaction_mm <- cj((con1), f1, id = ~ID, estimate = "mm", weights = ~scale_base_weights_f,by = ~Dissatisfaction, feature_labels=feature_names,
                         feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))
#plot
plot(dissatisfaction_mm, group="Dissatisfaction", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )

#save
#ggsave("dissatisfaction_mm.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 



    # subsets by constraint ---------------------------------------------------

#### WITHOUT INCLUSION PARL
#prepare
con1_minu_inclusion_parl_rand$Dissatisfaction <- con1_minu_inclusion_parl_rand$demo_satis_dicho_rec
dissatisfaction_mm <- cj((con1_minu_inclusion_parl_rand), f2, id = ~ID, estimate = "mm", weights = ~scale_base_weights_f,by = ~Dissatisfaction, feature_labels=feature_names,
                         feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))
#plot
plot(dissatisfaction_mm, group="Dissatisfaction", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )

#save
#ggsave("dissatisfaction_mm_without_inclusion_parl.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 



#### WITHOUT PARL
#prepare
con1_minu_parl$Dissatisfaction <- con1_minu_parl$demo_satis_dicho_rec
dissatisfaction_mm <- cj((con1_minu_parl), f2, id = ~ID, estimate = "mm", weights = ~scale_base_weights_f,by = ~Dissatisfaction, feature_labels=feature_names,
                         feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))
#plot
plot(dissatisfaction_mm, group="Dissatisfaction", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )

#save
#ggsave("dissatisfaction_mm_without_parl.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 









# Political Interest (MMs) ------------------------------------------------


    # all - no subsets by constraint ------------------------------------------
#### BASE

#prepare
con1$Interest<-con1$pol_int_dicho
interest_mm <- cj((con1), f1, id = ~ID, estimate = "mm", weights = ~scale_base_weights_f,by = ~Interest, feature_labels=feature_names,
                  feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))

#plot
plot(interest_mm, group="Interest", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )


#save
#ggsave("pol_int_mm.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 



    # subsets by constraint ---------------------------------------------------


#### WITHOUT INCLUSION PARLIAMENT
#prepare
con1_minu_inclusion_parl_rand$Interest<-con1_minu_inclusion_parl_rand$pol_int_dicho
interest_mm <- cj((con1_minu_inclusion_parl_rand), f2, id = ~ID, estimate = "mm", weights = ~scale_base_weights_f,by = ~Interest, feature_labels=feature_names,
                  feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))

#plot
plot(interest_mm, group="Interest", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )


#save
#ggsave("pol_int_mm_without_inclusion_parl.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 


#### WITHOUT PARLIAMENT
#prepare
con1_minu_parl$Interest<-con1_minu_parl$pol_int_dicho
interest_mm <- cj((con1_minu_parl), f2, id = ~ID, estimate = "mm", weights = ~scale_base_weights_f,by = ~Interest, feature_labels=feature_names,
                  feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))

#plot
plot(interest_mm, group="Interest", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )


#save
#ggsave("pol_int_mm_without_parl.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 






# Education (MMs) ---------------------------------------------------------


    # all - no subsets by constraint ------------------------------------------


#### BASE 

#prepare
con1$Education<-con1$gym


education_mm <- cj((con1), f1, id = ~ID, estimate = "mm", by = ~Education, feature_labels=feature_names,
                   weights = ~scale_base_weights_f,
                   feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))


#plot
plot(education_mm, group="Education", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )


#save
#ggsave("high_edu_mm_final.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 


    # subsets by constraint ---------------------------------------------------


#### WITHOUT INCLUSION PARLIAMENT
#prepare
con1_minu_inclusion_parl_rand$Education<-con1_minu_inclusion_parl_rand$gym


education_mm <- cj((con1_minu_inclusion_parl_rand), f2, id = ~ID, estimate = "mm", by = ~Education, feature_labels=feature_names,
                   weights = ~scale_base_weights_f,
                   feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))


#plot
plot(education_mm, group="Education", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )


#save
#ggsave("high_edu_mm_without_inclusion_parl_final.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 


#### WITHOUT PARLIAMENT
#prepare
con1_minu_parl$Education<-con1_minu_parl$gym


education_mm <- cj((con1_minu_parl), f2, id = ~ID, estimate = "mm", by = ~Education, feature_labels=feature_names,
                   weights = ~scale_base_weights_f,
                   feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))


#plot
plot(education_mm, group="Education", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )


#save
#ggsave("high_edu_mm_without__parl_final.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 







# Subgroup Marginal Means Robustness Checks -------------------------------

# Attention (MMs) ---------------------------------------------------

    # all - no subsets by constraint ------------------------------------------
#### BASE

#prepare
con1$Attention <- con1$attention
attention_mm <- cj((con1), f1, id = ~ID, estimate = "mm",weights = ~scale_base_weights_f, by = ~Attention, feature_labels=feature_names,
                   feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))

#plot
plot(attention_mm, group="Attention", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),)


#save
#ggsave("attention_mm.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 


    # subsets by constraint ---------------------------------------------------
#### WITHOUT INCLUSION PARL

#prepare
con1_minu_inclusion_parl_rand$Attention <- con1_minu_inclusion_parl_rand$attention
attention_mm <- cj((con1_minu_inclusion_parl_rand), f2, id = ~ID, estimate = "mm",weights = ~scale_base_weights_f, by = ~Attention, feature_labels=feature_names,
                   feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))

#plot
plot(attention_mm, group="Attention", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),)


#save
#ggsave("attention_mm_without_inclusion_parl.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 


#### WITHOUT INCLUSION PARL

#prepare
con1_minu_parl$Attention <- con1_minu_parl$attention
attention_mm <- cj((con1_minu_parl), f2, id = ~ID, estimate = "mm",weights = ~scale_base_weights_f, by = ~Attention, feature_labels=feature_names,
                   feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))

#plot
plot(attention_mm, group="Attention", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),)


#save
#ggsave("attention_mm_without_parl.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 






# Information Time (MMs) --------------------------------------------------------
    # all - no subsets by constraint ------------------------------------------
#### BASE

#prepare
con1$Information<-con1$info_text_dauer_dicho


infotime_mm <- cj((con1), f1, id = ~ID, estimate = "mm", weights = ~scale_base_weights_f,by = ~Information, feature_labels=feature_names,
                  feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))

#plot
plot(infotime_mm, group="Information", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )



#save
#ggsave("infotime_mm.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 


    # subsets by constraint ---------------------------------------------------


#### WITHOUT INCLUSION PARL

#prepare
con1_minu_inclusion_parl_rand$Information<-con1_minu_inclusion_parl_rand$info_text_dauer_dicho


infotime_mm <- cj((con1_minu_inclusion_parl_rand), f2, id = ~ID, estimate = "mm", weights = ~scale_base_weights_f,by = ~Information, feature_labels=feature_names,
                  feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))

#plot
plot(infotime_mm, group="Information", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )



#save
#ggsave("infotime_mm_without_inclusion_parl.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 


#### WITHOUT PARL

#prepare
con1_minu_parl$Information<-con1_minu_parl$info_text_dauer_dicho


infotime_mm <- cj((con1_minu_parl), f2, id = ~ID, estimate = "mm", weights = ~scale_base_weights_f,by = ~Information, feature_labels=feature_names,
                  feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))

#plot
plot(infotime_mm, group="Information", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )

#save
#ggsave("infotime_mm_without_parl.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 







# Conjoint Time (MMs) -----------------------------------------------------------
    # all - no subsets by constraint ------------------------------------------
#### BASE

#prepare
con1$Conjoint<-con1$timing_conjoint_num

conjointtime_mm <- cj((con1), f1, id = ~ID, estimate = "mm", weights = ~scale_base_weights_f,by = ~Conjoint, feature_labels=feature_names,
                      feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))

#plot
plot(conjointtime_mm, group="Conjoint", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )

#save
#ggsave("conjointtime_mm.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 


    # subsets by constraint ---------------------------------------------------


#### WITHOUT INCLUSION PARL

#prepare
con1_minu_inclusion_parl_rand$Conjoint<-con1_minu_inclusion_parl_rand$timing_conjoint_num

conjointtime_mm <- cj((con1_minu_inclusion_parl_rand), f2, id = ~ID, estimate = "mm", weights = ~scale_base_weights_f,by = ~Conjoint, feature_labels=feature_names,
                      feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))

#plot
plot(conjointtime_mm, group="Conjoint", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )

#save
ggsave("conjointtime_mm_without_inclusion_parl.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 


#### WITHOUT PARL

#prepare
con1_minu_parl$Conjoint<-con1_minu_parl$timing_conjoint_num

conjointtime_mm <- cj((con1_minu_parl), f2, id = ~ID, estimate = "mm", weights = ~scale_base_weights_f,by = ~Conjoint, feature_labels=feature_names,
                      feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"))

#plot
plot(conjointtime_mm, group="Conjoint", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )

#save
#ggsave("conjointtime_mm_without_parl.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 










# Winner vs. Loser (MMs) --------------------------------------------------
f3<- chosen ~ feature1*feature4+feature2+ feature3
f4<- chosen ~ feature1+ feature2+feature3+feature4

    # all - no subsets by constraint ------------------------------------------


#### BASE

#prepare
con1$Outcome_Favorability<-con1$feature6


of_mm <- cj((con1), f3, id = ~ID, estimate = "mm",weights = ~scale_base_weights_f, by = ~Outcome_Favorability, feature_labels=feature_names,
            feature_order = c("feature4", "feature3", "feature2", "feature1"))

#plot
plot(of_mm, group="Outcome_Favorability", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )


#save
#ggsave("outcome_fav_mm.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 




    # subsets by constraint ---------------------------------------------------

#### WITHOUT INCLUSION PARL

#prepare
con1_minu_inclusion_parl_rand$Outcome_Favorability<-con1_minu_inclusion_parl_rand$feature6


of_mm <- cj((con1_minu_inclusion_parl_rand), f4, id = ~ID, estimate = "mm",weights = ~scale_base_weights_f, by = ~Outcome_Favorability, feature_labels=feature_names,
            feature_order = c("feature4", "feature3", "feature2", "feature1"))

#plot
plot(of_mm, group="Outcome_Favorability", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )


#save
#ggsave("outcome_fav_mm_without_inclusion_parl.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 



#### WITHOUT PARL

#prepare
con1_minu_parl$Outcome_Favorability<-con1_minu_parl$feature6


of_mm <- cj((con1_minu_parl), f4, id = ~ID, estimate = "mm",weights = ~scale_base_weights_f, by = ~Outcome_Favorability, feature_labels=feature_names,
            feature_order = c("feature4", "feature3", "feature2", "feature1"))

#plot
plot(of_mm, group="Outcome_Favorability", vline=0.5, feature_headers = FALSE)+ggplot2::scale_color_manual(values=c("#1b9e77", "#e7298a"),na.translate=FALSE) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "#515151"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        legend.title = element_text( size=12), legend.text=element_text(size=12),
  )


#save
#ggsave("outcome_fav_mm_without_parl.png", width = 20, height=18, units = "cm", dpi=800, scale=1, bg="white") 








# RATING OUTCOME ----------------------------------------------------------


# AMCEs Rating Basemodel -------------------------------------------------------------------
#prepare
amces_rating_base <- rating ~ feature1 * feature4 + feature2 + feature3   +feature6
amce_rating <- cj(con1, amces_rating_base, id = ~ID,  estimate="amce",weights = ~scale_base_weights_f, feature_labels = feature_names)
amce_rating

#plot
amce_baseline_rating_plot<-plot(amce_rating,xlim = c(-0.4,0.4))+  
  theme(axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        axis.ticks = element_blank(), 
        legend.position = "none",
        panel.spacing.y = unit(1.5, "lines"),
        panel.border = element_rect(linetype = "blank"), 
        strip.background= element_rect(size = 2, colour="grey89", fill = "white"),
        strip.text.y.left = element_text(angle = 0),
        strip.placement = "outside") + 
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_brewer(palette="Dark2")

amce_baseline_rating_plot

#save
#ggsave("amce_rating_basemodel.png", width = 20, height=17, units = "cm", dpi=800, scale=1, bg="white") 






# Subgroup Analysis Rating  ------------------------------------------------------

amces_subgroup_rating <- rating ~ feature1 + feature2 + feature3  +feature4 +feature6




# Subgroup Overview Characteristics (AMCEs) -----------------------------------------------

    # all - no subsets by constraint ------------------------------------------
#### BASE

diff_amce_gym_rating <- cj(con1, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~gym)
diff_amce_interest_rating <- cj(con1, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~pol_int_dicho)
diff_amce_dissatisfied_rating <- cj(con1, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~demo_satis_dicho_rec)


diff_estimates_rating <-bind_rows(
  diff_amce_dissatisfied_rating %>% mutate(model = "Dissatisfied"),
  diff_amce_interest_rating %>% mutate(model = "Politically interested"),
  diff_amce_gym_rating %>% mutate(model = "Education high"))

diff_rating <- diff_estimates_rating %>% dplyr::select(BY:level, estimate:model)
diff_rating$model <- factor(diff_rating$model, levels = c("Dissatisfied",  "Politically interested","Education high"))


#plot overview
plot(diff_rating, feature_headers = FALSE, xlab = "\nEstimated Difference (AMCE) for different groups of adolescents\nRating Question") + 
  ggplot2::facet_wrap(~model, ncol = 3) +
  theme(
    panel.border = element_rect(linetype = "dotted", colour = "black"),
    panel.background = element_blank(), 
    axis.text.y = element_text(colour = "black", size = 12),
    axis.text.x = element_text(colour = "black", size = 12),
    axis.title = element_text(colour = "black", size = 12),
    panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
    strip.background = element_blank(),
    strip.text.x = element_text(colour = "black", size = 15),
    legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))


#ggsave("subgroup_overview_differences_amces_rating.png", width = 22, height=19, units = "cm", dpi=800, scale=1, bg="white") 



    # subsets by constraint ---------------------------------------------------
#### WITHOUT INCLUSION PARL

diff_amce_gym_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~gym)
diff_amce_interest_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~pol_int_dicho)
diff_amce_dissatisfied_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~demo_satis_dicho_rec)


diff_estimates_rating_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_dissatisfied_rating_minu_inclusion_parl_rand %>% mutate(model = "Dissatisfied"),
  diff_amce_interest_rating_minu_inclusion_parl_rand %>% mutate(model = "Politically interested"),
  diff_amce_gym_rating_minu_inclusion_parl_rand %>% mutate(model = "Education high"))

diff_rating_minu_inclusion_parl_rand <- diff_estimates_rating_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_rating_minu_inclusion_parl_rand$model <- factor(diff_rating_minu_inclusion_parl_rand$model, levels = c("Dissatisfied",  "Politically interested","Education high"))


#plot overview
plot(diff_rating_minu_inclusion_parl_rand, feature_headers = FALSE, xlab = "\nEstimated Difference (AMCE) for different groups of adolescents\nRating Question") + 
  ggplot2::facet_wrap(~model, ncol = 3) +
  theme(
    panel.border = element_rect(linetype = "dotted", colour = "black"),
    panel.background = element_blank(), 
    axis.text.y = element_text(colour = "black", size = 12),
    axis.text.x = element_text(colour = "black", size = 12),
    axis.title = element_text(colour = "black", size = 12),
    panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
    strip.background = element_blank(),
    strip.text.x = element_text(colour = "black", size = 15),
    legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))


#ggsave("subgroup_overview_differences_amces_rating_without_inclusion_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1, bg="white") 


#### WITHOUT PARL

diff_amce_gym_rating_minu_parl <- cj(con1_minu_parl, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~gym)
diff_amce_interest_rating_minu_parl <- cj(con1_minu_parl, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~pol_int_dicho)
diff_amce_dissatisfied_rating_minu_parl <- cj(con1_minu_parl, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~demo_satis_dicho_rec)


diff_estimates_rating_minu_parl <-bind_rows(
  diff_amce_dissatisfied_rating_minu_parl %>% mutate(model = "Dissatisfied"),
  diff_amce_interest_rating_minu_parl %>% mutate(model = "Politically interested"),
  diff_amce_gym_rating_minu_parl %>% mutate(model = "Education high"))

diff_rating_minu_parl <- diff_estimates_rating_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_rating_minu_parl$model <- factor(diff_rating_minu_parl$model, levels = c("Dissatisfied",  "Politically interested","Education high"))


#plot overview
plot(diff_rating_minu_parl, feature_headers = FALSE, xlab = "\nEstimated Difference (AMCE) for different groups of adolescents\nRating Question") + 
  ggplot2::facet_wrap(~model, ncol = 3) +
  theme(
    panel.border = element_rect(linetype = "dotted", colour = "black"),
    panel.background = element_blank(), 
    axis.text.y = element_text(colour = "black", size = 12),
    axis.text.x = element_text(colour = "black", size = 12),
    axis.title = element_text(colour = "black", size = 12),
    panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
    strip.background = element_blank(),
    strip.text.x = element_text(colour = "black", size = 15),
    legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))


#ggsave("subgroup_overview_differences_amces_rating_without_parl_oldweights.png", width = 22, height=19, units = "cm", dpi=800, scale=1, bg="white") 



# Subgroup Overview Robustness Checks (AMCEs) -----------------------------


    # all - no subsets by constraint ------------------------------------------
#### BASE

diff_amce_attention_rating <- cj(con1, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~attention)
diff_amce_timeinfo_rating <- cj(con1, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_amce_timeconjoint_rating <- cj(con1, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)

diff_estimates_rating_robust <-bind_rows(
  diff_amce_attention_rating %>% mutate(model = "Attentive"),
  diff_amce_timeinfo_rating %>% mutate(model = "Time Info."),
  diff_amce_timeconjoint_rating %>% mutate(model = "Time Conj."))


diff_rating_robust <- diff_estimates_rating_robust %>% dplyr::select(BY:level, estimate:model)
diff_rating_robust$model <- factor(diff_rating_robust$model, levels = c("Attentive",  "Time Info.",  "Time Conj." ))

#plot overview
plot(diff_rating_robust, feature_headers = FALSE, xlab = "\nEstimated Difference (AMCE) Robustness Checks\nRating Question") + 
  ggplot2::facet_wrap(~model, ncol = 3) +
  theme(
    panel.border = element_rect(linetype = "dotted", colour = "black"),
    panel.background = element_blank(), 
    axis.text.y = element_text(colour = "black", size = 12),
    axis.text.x = element_text(colour = "black", size = 12),
    axis.title = element_text(colour = "black", size = 12),
    panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
    strip.background = element_blank(),
    strip.text.x = element_text(colour = "black", size = 15),
    legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#ggsave("robustness_overview_differences_amces_rating.png", width = 22, height=19, units = "cm", dpi=800, scale=1, bg="white") 


    # subsets by constraint ---------------------------------------------------
#### WITHOUT INCLUSION PARL
diff_amce_attention_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~attention)
diff_amce_timeinfo_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_amce_timeconjoint_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)

diff_estimates_rating_robust_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_attention_rating_minu_inclusion_parl_rand %>% mutate(model = "Attentive"),
  diff_amce_timeinfo_rating_minu_inclusion_parl_rand %>% mutate(model = "Time Info."),
  diff_amce_timeconjoint_rating_minu_inclusion_parl_rand %>% mutate(model = "Time Conj."))


diff_rating_robust_minu_inclusion_parl_rand <- diff_estimates_rating_robust_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_rating_robust_minu_inclusion_parl_rand$model <- factor(diff_rating_robust_minu_inclusion_parl_rand$model, levels = c("Attentive",  "Time Info.",  "Time Conj." ))

#plot overview
plot(diff_rating_robust_minu_inclusion_parl_rand, feature_headers = FALSE, xlab = "\nEstimated Difference (AMCE) Robustness Checks\nRating Question") + 
  ggplot2::facet_wrap(~model, ncol = 3) +
  theme(
    panel.border = element_rect(linetype = "dotted", colour = "black"),
    panel.background = element_blank(), 
    axis.text.y = element_text(colour = "black", size = 12),
    axis.text.x = element_text(colour = "black", size = 12),
    axis.title = element_text(colour = "black", size = 12),
    panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
    strip.background = element_blank(),
    strip.text.x = element_text(colour = "black", size = 15),
    legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#ggsave("robustness_overview_differences_amces_rating_without_inclusion_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1, bg="white") 




#### WITHOUT PARL

diff_amce_attention_rating_minu_parl <- cj(con1_minu_parl, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~attention)
diff_amce_timeinfo_rating_minu_parl <- cj(con1_minu_parl, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_amce_timeconjoint_rating_minu_parl <- cj(con1_minu_parl, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)

diff_estimates_rating_robust_minu_parl <-bind_rows(
  diff_amce_attention_rating_minu_parl %>% mutate(model = "Attentive"),
  diff_amce_timeinfo_rating_minu_parl %>% mutate(model = "Time Info."),
  diff_amce_timeconjoint_rating_minu_parl %>% mutate(model = "Time Conj."))


diff_rating_robust_minu_parl <- diff_estimates_rating_robust_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_rating_robust_minu_parl$model <- factor(diff_rating_robust_minu_parl$model, levels = c("Attentive",  "Time Info.",  "Time Conj." ))

#plot overview
plot(diff_rating_robust_minu_parl, feature_headers = FALSE, xlab = "\nEstimated Difference (AMCE) Robustness Checks\nRating Question") + 
  ggplot2::facet_wrap(~model, ncol = 3) +
  theme(
    panel.border = element_rect(linetype = "dotted", colour = "black"),
    panel.background = element_blank(), 
    axis.text.y = element_text(colour = "black", size = 12),
    axis.text.x = element_text(colour = "black", size = 12),
    axis.title = element_text(colour = "black", size = 12),
    panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
    strip.background = element_blank(),
    strip.text.x = element_text(colour = "black", size = 15),
    legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("midnightblue", 9))

#ggsave("robustness_overview_differences_amces_rating_without_parl.png", width = 22, height=19, units = "cm", dpi=800, scale=1, bg="white") 





# winner loser difference -------------------------------------------------

    # all - no subsets by constraint ------------------------------------------
#### BASE

diff_winner_loser_rating <- cj((con1), rating ~ feature1 +feature2+ feature3+feature4, id = ~ID, estimate = "amce_diff", weights = ~scale_base_weights_f,
                               by = ~feature6, feature_labels=feature_names, feature_order = c("feature4", "feature3", "feature2", "feature1")) 



#Ploten Differences

winner_loser_diff_rating<-plot(diff_winner_loser_rating, feature_headers = FALSE, xlab = "Estimated Difference (AMCE) Winner vs Loser")  +
  theme(axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        axis.ticks = element_blank(), 
        legend.position = "none",
        panel.spacing.y = unit(1.5, "lines"),
        panel.border = element_rect(linetype = "blank"), 
        strip.background= element_rect(size = 2, colour="grey89", fill = "white"),
        strip.text.y.left = element_text(angle = 0),
        strip.placement = "outside") + 
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("darkblue", 9))+
  labs(title="Winner vs. Loser")
winner_loser_diff_rating


#ggsave("winner_loser_diff_rating.png", width = 15, height=16, units = "cm", dpi=800, scale=1, bg="white") 


    # subsets by constraint ---------------------------------------------------
#### WITHOUT INCLUSION PARL

diff_winner_loser_rating_minu_inclusion_parl_rand <- cj((con1_minu_inclusion_parl_rand), rating ~ feature1 +feature2+ feature3+feature4, id = ~ID, estimate = "amce_diff", weights = ~scale_base_weights_f,
                                                        by = ~feature6, feature_labels=feature_names, feature_order = c("feature4", "feature3", "feature2", "feature1")) 


diff_winner_loser_rating_minu_inclusion_parl_rand_plot<-plot(diff_winner_loser_rating_minu_inclusion_parl_rand, feature_headers = FALSE, xlab = "Estimated Difference (AMCE) Winner vs Loser")  +
  theme(axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        axis.ticks = element_blank(), 
        legend.position = "none",
        panel.spacing.y = unit(1.5, "lines"),
        panel.border = element_rect(linetype = "blank"), 
        strip.background= element_rect(size = 2, colour="grey89", fill = "white"),
        strip.text.y.left = element_text(angle = 0),
        strip.placement = "outside") + 
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("darkblue", 9))+
  labs(title="Winner vs. Loser")
diff_winner_loser_rating_minu_inclusion_parl_rand_plot


#ggsave("winner_loser_diff_rating_without_inclusion_parl.png", width = 15, height=16, units = "cm", dpi=800, scale=1, bg="white") 


#### WITHOUT PARL


diff_winner_loser_rating_minu_parl <- cj((con1_minu_parl), rating ~ feature1 +feature2+ feature3+feature4, id = ~ID, estimate = "amce_diff", weights = ~scale_base_weights_f,
                                         by = ~feature6, feature_labels=feature_names, feature_order = c("feature4", "feature3", "feature2", "feature1")) 


diff_winner_loser_rating_minu_parl_plot<-plot(diff_winner_loser_rating_minu_parl, feature_headers = FALSE, xlab = "Estimated Difference (AMCE) Winner vs Loser")  +
  theme(axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        axis.ticks = element_blank(), 
        legend.position = "none",
        panel.spacing.y = unit(1.5, "lines"),
        panel.border = element_rect(linetype = "blank"), 
        strip.background= element_rect(size = 2, colour="grey89", fill = "white"),
        strip.text.y.left = element_text(angle = 0),
        strip.placement = "outside") + 
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_manual(values = rep("darkblue", 9))+
  labs(title="Winner vs. Loser")
diff_winner_loser_rating_minu_parl_plot


#ggsave("winner_loser_diff_rating_without_parl.png", width = 15, height=16, units = "cm", dpi=800, scale=1, bg="white") 





# Baseline Marginal Means Rating -------------------------------------------------
mm_rating <- rating ~ feature1 * feature4 + feature2 + feature3   +feature6
mm_rating_1 <- cj(con1, mm, id = ~ID,  estimate="mm",weights = ~scale_base_weights_f, feature_labels = feature_names)
mm_rating_1


#plot baseline model
mm_plot_rating<-plot(mm_rating_1,xlim = c(0.4,0.6), vline=0.5)+  
  theme(axis.text.y = element_text(colour = "black", size = 12),
        axis.text.x = element_text(colour = "black", size = 12),
        axis.title = element_text(colour = "black", size = 12),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        axis.ticks = element_blank(), 
        legend.position = "none",
        panel.spacing.y = unit(1.5, "lines"),
        panel.border = element_rect(linetype = "blank"), 
        strip.background= element_rect(size = 2, colour="grey89", fill = "white"),
        strip.text.y.left = element_text(angle = 0),
        strip.placement = "outside") + 
  geom_linerange(aes(xmin=lower, xmax=upper), size = 0.8)+
  geom_point(size=2) +
  scale_color_brewer(palette="Dark2")

mm_plot_rating

#save baseline model
#ggsave("mm_basemodel_rating.png", width = 20, height=17, units = "cm", dpi=800, scale=1, bg="white") 






# Subgroups Übersicht AMCE & Marginal Means Choice and Rating -------------------------------

    # all - no subsets by constraint --------------------------------------------------------------------
#repeat relevant calculations


#AMCE choice
amces_subgroup_choice <- chosen ~ feature1 + feature2 + feature3  +feature4 +feature6
#base
diff_amce_gym_choice <- cj(con1, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~gym)
diff_amce_interest_choice <- cj(con1, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~pol_int_dicho)
diff_amce_dissatisfied_choice <- cj(con1, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~demo_satis_dicho_rec)

#MM choice
mm_subgroup <- chosen ~ feature1+ feature2 + feature3 +feature4 +feature6

diff_mm_gym_choice <- cj(con1, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~Education)
diff_mm_interest_choice <- cj(con1, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~Interest)
diff_mm_dissatisfied_choice <- cj(con1, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~Dissatisfaction)



#AMCE rating
amces_subgroup_rating <- rating ~ feature1 + feature2 + feature3  +feature4 +feature6
#### BASE
diff_amce_gym_rating <- cj(con1, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~gym)
diff_amce_interest_rating <- cj(con1, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~pol_int_dicho)
diff_amce_dissatisfied_rating <- cj(con1, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~demo_satis_dicho_rec)


#MM rating
mm_subgroup_rating <- rating ~ feature1 + feature2 + feature3 +feature4 +feature6

diff_mm_gym_rating <- cj(con1, mm_subgroup_rating, id = ~ID, estimate="mm_diff", weights = ~scale_base_weights_f,feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~Education)
diff_mm_interest_rating <- cj(con1, mm_subgroup_rating, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~Interest)
diff_mm_dissatisfied_rating <- cj(con1, mm_subgroup_rating, id = ~ID, estimate="mm_diff", weights = ~scale_base_weights_f,feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~Dissatisfaction)



#### DISSATISFACTION

#Dissatisfaction: AMCE choice; MM choice; AMCE rating; MM rating
#AMCE choice: diff_amce_dissatisfied_choice; MM choice: diff_mm_dissatisfied_choice
#AMCE rating: diff_amce_dissatisfied_rating; MM rating: diff_mm_dissatisfied_rating


overview_choice_satisfaction <-bind_rows(
  diff_amce_dissatisfied_choice %>% mutate(model = "AMCE (Choice)"),
  diff_mm_dissatisfied_choice %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_satis <- overview_choice_satisfaction %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_satis$model <- factor(diff_overview_choice_satis$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_satisfaction_plot<-plot(diff_overview_choice_satis, feature_headers = FALSE, xlab="", xlim = c(-0.2,0.2)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_choice_satisfaction_plot


overview_rating_satisfaction <-bind_rows(
  diff_amce_dissatisfied_rating %>% mutate(model = "AMCE (Rating)"),
  diff_mm_dissatisfied_rating %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_satis <- overview_rating_satisfaction %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_satis$model <- factor(diff_overview_rating_satis$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_satisfaction_plot<-plot(diff_overview_rating_satis, feature_headers = FALSE, xlab="",xlim = c(-0.4,0.4)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_satisfaction_plot




combined_satis <- ggarrange(overview_choice_satisfaction_plot + theme(axis.ticks.y = element_blank(),
                                                                      plot.margin = margin(r = 1) ),
                            overview_rating_satisfaction_plot + theme(
                              
                              axis.ticks.y = element_blank(),
                              axis.title.y = element_blank(),
                              axis.text.y=element_blank(),
                              plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_satis

annotate_figure(combined_satis, bottom = text_grob("Estimated differences for 'political dissatisfaction'\ncompared to 'political satisfaction'", 
                                                   color = "black", size = 24, hjust= 0.27))

ggsave("subgroup_overview_dissatisfaction_mm_amce_rating_choice.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 


#### POLITICAL INTEREST

#AMCE choice: diff_amce_interest_choice; MM choice: diff_mm_interest_choice
#AMCE rating: diff_amce_interest_rating; MM rating: diff_mm_interest_rating

overview_choice_interest <-bind_rows(
  diff_amce_interest_choice %>% mutate(model = "AMCE (Choice)"),
  diff_mm_interest_choice %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_interest <- overview_choice_interest %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_interest$model <- factor(diff_overview_choice_interest$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_interest_plot<-plot(diff_overview_choice_interest, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_choice_interest_plot


overview_rating_interest <-bind_rows(
  diff_amce_interest_rating %>% mutate(model = "AMCE (Rating)"),
  diff_mm_interest_rating %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_interest <- overview_rating_interest %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_interest$model <- factor(diff_overview_rating_interest$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_interest_plot<-plot(diff_overview_rating_interest, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_interest_plot




combined_interest <- ggarrange(overview_choice_interest_plot + theme(axis.ticks.y = element_blank(),
                                                                     plot.margin = margin(r = 1) ),
                               overview_rating_interest_plot + theme(
                                 axis.ticks.y = element_blank(),
                                 axis.title.y = element_blank(),
                                 axis.text.y=element_blank(),
                                 plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_interest

annotate_figure(combined_interest, bottom = text_grob("Estimated differences for 'high Interest'\ncompared to 'low Interest'", 
                                                      color = "black", size = 24, hjust= 0.27))


#ggsave("subgroup_overview_pol_int_mm_amce_rating_choice.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 

#### EDUCATION

#AMCE choice: diff_amce_gym_choice; MM choice: diff_mm_gym_choice
#AMCE rating: diff_amce_gym_rating; MM rating: diff_mm_gym_rating

overview_choice_gym <-bind_rows(
  diff_amce_gym_choice %>% mutate(model = "AMCE (Choice)"),
  diff_mm_gym_choice %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_gym <- overview_choice_gym %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_gym$model <- factor(diff_overview_choice_gym$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_gym_plot<-plot(diff_overview_choice_gym, feature_headers = FALSE, xlab="", xlim = c(-0.15,0.15)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_choice_gym_plot


overview_rating_gym <-bind_rows(
  diff_amce_gym_rating %>% mutate(model = "AMCE (Rating)"),
  diff_mm_gym_rating %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_gym <- overview_rating_gym %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_gym$model <- factor(diff_overview_rating_gym$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_gym_plot<-plot(diff_overview_rating_gym, feature_headers = FALSE, xlab="", xlim = c(-0.4,0.4)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_gym_plot




combined_gym <- ggarrange(overview_choice_gym_plot + theme(axis.ticks.y = element_blank(),
                                                           plot.margin = margin(r = 1) ),
                          overview_rating_gym_plot + theme(
                            axis.ticks.y = element_blank(),
                            axis.title.y = element_blank(),
                            axis.text.y=element_blank(),
                            plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_gym

annotate_figure(combined_gym, bottom = text_grob("Estimated differences for 'Education high'\ncompared to 'Education low'", 
                                                 color = "black", size = 24, hjust= 0.27))



ggsave("subgroup_overview_gym_mm_amce_rating_choice.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 




    # subsets by constraint --------------------------------------------------------------------

#without inclusion parliament
#repeat relevant calculations

#AMCEs choice
amces_subgroup_choice <- chosen ~ feature1 + feature2 + feature3  +feature4 +feature6
diff_amce_gym_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~gym)
diff_amce_interest_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~pol_int_dicho)
diff_amce_dissatisfied_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~demo_satis_dicho_rec)
diff_amce_gym_choice_minu_inclusion_parl_rand
#AMCEs rating
amces_subgroup_rating <- rating ~ feature1 + feature2 + feature3  +feature4 +feature6
diff_amce_gym_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~gym)
diff_amce_interest_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~pol_int_dicho)
diff_amce_dissatisfied_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~demo_satis_dicho_rec)


#MM choice
mm_subgroup <- chosen ~ feature1+ feature2 + feature3 +feature4 +feature6

diff_mm_gym_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~Education)
diff_mm_interest_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~Interest)
diff_mm_dissatisfied_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~Dissatisfaction)



#MM rating
mm_subgroup_rating <- rating ~ feature1 + feature2 + feature3 +feature4 +feature6

diff_mm_gym_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, mm_subgroup_rating, id = ~ID, estimate="mm_diff", weights = ~scale_base_weights_f,feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~Education)
diff_mm_interest_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, mm_subgroup_rating, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~Interest)
diff_mm_dissatisfied_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, mm_subgroup_rating, id = ~ID, estimate="mm_diff", weights = ~scale_base_weights_f,feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~Dissatisfaction)


#### DISSATISFACTION


overview_choice_satisfaction_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_dissatisfied_choice_minu_inclusion_parl_rand %>% mutate(model = "AMCE (Choice)"),
  diff_mm_dissatisfied_choice_minu_inclusion_parl_rand %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_satis_minu_inclusion_parl_rand <- overview_choice_satisfaction_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_satis_minu_inclusion_parl_rand$model <- factor(diff_overview_choice_satis_minu_inclusion_parl_rand$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_satisfaction_plot_minu_inclusion_parl_rand<-plot(diff_overview_choice_satis_minu_inclusion_parl_rand, feature_headers = FALSE, xlab="", xlim = c(-0.2,0.2)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_choice_satisfaction_plot_minu_inclusion_parl_rand


overview_rating_satisfaction_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_dissatisfied_rating_minu_inclusion_parl_rand %>% mutate(model = "AMCE (Rating)"),
  diff_mm_dissatisfied_rating_minu_inclusion_parl_rand %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_satis_minu_inclusion_parl_rand <- overview_rating_satisfaction_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_satis_minu_inclusion_parl_rand$model <- factor(diff_overview_rating_satis_minu_inclusion_parl_rand$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_satisfaction_plot_minu_inclusion_parl_rand<-plot(diff_overview_rating_satis_minu_inclusion_parl_rand, feature_headers = FALSE, xlab="",xlim = c(-0.4,0.4)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_satisfaction_plot_minu_inclusion_parl_rand




combined_satis_minu_inclusion_parl_rand <- ggarrange(overview_choice_satisfaction_plot_minu_inclusion_parl_rand + theme(axis.ticks.y = element_blank(),
                                                                                                                        plot.margin = margin(r = 1) ),
                                                     overview_rating_satisfaction_plot_minu_inclusion_parl_rand + theme(
                                                       
                                                       axis.ticks.y = element_blank(),
                                                       axis.title.y = element_blank(),
                                                       axis.text.y=element_blank(),
                                                       plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_satis_minu_inclusion_parl_rand

annotate_figure(combined_satis_minu_inclusion_parl_rand, bottom = text_grob("Estimated differences for 'political dissatisfaction'\ncompared to 'political satisfaction'", 
                                                                            color = "black", size = 24, hjust= 0.27))

#ggsave("subgroup_overview_dissatisfaction_mm_amce_rating_choice_without_inclusion_parl.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 


#### POL INTEREST

overview_choice_interest_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_interest_choice_minu_inclusion_parl_rand %>% mutate(model = "AMCE (Choice)"),
  diff_mm_interest_choice_minu_inclusion_parl_rand %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_interest_minu_inclusion_parl_rand <- overview_choice_interest_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_interest_minu_inclusion_parl_rand$model <- factor(diff_overview_choice_interest_minu_inclusion_parl_rand$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_interest_plot_minu_inclusion_parl_rand<-plot(diff_overview_choice_interest_minu_inclusion_parl_rand, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_choice_interest_plot_minu_inclusion_parl_rand


overview_rating_interest_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_interest_rating_minu_inclusion_parl_rand %>% mutate(model = "AMCE (Rating)"),
  diff_mm_interest_rating_minu_inclusion_parl_rand %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_interest_minu_inclusion_parl_rand <- overview_rating_interest_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_interest_minu_inclusion_parl_rand$model <- factor(diff_overview_rating_interest_minu_inclusion_parl_rand$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_interest_plot_minu_inclusion_parl_rand<-plot(diff_overview_rating_interest_minu_inclusion_parl_rand, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_interest_plot_minu_inclusion_parl_rand




combined_interest_minu_inclusion_parl_rand <- ggarrange(overview_choice_interest_plot_minu_inclusion_parl_rand + theme(axis.ticks.y = element_blank(),
                                                                                                                       plot.margin = margin(r = 1) ),
                                                        overview_rating_interest_plot_minu_inclusion_parl_rand + theme(
                                                          axis.ticks.y = element_blank(),
                                                          axis.title.y = element_blank(),
                                                          axis.text.y=element_blank(),
                                                          plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_interest_minu_inclusion_parl_rand

annotate_figure(combined_interest_minu_inclusion_parl_rand, bottom = text_grob("Estimated differences for 'high Interest'\ncompared to 'low Interest'", 
                                                                               color = "black", size = 24, hjust= 0.27))


#ggsave("subgroup_overview_pol_int_mm_amce_rating_choice_without_inclusion_parl.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 

#### EDUCATION


overview_choice_gym_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_gym_choice_minu_inclusion_parl_rand %>% mutate(model = "AMCE (Choice)"),
  diff_mm_gym_choice_minu_inclusion_parl_rand %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_gym_minu_inclusion_parl_rand <- overview_choice_gym_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_gym_minu_inclusion_parl_rand$model <- factor(diff_overview_choice_gym_minu_inclusion_parl_rand$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_gym_plot_minu_inclusion_parl_rand<-plot(diff_overview_choice_gym_minu_inclusion_parl_rand, feature_headers = FALSE, xlab="", xlim = c(-0.15,0.15)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_choice_gym_plot_minu_inclusion_parl_rand


overview_rating_gym_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_gym_rating_minu_inclusion_parl_rand %>% mutate(model = "AMCE (Rating)"),
  diff_mm_gym_rating_minu_inclusion_parl_rand %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_gym_minu_inclusion_parl_rand <- overview_rating_gym_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_gym_minu_inclusion_parl_rand$model <- factor(diff_overview_rating_gym_minu_inclusion_parl_rand$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_gym_plot_minu_inclusion_parl_rand<-plot(diff_overview_rating_gym_minu_inclusion_parl_rand, feature_headers = FALSE, xlab="", xlim = c(-0.4,0.4)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_gym_plot_minu_inclusion_parl_rand




combined_gym_minu_inclusion_parl_rand <- ggarrange(overview_choice_gym_plot_minu_inclusion_parl_rand + theme(axis.ticks.y = element_blank(),
                                                                                                             plot.margin = margin(r = 1) ),
                                                   overview_rating_gym_plot_minu_inclusion_parl_rand + theme(
                                                     axis.ticks.y = element_blank(),
                                                     axis.title.y = element_blank(),
                                                     axis.text.y=element_blank(),
                                                     plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_gym_minu_inclusion_parl_rand

annotate_figure(combined_gym_minu_inclusion_parl_rand, bottom = text_grob("Estimated differences for 'Education high'\ncompared to 'Education low'", 
                                                                          color = "black", size = 24, hjust= 0.27))



ggsave("subgroup_overview_gym_mm_amce_rating_choice_without_inclusion_parl.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 






#without parliament
#repeat relevant calculations

#AMCEs choice
amces_subgroup_choice <- chosen ~ feature1 + feature2 + feature3  +feature4 +feature6
diff_amce_gym_choice_minu_parl <- cj(con1_minu_parl, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~gym)
diff_amce_interest_choice_minu_parl <- cj(con1_minu_parl, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~pol_int_dicho)
diff_amce_dissatisfied_choice_minu_parl <- cj(con1_minu_parl, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~demo_satis_dicho_rec)

#AMCEs rating
amces_subgroup_rating <- rating ~ feature1 + feature2 + feature3  +feature4 +feature6
diff_amce_gym_rating_minu_parl <- cj(con1_minu_parl, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~gym)
diff_amce_interest_rating_minu_parl <- cj(con1_minu_parl, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~pol_int_dicho)
diff_amce_dissatisfied_rating_minu_parl <- cj(con1_minu_parl, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~demo_satis_dicho_rec)


#MM choice
mm_subgroup <- chosen ~ feature1+ feature2 + feature3 +feature4 +feature6

diff_mm_gym_choice_minu_parl <- cj(con1_minu_parl, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~Education)
diff_mm_interest_choice_minu_parl <- cj(con1_minu_parl, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~Interest)
diff_mm_dissatisfied_choice_minu_parl <- cj(con1_minu_parl, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~Dissatisfaction)



#MM rating
mm_subgroup_rating <- rating ~ feature1 + feature2 + feature3 +feature4 +feature6

diff_mm_gym_rating_minu_parl <- cj(con1_minu_parl, mm_subgroup_rating, id = ~ID, estimate="mm_diff", weights = ~scale_base_weights_f,feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~Education)
diff_mm_interest_rating_minu_parl <- cj(con1_minu_parl, mm_subgroup_rating, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~Interest)
diff_mm_dissatisfied_rating_minu_parl <- cj(con1_minu_parl, mm_subgroup_rating, id = ~ID, estimate="mm_diff", weights = ~scale_base_weights_f,feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~Dissatisfaction)


#### DISSATISFACTION


overview_choice_satisfaction_minu_parl <-bind_rows(
  diff_amce_dissatisfied_choice_minu_parl %>% mutate(model = "AMCE (Choice)"),
  diff_mm_dissatisfied_choice_minu_parl %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_satis_minu_parl <- overview_choice_satisfaction_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_satis_minu_parl$model <- factor(diff_overview_choice_satis_minu_parl$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_satisfaction_plot_minu_parl<-plot(diff_overview_choice_satis_minu_parl, feature_headers = FALSE, xlab="", xlim = c(-0.2,0.2)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_choice_satisfaction_plot_minu_parl


overview_rating_satisfaction_minu_parl <-bind_rows(
  diff_amce_dissatisfied_rating_minu_parl %>% mutate(model = "AMCE (Rating)"),
  diff_mm_dissatisfied_rating_minu_parl %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_satis_minu_parl <- overview_rating_satisfaction_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_satis_minu_parl$model <- factor(diff_overview_rating_satis_minu_parl$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_satisfaction_plot_minu_parl<-plot(diff_overview_rating_satis_minu_parl, feature_headers = FALSE, xlab="",xlim = c(-0.4,0.4)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_satisfaction_plot_minu_parl




combined_satis_minu_parl <- ggarrange(overview_choice_satisfaction_plot_minu_parl + theme(axis.ticks.y = element_blank(),
                                                                                          plot.margin = margin(r = 1) ),
                                      overview_rating_satisfaction_plot_minu_parl + theme(
                                        
                                        axis.ticks.y = element_blank(),
                                        axis.title.y = element_blank(),
                                        axis.text.y=element_blank(),
                                        plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_satis_minu_parl

annotate_figure(combined_satis_minu_parl, bottom = text_grob("Estimated differences for 'political dissatisfaction'\ncompared to 'political satisfaction'", 
                                                             color = "black", size = 24, hjust= 0.27))

#ggsave("subgroup_overview_dissatisfaction_mm_amce_rating_choice_without_parl.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 


#### POL INTEREST

overview_choice_interest_minu_parl <-bind_rows(
  diff_amce_interest_choice_minu_parl %>% mutate(model = "AMCE (Choice)"),
  diff_mm_interest_choice_minu_parl %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_interest_minu_parl <- overview_choice_interest_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_interest_minu_parl$model <- factor(diff_overview_choice_interest_minu_parl$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_interest_plot_minu_parl<-plot(diff_overview_choice_interest_minu_parl, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_choice_interest_plot_minu_parl


overview_rating_interest_minu_parl <-bind_rows(
  diff_amce_interest_rating_minu_parl %>% mutate(model = "AMCE (Rating)"),
  diff_mm_interest_rating_minu_parl %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_interest_minu_parl <- overview_rating_interest_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_interest_minu_parl$model <- factor(diff_overview_rating_interest_minu_parl$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_interest_plot_minu_parl<-plot(diff_overview_rating_interest_minu_parl, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_interest_plot_minu_parl




combined_interest_minu_parl <- ggarrange(overview_choice_interest_plot_minu_parl + theme(axis.ticks.y = element_blank(),
                                                                                         plot.margin = margin(r = 1) ),
                                         overview_rating_interest_plot_minu_parl + theme(
                                           axis.ticks.y = element_blank(),
                                           axis.title.y = element_blank(),
                                           axis.text.y=element_blank(),
                                           plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_interest_minu_parl

annotate_figure(combined_interest_minu_parl, bottom = text_grob("Estimated differences for 'high Interest'\ncompared to 'low Interest'", 
                                                                color = "black", size = 24, hjust= 0.27))


#ggsave("subgroup_overview_pol_int_mm_amce_rating_choice_without_parl_oldweights.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 

#### EDUCATION


overview_choice_gym_minu_parl <-bind_rows(
  diff_amce_gym_choice_minu_parl %>% mutate(model = "AMCE (Choice)"),
  diff_mm_gym_choice_minu_parl %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_gym_minu_parl <- overview_choice_gym_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_gym_minu_parl$model <- factor(diff_overview_choice_gym_minu_parl$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_gym_plot_minu_parl<-plot(diff_overview_choice_gym_minu_parl, feature_headers = FALSE, xlab="", xlim = c(-0.15,0.15)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_choice_gym_plot_minu_parl


overview_rating_gym_minu_parl <-bind_rows(
  diff_amce_gym_rating_minu_parl %>% mutate(model = "AMCE (Rating)"),
  diff_mm_gym_rating_minu_parl %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_gym_minu_parl <- overview_rating_gym_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_gym_minu_parl$model <- factor(diff_overview_rating_gym_minu_parl$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_gym_plot_minu_parl<-plot(diff_overview_rating_gym_minu_parl, feature_headers = FALSE, xlab="", xlim = c(-0.4,0.4)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_gym_plot_minu_parl




combined_gym_minu_parl <- ggarrange(overview_choice_gym_plot_minu_parl + theme(axis.ticks.y = element_blank(),
                                                                               plot.margin = margin(r = 1) ),
                                    overview_rating_gym_plot_minu_parl + theme(
                                      axis.ticks.y = element_blank(),
                                      axis.title.y = element_blank(),
                                      axis.text.y=element_blank(),
                                      plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_gym_minu_parl

annotate_figure(combined_gym_minu_parl, bottom = text_grob("Estimated differences for 'Education high'\ncompared to 'Education low'", 
                                                           color = "black", size = 24, hjust= 0.27))



ggsave("subgroup_overview_gym_mm_amce_rating_choice_without_parl.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 











# Robustness Check Übersicht AMCE & Marginal Means Choice and Rating ----------------------------------------------
    # all - no subsets by constraint --------------------------------------------------------------------

#calculate relevant again

#AMCEs choice
amces_subgroup_choice <- chosen ~ feature1 + feature2 + feature3  +feature4 +feature6
diff_amce_attention_choice <- cj(con1, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~attention)
diff_amce_timeinfo_choice <- cj(con1, amces_subgroup_choice, id = ~ID, estimate="amce_diff", weights = ~scale_base_weights_f,feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_amce_timeconjoint_choice <- cj(con1, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)


#AMCEs rating
amces_subgroup_rating <- rating ~ feature1 + feature2 + feature3  +feature4 +feature6
diff_amce_attention_rating <- cj(con1, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~attention)
diff_amce_timeinfo_rating <- cj(con1, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_amce_timeconjoint_rating <- cj(con1, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)

#MMs choice
mm_subgroup <- chosen ~ feature1+ feature2 + feature3 +feature4 +feature6
diff_mm_attention_choice <- cj(con1, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~attention)
diff_mm_timeinfo_choice <- cj(con1, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_mm_timeconjoint_choice <- cj(con1, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)


#MMs rating
mm_subgroup_rating <- rating ~ feature1+ feature2 + feature3 +feature4 +feature6
diff_mm_attention_rating <- cj(con1, mm_subgroup_rating, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~attention)
diff_mm_timeinfo_rating <- cj(con1, mm_subgroup_rating, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_mm_timeconjoint_rating <- cj(con1, mm_subgroup_rating, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)




#### ATTENTION

overview_choice_attention <-bind_rows(
  diff_amce_attention_choice %>% mutate(model = "AMCE (Choice)"),
  diff_mm_attention_choice %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_attention <- overview_choice_attention %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_attention$model <- factor(diff_overview_choice_attention$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_attention_plot<-plot(diff_overview_choice_attention, feature_headers = FALSE, xlab="", xlim = c(-0.15,0.15)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))
overview_choice_attention_plot


overview_rating_attention <-bind_rows(
  diff_amce_attention_rating %>% mutate(model = "AMCE (Rating)"),
  diff_mm_attention_rating %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_attention <- overview_rating_attention %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_attention$model <- factor(diff_overview_rating_attention$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_attention_plot<-plot(diff_overview_rating_attention, feature_headers = FALSE, xlab="", xlim = c(-0.3,0.3)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_attention_plot



combined_attention <- ggarrange(overview_choice_attention_plot + theme(axis.ticks.y = element_blank(),
                                                                       plot.margin = margin(r = 1) ),
                                overview_rating_attention_plot + theme(
                                  axis.ticks.y = element_blank(),
                                  axis.title.y = element_blank(),
                                  axis.text.y=element_blank(),
                                  plot.margin = margin(l = 1)), widths = c(1.5,1))

combined_attention

annotate_figure(combined_attention, bottom = text_grob("Estimated differences for 'high Attention'\ncompared to 'low Attention'", 
                                                       color = "black", size = 24, hjust= 0.27))


#ggsave("subgroup_overview_att_mm_amce_rating_choice.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 



#### TIMEINFO

overview_choice_timeinfo <-bind_rows(
  diff_amce_timeinfo_choice %>% mutate(model = "AMCE (Choice)"),
  diff_mm_timeinfo_choice %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_timeinfo <- overview_choice_timeinfo %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_timeinfo$model <- factor(diff_overview_choice_timeinfo$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_timeinfo_plot<-plot(diff_overview_choice_timeinfo, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_choice_timeinfo_plot


overview_rating_timeinfo <-bind_rows(
  diff_amce_timeinfo_rating %>% mutate(model = "AMCE (Rating)"),
  diff_mm_timeinfo_rating %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_timeinfo <- overview_rating_timeinfo %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_timeinfo$model <- factor(diff_overview_rating_timeinfo$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_timeinfo_plot<-plot(diff_overview_rating_timeinfo, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))
overview_rating_timeinfo_plot




combined_timeinfo <- ggarrange(overview_choice_timeinfo_plot + theme(axis.ticks.y = element_blank(),
                                                                     plot.margin = margin(r = 1) ),
                               overview_rating_timeinfo_plot + theme(
                                 axis.ticks.y = element_blank(),
                                 axis.title.y = element_blank(),
                                 axis.text.y=element_blank(),
                                 plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_timeinfo

annotate_figure(combined_timeinfo, bottom = text_grob("Estimated differences for 'more Time Info'\ncompared to 'less Time Info'", 
                                                      color = "black", size = 24, hjust= 0.27))


#ggsave("subgroup_overview_timeinfo_mm_amce_rating_choice.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 



#### TIME CONJOINT

overview_choice_timeconjoint <-bind_rows(
  diff_amce_timeconjoint_choice %>% mutate(model = "AMCE (Choice)"),
  diff_mm_timeconjoint_choice %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_timeconjoint <- overview_choice_timeconjoint %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_timeconjoint$model <- factor(diff_overview_choice_timeconjoint$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_timeconjoint_plot<-plot(diff_overview_choice_timeconjoint, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))


overview_choice_timeconjoint_plot



overview_rating_timeconjoint <-bind_rows(
  diff_amce_timeconjoint_rating %>% mutate(model = "AMCE (Rating)"),
  diff_mm_timeconjoint_rating %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_timeconjoint <- overview_rating_timeconjoint %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_timeconjoint$model <- factor(diff_overview_rating_timeconjoint$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_timeconjoint_plot<-plot(diff_overview_rating_timeconjoint, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))
overview_rating_timeconjoint_plot




combined_timeconjoint <- ggarrange(overview_choice_timeconjoint_plot + theme(axis.ticks.y = element_blank(),
                                                                             plot.margin = margin(r = 1) ),
                                   overview_rating_timeconjoint_plot + theme(
                                     axis.ticks.y = element_blank(),
                                     axis.title.y = element_blank(),
                                     axis.text.y=element_blank(),
                                     plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_timeconjoint

annotate_figure(combined_timeconjoint, bottom = text_grob("Estimated differences for 'more Time Conjoint'\ncompared to 'less Time Conjoint'", 
                                                          color = "black", size = 24, hjust= 0.27))


#ggsave("subgroup_overview_timeconjoint_mm_amce_rating_choice.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 



    # subsets by constraint --------------------------------------------------

#without inclusion parliament
#calculate relevant again

#AMCEs choice
amces_subgroup_choice <- chosen ~ feature1 + feature2 + feature3  +feature4 +feature6
diff_amce_attention_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~attention)
diff_amce_timeinfo_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_choice, id = ~ID, estimate="amce_diff", weights = ~scale_base_weights_f,feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_amce_timeconjoint_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)


#AMCEs rating
amces_subgroup_rating_minu_inclusion_parl_rand <- rating ~ feature1 + feature2 + feature3  +feature4 +feature6
diff_amce_attention_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~attention)
diff_amce_timeinfo_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_amce_timeconjoint_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)

#MMs choice
mm_subgroup <- chosen ~ feature1+ feature2 + feature3 +feature4 +feature6
diff_mm_attention_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~attention)
diff_mm_timeinfo_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_mm_timeconjoint_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)


#MMs rating
mm_subgroup_rating_minu_inclusion_parl_rand <- rating ~ feature1+ feature2 + feature3 +feature4 +feature6
diff_mm_attention_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, mm_subgroup_rating, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~attention)
diff_mm_timeinfo_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, mm_subgroup_rating, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_mm_timeconjoint_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, mm_subgroup_rating, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)




#### ATTENTION

overview_choice_attention_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_attention_choice_minu_inclusion_parl_rand %>% mutate(model = "AMCE (Choice)"),
  diff_mm_attention_choice_minu_inclusion_parl_rand %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_attention_minu_inclusion_parl_rand <- overview_choice_attention_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_attention_minu_inclusion_parl_rand$model <- factor(diff_overview_choice_attention_minu_inclusion_parl_rand$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_attention_plot_minu_inclusion_parl_rand<-plot(diff_overview_choice_attention_minu_inclusion_parl_rand, feature_headers = FALSE, xlab="", xlim = c(-0.15,0.15)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))
overview_choice_attention_plot_minu_inclusion_parl_rand


overview_rating_attention_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_attention_rating_minu_inclusion_parl_rand %>% mutate(model = "AMCE (Rating)"),
  diff_mm_attention_rating_minu_inclusion_parl_rand %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_attention_minu_inclusion_parl_rand <- overview_rating_attention_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_attention_minu_inclusion_parl_rand$model <- factor(diff_overview_rating_attention_minu_inclusion_parl_rand$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_attention_plot_minu_inclusion_parl_rand<-plot(diff_overview_rating_attention_minu_inclusion_parl_rand, feature_headers = FALSE, xlab="", xlim = c(-0.3,0.3)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_attention_plot_minu_inclusion_parl_rand



combined_attention_minu_inclusion_parl_rand <- ggarrange(overview_choice_attention_plot_minu_inclusion_parl_rand + theme(axis.ticks.y = element_blank(),
                                                                                                                         plot.margin = margin(r = 1) ),
                                                         overview_rating_attention_plot_minu_inclusion_parl_rand + theme(
                                                           axis.ticks.y = element_blank(),
                                                           axis.title.y = element_blank(),
                                                           axis.text.y=element_blank(),
                                                           plot.margin = margin(l = 1)), widths = c(1.5,1))

combined_attention_minu_inclusion_parl_rand

annotate_figure(combined_attention_minu_inclusion_parl_rand, bottom = text_grob("Estimated differences for 'high Attention'\ncompared to 'low Attention'", 
                                                                                color = "black", size = 24, hjust= 0.27))


#ggsave("subgroup_overview_att_mm_amce_rating_choice_without_inclusion_parl.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 



#### TIMEINFO

overview_choice_timeinfo_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_timeinfo_choice_minu_inclusion_parl_rand %>% mutate(model = "AMCE (Choice)"),
  diff_mm_timeinfo_choice_minu_inclusion_parl_rand %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_timeinfo_minu_inclusion_parl_rand <- overview_choice_timeinfo_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_timeinfo_minu_inclusion_parl_rand$model <- factor(diff_overview_choice_timeinfo_minu_inclusion_parl_rand$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_timeinfo_plot_minu_inclusion_parl_rand<-plot(diff_overview_choice_timeinfo_minu_inclusion_parl_rand, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_choice_timeinfo_plot_minu_inclusion_parl_rand


overview_rating_timeinfo_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_timeinfo_rating_minu_inclusion_parl_rand %>% mutate(model = "AMCE (Rating)"),
  diff_mm_timeinfo_rating_minu_inclusion_parl_rand %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_timeinfo_minu_inclusion_parl_rand <- overview_rating_timeinfo_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_timeinfo_minu_inclusion_parl_rand$model <- factor(diff_overview_rating_timeinfo_minu_inclusion_parl_rand$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_timeinfo_plot_minu_inclusion_parl_rand<-plot(diff_overview_rating_timeinfo_minu_inclusion_parl_rand, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))
overview_rating_timeinfo_plot_minu_inclusion_parl_rand




combined_timeinfo_minu_inclusion_parl_rand <- ggarrange(overview_choice_timeinfo_plot_minu_inclusion_parl_rand + theme(axis.ticks.y = element_blank(),
                                                                                                                       plot.margin = margin(r = 1) ),
                                                        overview_rating_timeinfo_plot_minu_inclusion_parl_rand + theme(
                                                          axis.ticks.y = element_blank(),
                                                          axis.title.y = element_blank(),
                                                          axis.text.y=element_blank(),
                                                          plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_timeinfo_minu_inclusion_parl_rand

annotate_figure(combined_timeinfo_minu_inclusion_parl_rand, bottom = text_grob("Estimated differences for 'more Time Info'\ncompared to 'less Time Info'", 
                                                                               color = "black", size = 24, hjust= 0.27))


#ggsave("subgroup_overview_timeinfo_mm_amce_rating_choice_without_inclusion_parl.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 



#### TIME CONJOINT

overview_choice_timeconjoint_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_timeconjoint_choice_minu_inclusion_parl_rand %>% mutate(model = "AMCE (Choice)"),
  diff_mm_timeconjoint_choice_minu_inclusion_parl_rand %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_timeconjoint_minu_inclusion_parl_rand <- overview_choice_timeconjoint_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_timeconjoint_minu_inclusion_parl_rand$model <- factor(diff_overview_choice_timeconjoint_minu_inclusion_parl_rand$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_timeconjoint_plot_minu_inclusion_parl_rand<-plot(diff_overview_choice_timeconjoint_minu_inclusion_parl_rand, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))


overview_choice_timeconjoint_plot_minu_inclusion_parl_rand



overview_rating_timeconjoint_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_timeconjoint_rating_minu_inclusion_parl_rand %>% mutate(model = "AMCE (Rating)"),
  diff_mm_timeconjoint_rating_minu_inclusion_parl_rand %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_timeconjoint_minu_inclusion_parl_rand <- overview_rating_timeconjoint_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_timeconjoint_minu_inclusion_parl_rand$model <- factor(diff_overview_rating_timeconjoint_minu_inclusion_parl_rand$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_timeconjoint_plot_minu_inclusion_parl_rand<-plot(diff_overview_rating_timeconjoint_minu_inclusion_parl_rand, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_timeconjoint_plot_minu_inclusion_parl_rand




combined_timeconjoint_minu_inclusion_parl_rand <- ggarrange(overview_choice_timeconjoint_plot_minu_inclusion_parl_rand + theme(axis.ticks.y = element_blank(),
                                                                                                                               plot.margin = margin(r = 1) ),
                                                            overview_rating_timeconjoint_plot_minu_inclusion_parl_rand + theme(
                                                              axis.ticks.y = element_blank(),
                                                              axis.title.y = element_blank(),
                                                              axis.text.y=element_blank(),
                                                              plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_timeconjoint_minu_inclusion_parl_rand

annotate_figure(combined_timeconjoint_minu_inclusion_parl_rand, bottom = text_grob("Estimated differences for 'more Time Conjoint'\ncompared to 'less Time Conjoint'", 
                                                                                   color = "black", size = 24, hjust= 0.27))


#ggsave("subgroup_overview_timeconjoint_mm_amce_rating_choice_without_inclusion_parl.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 

#without parliament
#calculate relevant again

#AMCEs choice
amces_subgroup_choice <- chosen ~ feature1 + feature2 + feature3  +feature4 +feature6
diff_amce_attention_choice_minu_parl <- cj(con1_minu_parl, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~attention)
diff_amce_timeinfo_choice_minu_parl <- cj(con1_minu_parl, amces_subgroup_choice, id = ~ID, estimate="amce_diff", weights = ~scale_base_weights_f,feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_amce_timeconjoint_choice_minu_parl <- cj(con1_minu_parl, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)


#AMCEs rating
amces_subgroup_rating_minu_parl <- rating ~ feature1 + feature2 + feature3  +feature4 +feature6
diff_amce_attention_rating_minu_parl <- cj(con1_minu_parl, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~attention)
diff_amce_timeinfo_rating_minu_parl <- cj(con1_minu_parl, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_amce_timeconjoint_rating_minu_parl <- cj(con1_minu_parl, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)

#MMs choice
mm_subgroup <- chosen ~ feature1+ feature2 + feature3 +feature4 +feature6
diff_mm_attention_choice_minu_parl <- cj(con1_minu_parl, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~attention)
diff_mm_timeinfo_choice_minu_parl <- cj(con1_minu_parl, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_mm_timeconjoint_choice_minu_parl <- cj(con1_minu_parl, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)


#MMs rating
mm_subgroup_rating_minu_parl <- rating ~ feature1+ feature2 + feature3 +feature4 +feature6
diff_mm_attention_rating_minu_parl <- cj(con1_minu_parl, mm_subgroup_rating, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"), by = ~attention)
diff_mm_timeinfo_rating_minu_parl <- cj(con1_minu_parl, mm_subgroup_rating, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),  by = ~info_text_dauer_dicho)
diff_mm_timeconjoint_rating_minu_parl <- cj(con1_minu_parl, mm_subgroup_rating, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names, feature_order = c("feature6", "feature4", "feature3", "feature2", "feature1"),by = ~timing_conjoint_num)




#### ATTENTION

overview_choice_attention_minu_parl <-bind_rows(
  diff_amce_attention_choice_minu_parl %>% mutate(model = "AMCE (Choice)"),
  diff_mm_attention_choice_minu_parl %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_attention_minu_parl <- overview_choice_attention_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_attention_minu_parl$model <- factor(diff_overview_choice_attention_minu_parl$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_attention_plot_minu_parl<-plot(diff_overview_choice_attention_minu_parl, feature_headers = FALSE, xlab="", xlim = c(-0.15,0.15)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))
overview_choice_attention_plot_minu_parl


overview_rating_attention_minu_parl <-bind_rows(
  diff_amce_attention_rating_minu_parl %>% mutate(model = "AMCE (Rating)"),
  diff_mm_attention_rating_minu_parl %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_attention_minu_parl <- overview_rating_attention_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_attention_minu_parl$model <- factor(diff_overview_rating_attention_minu_parl$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_attention_plot_minu_parl<-plot(diff_overview_rating_attention_minu_parl, feature_headers = FALSE, xlab="", xlim = c(-0.3,0.3)) + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_attention_plot_minu_parl



combined_attention_minu_parl <- ggarrange(overview_choice_attention_plot_minu_parl + theme(axis.ticks.y = element_blank(),
                                                                                           plot.margin = margin(r = 1) ),
                                          overview_rating_attention_plot_minu_parl + theme(
                                            axis.ticks.y = element_blank(),
                                            axis.title.y = element_blank(),
                                            axis.text.y=element_blank(),
                                            plot.margin = margin(l = 1)), widths = c(1.5,1))

combined_attention_minu_parl

annotate_figure(combined_attention_minu_parl, bottom = text_grob("Estimated differences for 'high Attention'\ncompared to 'low Attention'", 
                                                                 color = "black", size = 24, hjust= 0.27))


#ggsave("subgroup_overview_att_mm_amce_rating_choice_without_parl.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 



#### TIMEINFO

overview_choice_timeinfo_minu_parl <-bind_rows(
  diff_amce_timeinfo_choice_minu_parl %>% mutate(model = "AMCE (Choice)"),
  diff_mm_timeinfo_choice_minu_parl %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_timeinfo_minu_parl <- overview_choice_timeinfo_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_timeinfo_minu_parl$model <- factor(diff_overview_choice_timeinfo_minu_parl$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_timeinfo_plot_minu_parl<-plot(diff_overview_choice_timeinfo_minu_parl, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_choice_timeinfo_plot_minu_parl


overview_rating_timeinfo_minu_parl <-bind_rows(
  diff_amce_timeinfo_rating_minu_parl %>% mutate(model = "AMCE (Rating)"),
  diff_mm_timeinfo_rating_minu_parl %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_timeinfo_minu_parl <- overview_rating_timeinfo_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_timeinfo_minu_parl$model <- factor(diff_overview_rating_timeinfo_minu_parl$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_timeinfo_plot_minu_parl<-plot(diff_overview_rating_timeinfo_minu_parl, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))
overview_rating_timeinfo_plot_minu_parl




combined_timeinfo_minu_parl <- ggarrange(overview_choice_timeinfo_plot_minu_parl + theme(axis.ticks.y = element_blank(),
                                                                                         plot.margin = margin(r = 1) ),
                                         overview_rating_timeinfo_plot_minu_parl + theme(
                                           axis.ticks.y = element_blank(),
                                           axis.title.y = element_blank(),
                                           axis.text.y=element_blank(),
                                           plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_timeinfo_minu_parl

annotate_figure(combined_timeinfo_minu_parl, bottom = text_grob("Estimated differences for 'more Time Info'\ncompared to 'less Time Info'", 
                                                                color = "black", size = 24, hjust= 0.27))


#ggsave("subgroup_overview_timeinfo_mm_amce_rating_choice_without_parl.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 



#### TIME CONJOINT

overview_choice_timeconjoint_minu_parl <-bind_rows(
  diff_amce_timeconjoint_choice_minu_parl %>% mutate(model = "AMCE (Choice)"),
  diff_mm_timeconjoint_choice_minu_parl %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_timeconjoint_minu_parl <- overview_choice_timeconjoint_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_timeconjoint_minu_parl$model <- factor(diff_overview_choice_timeconjoint_minu_parl$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_timeconjoint_plot_minu_parl<-plot(diff_overview_choice_timeconjoint_minu_parl, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))


overview_choice_timeconjoint_plot_minu_parl



overview_rating_timeconjoint_minu_parl <-bind_rows(
  diff_amce_timeconjoint_rating_minu_parl %>% mutate(model = "AMCE (Rating)"),
  diff_mm_timeconjoint_rating_minu_parl %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_timeconjoint_minu_parl <- overview_rating_timeconjoint_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_timeconjoint_minu_parl$model <- factor(diff_overview_rating_timeconjoint_minu_parl$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_timeconjoint_plot_minu_parl<-plot(diff_overview_rating_timeconjoint_minu_parl, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_timeconjoint_plot_minu_parl




combined_timeconjoint_minu_parl <- ggarrange(overview_choice_timeconjoint_plot_minu_parl + theme(axis.ticks.y = element_blank(),
                                                                                                 plot.margin = margin(r = 1) ),
                                             overview_rating_timeconjoint_plot_minu_parl + theme(
                                               axis.ticks.y = element_blank(),
                                               axis.title.y = element_blank(),
                                               axis.text.y=element_blank(),
                                               plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_timeconjoint_minu_parl

annotate_figure(combined_timeconjoint_minu_parl, bottom = text_grob("Estimated differences for 'more Time Conjoint'\ncompared to 'less Time Conjoint'", 
                                                                    color = "black", size = 24, hjust= 0.27))


#ggsave("subgroup_overview_timeconjoint_mm_amce_rating_choice_without_parl.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 








# Winner vs. Loser --------------------------------------------------------
    # all - no subsets by constraint ------------------------------------------


#calculate relevant again

#AMCEs choice
amces_subgroup_choice <- chosen ~ feature1 + feature2 + feature3  +feature4 
diff_amce_winner_loser_choice <- cj(con1, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature4", "feature3", "feature2", "feature1"),  by = ~Outcome_Favorability)


#AMCEs rating
amces_subgroup_rating <- rating ~ feature1 + feature2 + feature3  +feature4 
diff_amce_winner_loser_rating <- cj(con1, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature4", "feature3", "feature2", "feature1"),  by = ~Outcome_Favorability)

#MMs choice
mm_subgroup <- chosen ~ feature1+ feature2 + feature3 +feature4 
diff_mm_winner_loser_choice <- cj(con1, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature4", "feature3", "feature2", "feature1"), by = ~Outcome_Favorability)


#MMs rating
mm_subgroup_rating <- rating ~ feature1+ feature2 + feature3 +feature4 
diff_mm_winner_loser_rating <- cj(con1, mm_subgroup_rating, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature4", "feature3", "feature2", "feature1"), by = ~Outcome_Favorability)



overview_choice_winnerloser <-bind_rows(
  diff_amce_winner_loser_choice %>% mutate(model = "AMCE (Choice)"),
  diff_mm_winner_loser_choice %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_winnerloser <- overview_choice_winnerloser %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_winnerloser$model <- factor(diff_overview_choice_winnerloser$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_winnerloser_plot<-plot(diff_overview_choice_winnerloser, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_choice_winnerloser_plot


overview_rating_winnerloser <-bind_rows(
  diff_amce_winner_loser_rating %>% mutate(model = "AMCE (Rating)"),
  diff_mm_winner_loser_rating %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_winnerloser <- overview_rating_winnerloser %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_winnerloser$model <- factor(diff_overview_rating_winnerloser$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_winnerloser_plot<-plot(diff_overview_rating_winnerloser, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_winnerloser_plot




combined_winnerloser <- ggarrange(overview_choice_winnerloser_plot + theme(axis.ticks.y = element_blank(),
                                                                           plot.margin = margin(r = 1) ),
                                  overview_rating_winnerloser_plot + theme(
                                    axis.ticks.y = element_blank(),
                                    axis.title.y = element_blank(),
                                    axis.text.y=element_blank(),
                                    plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_winnerloser

annotate_figure(combined_winnerloser, bottom = text_grob("Estimated differences for 'Preference Match'\ncompared to 'Preference Mismatch'", 
                                                         color = "black", size = 14, hjust= 0.16))

ggsave("subgroup_overview_winnerloser_mm_amce_rating_choice.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 



    # subsets by constraint --------------------------------------------------

#without inclusion parliament
#calculate relevant again

#AMCEs choice
amces_subgroup_choice <- chosen ~ feature1 + feature2 + feature3  +feature4 
diff_amce_winner_loser_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature4", "feature3", "feature2", "feature1"),  by = ~Outcome_Favorability)


#AMCEs rating
amces_subgroup_rating <- rating ~ feature1 + feature2 + feature3  +feature4 
diff_amce_winner_loser_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature4", "feature3", "feature2", "feature1"),  by = ~Outcome_Favorability)

#MMs choice
mm_subgroup <- chosen ~ feature1+ feature2 + feature3 +feature4 
diff_mm_winner_loser_choice_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature4", "feature3", "feature2", "feature1"), by = ~Outcome_Favorability)


#MMs rating
mm_subgroup_rating <- rating ~ feature1+ feature2 + feature3 +feature4 
diff_mm_winner_loser_rating_minu_inclusion_parl_rand <- cj(con1_minu_inclusion_parl_rand, mm_subgroup_rating, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature4", "feature3", "feature2", "feature1"), by = ~Outcome_Favorability)



overview_choice_winnerloser_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_winner_loser_choice_minu_inclusion_parl_rand %>% mutate(model = "AMCE (Choice)"),
  diff_mm_winner_loser_choice_minu_inclusion_parl_rand %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_winnerloser_minu_inclusion_parl_rand <- overview_choice_winnerloser_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_winnerloser_minu_inclusion_parl_rand$model <- factor(diff_overview_choice_winnerloser_minu_inclusion_parl_rand$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_winnerloser_plot_minu_inclusion_parl_rand<-plot(diff_overview_choice_winnerloser_minu_inclusion_parl_rand, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))
overview_choice_winnerloser_plot_minu_inclusion_parl_rand


overview_rating_winnerloser_minu_inclusion_parl_rand <-bind_rows(
  diff_amce_winner_loser_rating_minu_inclusion_parl_rand %>% mutate(model = "AMCE (Rating)"),
  diff_mm_winner_loser_rating_minu_inclusion_parl_rand %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_winnerloser_minu_inclusion_parl_rand <- overview_rating_winnerloser_minu_inclusion_parl_rand %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_winnerloser_minu_inclusion_parl_rand$model <- factor(diff_overview_rating_winnerloser_minu_inclusion_parl_rand$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_winnerloser_plot_minu_inclusion_parl_rand<-plot(diff_overview_rating_winnerloser_minu_inclusion_parl_rand, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_winnerloser_plot_minu_inclusion_parl_rand




combined_winnerloser_minu_inclusion_parl_rand <- ggarrange(overview_choice_winnerloser_plot_minu_inclusion_parl_rand + theme(axis.ticks.y = element_blank(),
                                                                                                                             plot.margin = margin(r = 1) ),
                                                           overview_rating_winnerloser_plot_minu_inclusion_parl_rand + theme(
                                                             axis.ticks.y = element_blank(),
                                                             axis.title.y = element_blank(),
                                                             axis.text.y=element_blank(),
                                                             plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_winnerloser_minu_inclusion_parl_rand

annotate_figure(combined_winnerloser_minu_inclusion_parl_rand, bottom = text_grob("Estimated differences for 'Preference Match'\ncompared to 'Preference Mismatch'", 
                                                                                  color = "black", size = 14, hjust= 0.16))

#ggsave("subgroup_overview_winnerloser_mm_amce_rating_choice_without_inclusion_parl.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 


#without parliament
#calculate relevant again

#AMCEs choice
amces_subgroup_choice <- chosen ~ feature1 + feature2 + feature3  +feature4 
diff_amce_winner_loser_choice_minu_parl <- cj(con1_minu_parl, amces_subgroup_choice, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature4", "feature3", "feature2", "feature1"),  by = ~Outcome_Favorability)


#AMCEs rating
amces_subgroup_rating <- rating ~ feature1 + feature2 + feature3  +feature4 
diff_amce_winner_loser_rating_minu_parl <- cj(con1_minu_parl, amces_subgroup_rating, id = ~ID, estimate="amce_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature4", "feature3", "feature2", "feature1"),  by = ~Outcome_Favorability)

#MMs choice
mm_subgroup <- chosen ~ feature1+ feature2 + feature3 +feature4 
diff_mm_winner_loser_choice_minu_parl <- cj(con1_minu_parl, mm_subgroup, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature4", "feature3", "feature2", "feature1"), by = ~Outcome_Favorability)


#MMs rating
mm_subgroup_rating <- rating ~ feature1+ feature2 + feature3 +feature4 
diff_mm_winner_loser_rating_minu_parl <- cj(con1_minu_parl, mm_subgroup_rating, id = ~ID, estimate="mm_diff",weights = ~scale_base_weights_f, feature_labels = feature_names,feature_order = c("feature4", "feature3", "feature2", "feature1"), by = ~Outcome_Favorability)



overview_choice_winnerloser_minu_parl <-bind_rows(
  diff_amce_winner_loser_choice_minu_parl %>% mutate(model = "AMCE (Choice)"),
  diff_mm_winner_loser_choice_minu_parl %>% mutate(model = "Marginal Means (Choice)"))

diff_overview_choice_winnerloser_minu_parl <- overview_choice_winnerloser_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_overview_choice_winnerloser_minu_parl$model <- factor(diff_overview_choice_winnerloser_minu_parl$model, levels = c("AMCE (Choice)",  "Marginal Means (Choice)"))

overview_choice_winnerloser_plot_minu_parl<-plot(diff_overview_choice_winnerloser_minu_parl, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_choice_winnerloser_plot_minu_parl


overview_rating_winnerloser_minu_parl <-bind_rows(
  diff_amce_winner_loser_rating_minu_parl %>% mutate(model = "AMCE (Rating)"),
  diff_mm_winner_loser_rating_minu_parl %>% mutate(model = "Marginal Means (Rating)"))

diff_overview_rating_winnerloser_minu_parl <- overview_rating_winnerloser_minu_parl %>% dplyr::select(BY:level, estimate:model)
diff_overview_rating_winnerloser_minu_parl$model <- factor(diff_overview_rating_winnerloser_minu_parl$model, levels = c("AMCE (Rating)",  "Marginal Means (Rating)"))

overview_rating_winnerloser_plot_minu_parl<-plot(diff_overview_rating_winnerloser_minu_parl, feature_headers = FALSE, xlab="") + 
  ggplot2::facet_wrap(~model, ncol = 2) +
  theme(panel.border = element_rect(linetype = "dotted", colour = "black"),
        panel.background = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 18),
        axis.text.x = element_text(colour = "black", size = 18),
        axis.title = element_text(colour = "black", size = 18),
        panel.grid.major = element_line(colour = "grey", linetype = "dotted"),
        strip.background = element_blank(),
        strip.text.x = element_text(colour = "black", size = 19),
        legend.position = "none")+
  geom_linerange(aes(xmin=lower, xmax=upper), size = 1)+
  geom_point(size=2.5) +
  scale_color_manual(values = rep("midnightblue", 9))

overview_rating_winnerloser_plot_minu_parl




combined_winnerloser_minu_parl <- ggarrange(overview_choice_winnerloser_plot_minu_parl + theme(axis.ticks.y = element_blank(),
                                                                                               plot.margin = margin(r = 1) ),
                                            overview_rating_winnerloser_plot_minu_parl + theme(
                                              axis.ticks.y = element_blank(),
                                              axis.title.y = element_blank(),
                                              axis.text.y=element_blank(),
                                              plot.margin = margin(l = 1)), widths = c(1.48,1))

combined_winnerloser_minu_parl

annotate_figure(combined_winnerloser_minu_parl, bottom = text_grob("Estimated differences for 'Preference Match'\ncompared to 'Preference Mismatch'", 
                                                                   color = "black", size = 14, hjust= 0.16))

#ggsave("subgroup_overview_winnerloser_mm_amce_rating_choice_without_parl.png", width = 28, height=24, units = "cm", dpi=800, scale=1.5, bg="white") 






# #### SUPPLEMENTARY INFORMATION used R version and package versions #### -------------------------------------
# R version 4.3.3 (2024-02-29 ucrt)
# Platform: x86_64-w64-mingw32/x64 (64-bit)
# Running under: Windows 10 x64 (build 19045)
# 
# Matrix products: default
# 
# 
# locale:
#   [1] LC_COLLATE=German_Germany.utf8  LC_CTYPE=German_Germany.utf8    LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C                   
# [5] LC_TIME=German_Germany.utf8    
# 
# time zone: Europe/Berlin
# tzcode source: internal
# 
# attached base packages:
#   [1] grid      stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#   [1] expss_0.11.6      maditr_0.8.4      GDAtools_2.0.1    janitor_2.2.0     openxlsx_4.2.5.2  cregg_0.4.0       corrplot_0.92     DescTools_0.99.52
# [9] sjstats_0.18.2    lme4_1.1-35.1     margins_0.3.26    labelled_2.12.0   knitr_1.45        psych_2.3.9       car_3.1-2         carData_3.0-5    
# [17] ggeffects_1.3.2   haven_2.5.3       lubridate_1.9.3   forcats_1.0.0     stringr_1.5.0     dplyr_1.1.3       purrr_1.0.2       readr_2.1.4      
# [25] tibble_3.2.1      tidyverse_2.0.0   texreg_1.39.3     sjPlot_2.8.15     broom_1.0.5       readxl_1.4.3      survey_4.2-1      survival_3.5-8   
# [33] Matrix_1.6-1.1    stargazer_5.2.3   jtools_2.2.2      ggstance_0.3.6    FactoMineR_2.9    factoextra_1.0.7  MASS_7.3-60.0.1   magrittr_2.0.3   
# [41] sjmisc_2.8.9      lsr_0.5.2         nnet_7.3-19       tidycomm_0.2.1    ggpubr_0.6.0      ggalluvial_0.12.5 qdapTools_1.3.7   tidyr_1.3.0      
# [49] ggplot2_3.4.4     pacman_0.5.1     
# 
# loaded via a namespace (and not attached):
#   [1] splines_4.3.3        bitops_1.0-7         cellranger_1.1.0     XML_3.99-0.15        rpart_4.1.23         lifecycle_1.0.4      rstatix_0.7.2       
# [8] lattice_0.22-5       insight_0.19.6       flashClust_1.01-2    rockchalk_1.8.157    backports_1.4.1      Hmisc_5.1-1          rmarkdown_2.25      
# [15] qgraph_1.9.8         zip_2.3.0            gld_2.6.6            RColorBrewer_1.1-3   pbapply_1.7-2        DBI_1.1.3            minqa_1.2.6         
# [22] abind_1.4-5          expm_0.999-9         quadprog_1.5-8       RCurl_1.98-1.14      sandwich_3.0-2       prediction_0.3.14    ggrepel_0.9.4       
# [29] arm_1.13-1           performance_0.10.8   DT_0.31              tidyselect_1.2.0     farver_2.1.1         matrixStats_1.2.0    stats4_4.3.3        
# [36] base64enc_0.1-3      e1071_1.7-14         Formula_1.2-5        emmeans_1.8.9        tools_4.3.3          chron_2.3-61         Rcpp_1.0.11         
# [43] glue_1.6.2           mnormt_2.1.1         careless_1.2.2       gridExtra_2.3        xfun_0.41            withr_2.5.2          descriptio_1.2      
# [50] fastmap_1.1.1        mitools_2.4          boot_1.3-29          fansi_1.0.5          digest_0.6.33        mi_1.1               timechange_0.2.0    
# [57] R6_2.5.1             estimability_1.4.1   colorspace_2.1-0     gtools_3.9.4         jpeg_0.1-10          utf8_1.2.4           generics_0.1.3      
# [64] data.table_1.14.8    corpcor_1.6.10       class_7.3-22         httr_1.4.7           htmlwidgets_1.6.2    scatterplot3d_0.3-44 pkgconfig_2.0.3     
# [71] Exact_3.2            sem_3.1-15           gtable_0.3.4         lmtest_0.9-40        htmltools_0.5.7      lavaan_0.6-16        semPlot_1.1.6       
# [78] multcompView_0.1-9   scales_1.2.1         lmom_3.0             leaps_3.1            png_0.1-8            snakecase_0.11.1     rstudioapi_0.15.0   
# [85] tzdb_0.4.0           reshape2_1.4.4       coda_0.19-4          checkmate_2.3.0      nlme_3.1-164         nloptr_2.0.3         zoo_1.8-12          
# [92] proxy_0.4-27         sjlabelled_1.2.0     rootSolve_1.8.2.4    parallel_4.3.3       foreign_0.8-86       pillar_1.9.0         vctrs_0.6.4         
# [99] OpenMx_2.21.8        xtable_1.8-4         cluster_2.1.6        htmlTable_2.4.2      evaluate_0.23        pbivnorm_0.6.0       mvtnorm_1.2-3       
# [106] cli_3.6.1            kutils_1.73          compiler_4.3.3       rlang_1.1.2          crayon_1.5.2         ggsignif_0.6.4       labeling_0.4.3      
# [113] modelr_0.1.11        fdrtool_1.2.17       plyr_1.8.9           stringi_1.7.12       pander_0.6.5         munsell_0.5.0        lisrelToR_0.1.5     
# [120] bayestestR_0.13.1    hms_1.1.3            glasso_1.11          igraph_1.5.1         RcppParallel_5.1.7  
# > 