rm(list = ls())
gc()

fpackage.check <- function(packages) {
    lapply(packages, FUN = function(x) {
        if (!require(x, character.only = TRUE)) {
            install.packages(x, dependencies = TRUE)
            library(x, character.only = TRUE)
        }
    })
}

fsave <- function(x, location = "./data/processed/") {
    if (location == "./data/processed/") {
        ifelse(!dir.exists("data"), dir.create("data"), FALSE)
        ifelse(!dir.exists("data/processed"), dir.create("data/processed"), FALSE)
    }

    object = deparse(substitute(x))
    datename <- substr(gsub("[:-]", "", Sys.time()), 1, 8)
    totalname <- paste(location, object, "_", datename, ".rda", sep = "")
    assign(eval(object), x, envir = .GlobalEnv)
    print(paste("SAVED: ", totalname, sep = ""))
    save(list = object, file = totalname)
}

fload <- function(fileName) {
    load(fileName)
    get(ls()[ls() != "fileName"])
}

fshowdf <- function(x, caption = NULL, digits = 2, footnote = NULL, ...) {
    knitr::kable(x, digits = digits, "html", caption = caption, ...) %>%
        kableExtra::footnote(general = footnote) %>%
        kableExtra::kable_styling(bootstrap_options = c("striped", "hover")) %>%
        kableExtra::scroll_box(width = "100%", height = "300px")
}

fcreate_exposure_means <- function(dat, dim_var, comp_var, cont_var, g1, g2, setting, datasource) {
    ii_dim_mean <- dat %>%
        filter(!is.na({
            {
                comp_var
            }
        }) & !is.na({
            {
                cont_var
            }
        })) %>%
        summarize(ii_comp_mean = mean({
            {
                comp_var
            }
        }) * 0.01, ii_cont_mean = mean({
            {
                cont_var
            }
        }) * 0.01)
    ii_dim_g1 <- dat %>%
        filter(!is.na({
            {
                comp_var
            }
        }) & !is.na({
            {
                cont_var
            }
        }) & {
            {
                dim_var
            }
        } == g1) %>%
        summarize(ii_comp_g1 = mean({
            {
                comp_var
            }
        }) * 0.01, ii_cont_g1 = mean({
            {
                cont_var
            }
        }) * 0.01)
    ii_dim_g2 <- dat %>%
        filter(!is.na({
            {
                comp_var
            }
        }) & !is.na({
            {
                cont_var
            }
        }) & {
            {
                dim_var
            }
        } == g2) %>%
        summarize(ii_comp_g2 = mean({
            {
                comp_var
            }
        }) * 0.01, ii_cont_g2 = mean({
            {
                cont_var
            }
        }) * 0.01)
    ii_setting <- paste0(setting)
    ii_dimension <- paste0("Share ", g1)
    ii_datatype <- paste0(datasource)
    ii <- cbind(ii_dimension, ii_setting, ii_datatype, ii_dim_mean, ii_dim_g1, ii_dim_g2)
}

fcreate_exposure_means_robust <- function(dat, dim_var, comp_var, comp_buurt, cont_var, g1, g2, setting) {
    ii_dim_mean <- dat %>%
        filter(!is.na({
            {
                comp_var
            }
        }) & !is.na({
            {
                cont_var
            }
        }) & !is.na({
            {
                comp_buurt
            }
        })) %>%
        summarize(ii_comp_mean = mean({
            {
                comp_var
            }
        }) * 0.01, ii_compbuurt_mean = mean({
            {
                comp_buurt
            }
        }) * 0.01)
    ii_dim_g1 <- dat %>%
        filter(!is.na({
            {
                comp_var
            }
        }) & !is.na({
            {
                cont_var
            }
        }) & !is.na({
            {
                comp_buurt
            }
        }) & {
            {
                dim_var
            }
        } == g1) %>%
        summarize(ii_comp_g1 = mean({
            {
                comp_var
            }
        }) * 0.01, ii_compbuurt_g1 = mean({
            {
                comp_buurt
            }
        }) * 0.01)
    ii_dim_g2 <- dat %>%
        filter(!is.na({
            {
                comp_var
            }
        }) & !is.na({
            {
                cont_var
            }
        }) & !is.na({
            {
                comp_buurt
            }
        }) & {
            {
                dim_var
            }
        } == g2) %>%
        summarize(ii_comp_g2 = mean({
            {
                comp_var
            }
        }) * 0.01, ii_compbuurt_g2 = mean({
            {
                comp_buurt
            }
        }) * 0.01)
    ii_setting <- paste0(setting)
    ii_dimension <- paste0("Share ", g1)
    ii <- cbind(ii_dimension, ii_setting, ii_dim_mean, ii_dim_g1, ii_dim_g2)
}

fcustom_legend <- function(text_entry1, text_entry2) {
    legend <- grobTree(pointsGrob(x = unit(c(0.02, 0.07), "npc"), y = unit(c(0.8, 0.8), "npc"), pch = c(22,
        19), gp = gpar(col = "#56B4E9", fill = "#56B4E9")), textGrob(text_entry1, x = unit(0.11, "npc"),
        y = unit(0.8, "npc"), just = "left", gp = gpar(fontsize = 11)), pointsGrob(x = unit(0.02, "npc"),
        y = unit(0.5, "npc"), pch = 22, gp = gpar(col = "#999999", fill = "#999999")), textGrob("Average intra-setting exposure among\nall respondents",
        x = unit(0.11, "npc"), y = unit(0.5, "npc"), just = "left", gp = gpar(fontsize = 11)), pointsGrob(x = unit(c(0.02,
        0.07), "npc"), y = unit(c(0.2, 0.2), "npc"), pch = c(22, 19), gp = gpar(col = "#E69F00", fill = "#E69F00")),
        textGrob(text_entry2, x = unit(0.11, "npc"), y = unit(0.2, "npc"), just = "left", gp = gpar(fontsize = 11)))
    return(legend)
}

packages = c("tidyverse", "cregg", "ggplot2", "ggtext", "cowplot", "patchwork", "ggpubr", "ggh4x", "gridtext",
    "grid")
fpackage.check(packages)
rm(packages)

#create subgroups
df$Tert <- NA_real_
df$Tert[df$edu == "college degree"] <- 1L
df$Tert[df$edu == "no college degree"] <- 2L
df$Tert <- factor(df$Tert, 1:2, c("With college degree", "Without college degree"))

df$Age <- NA_real_
df$Age[df$age == "50 or older"] <- 1L
df$Age[df$age == "under 50" ] <- 2L
df$Age <- factor(df$Age, 1:2, c("50 years or older", "Below 50 years"))

df$TM <- NA_real_
df$TM[df$eth_dtm == "Dutch"] <- 1L
df$TM[df$eth_dtm == "Turkish/Moroccan"] <- 2L
df$TM <- factor(df$TM, 1:2, c("No migration background", "Turkish/Moroccan background" ))
df$TM <- factor(df$TM, levels = c( "Turkish/Moroccan background", "No migration background"))

#outcome and conjoint features, including only composition features we are interested in;
#and the subgrouping variable

#conditional marginal means (composition dimension; disaggregated by respondent value for this dimension)

#1. nbh
nmmeduc <- cregg::cj(df[!is.na(df$Tert), ], enbh_choice ~ enbh_comp_educ, id = ~id, estimate = "mm", by = ~Tert)
nmmage <- cregg::cj(df[!is.na(df$Age), ], enbh_choice ~ enbh_comp_age, id = ~id, estimate = "mm", by = ~Age)
nmmeth <- cregg::cj(df[!is.na(df$TM), ], enbh_choice ~ enbh_comp_eth, id = ~id, estimate = "mm", by = ~TM)

#2. org
cmmeduc <- cregg::cj(df[!is.na(df$Tert), ], eciv_choice ~ eciv_comp_educ, id = ~id, estimate = "mm", by = ~Tert)
cmmage <- cregg::cj(df[!is.na(df$Age), ], eciv_choice ~ eciv_comp_age, id = ~id, estimate = "mm", by = ~Age)
cmmeth <- cregg::cj(df[!is.na(df$TM), ], eciv_choice ~ eciv_comp_eth, id = ~id, estimate = "mm", by = ~TM)

#rowbind the dataframes
#leave out last column, reflecting the subgrouping variable
mms <- rbind(nmmage[-12],nmmeth[-12],nmmeduc[-12],
             cmmage[-12],cmmeth[-12],cmmeduc[-12])

#rename labels
#1. of outcome (setting)
mms$Outcome <- factor(ifelse(mms$outcome == "eciv_choice", "Civic organizations", "Neighborhoods"),
                       levels = c("Neighborhoods", "Civic organizations"))

#2. of (demographic composition) features
mms$Feature <- factor(sub(".*_", "", mms$feature), 
                        levels = c("age", "eth", "educ"),
                        labels = c("Share of people aged 50 years or older", "Share of people with Turk/Moroc background",  "Share of people with a college degree"))

mms$Subgroup <- mms$BY
mms$Subgroup <- factor(mms$Subgroup, levels = c("50 years or older", "Below 50 years", "Turkish/Moroccan background",   "No migration background" , "With college degree","Without college degree"  )) 

# nice color palette
cbPalette <- c("#999999","#E69F00","#56B4E9","#D55E00","#F0E442","#0072B2")

p2 <- plot(mms, vline=1/2, feature_headers = NA, size = 2, group = "Subgroup", xlab = "Marginal mean") +
  facet_grid2( Outcome ~ Feature, scales = "free", independent = "y", switch = "y") +
  geom_text( aes(label = sprintf("%0.2f (%0.2f)", estimate, std.error)), size = 3, colour = "black", vjust = -1, position = position_dodge(width = 0.9)) +
  scale_color_manual(values = cbPalette) +
  facetted_pos_scales(
    y = list(
      Feature == "Share of people with a college degree" ~ scale_y_discrete(
        labels = c("three quarters (75%)", "half (50%)", "a quarter (25%)")),
      Feature == "Share of people aged 50 years or older" ~ scale_y_discrete(
        labels = c("half (50%)", "just under half (40%)", "a quarter (25%)")),
      Feature == "Share of people with Turk/Moroc background" ~ scale_y_discrete(
        labels = c("a quarter (25%)", "minority (10%)", "none (0%)"))
      ),
    #freeing up the x scales does not work for some reason... so i set them manually:
    x = list(
      Feature == "Share of people with a college degree" ~ scale_x_continuous(
        limits = c(.40,.60), labels = scales::number_format(accuracy = 0.01)),
      Feature == "Share of people aged 50 years or older" ~ scale_x_continuous(
        limits = c(.45,.60), labels = scales::number_format(accuracy = 0.01)),
      Feature == "Share of people with Turk/Moroc background" ~ scale_x_continuous(
        limits = c(.30,.70), labels = scales::number_format(accuracy = 0.01))
      )) +
  theme(legend.position = "none",
       # strip.text.x = element_markdown(),
        strip.placement = "outside")

#split legend into multiple parts
split1 <- mms[mms$Subgroup %in% c("50 years or older", "Below 50 years"), ]
split2 <- mms[mms$Subgroup %in% c("Turkish/Moroccan background", "No migration background"), ]
split3 <- mms[mms$Subgroup %in% c("With college degree", "Without college degree"), ]

#reverse order factor.... bug in ggplot2 makes the order of the subgroups in the plot and the legend reversed
split1$Subgroup <- factor(split1$Subgroup, levels = rev(levels(split1$Subgroup)))
split2$Subgroup <- factor(split2$Subgroup, levels = rev(levels(split2$Subgroup)))
split3$Subgroup <- factor(split3$Subgroup, levels = rev(levels(split3$Subgroup)))

#plot legends
split1p <- ggplot(split1, aes(x = estimate, y = level, color=Subgroup)) + geom_point() + labs(color = "Respondent subgroups:") + theme_bw() + scale_color_manual(values = cbPalette[c(2,1)]) 
split2p <- ggplot(split2, aes(x = estimate, y = level, color=Subgroup)) + geom_point() + labs(color = "Respondent subgroups:") + theme_bw() + scale_color_manual(values = cbPalette[c(4,3)])
split3p <- ggplot(split3, aes(x = estimate, y = level, color=Subgroup)) + geom_point() + labs(color = "Respondent subgroups:") + theme_bw() + scale_color_manual(values = cbPalette[c(6,5)])

#get legends
leg1 <- get_legend(split1p)
leg2 <- get_legend(split2p)
leg3 <- get_legend(split3p)

#combine legends
leg123 <- plot_grid(leg1,leg2,leg3, nrow=1, align="h")

#create a blank plot for legend alignment
blank_p <- plot_spacer()

#combine all in one
p3 <- plot_grid(p2, leg123, align = "h", nrow=2, rel_heights = c(7,1))
plot(p3)

#ggsave("./figures/cmms_liss.png", p3, width=13) 

#test of preference heterogeneity (nested model comparison test
cregg::cj_anova(df[!is.na(df$Tert), ], eciv_choice ~ eciv_comp_educ, by = ~Tert)

cregg::cj_anova(df[!is.na(df$Tert), ], eciv_choice ~ eciv_comp_eth, by = ~Tert)

cregg::cj_anova(df[!is.na(df$Tert), ], eciv_choice ~ eciv_comp_age, by = ~Tert)

# function to calculate (subgroup) IRR
f_irr <- function(df) {
    irr_cd <- 1 - mean(df$repchoice[df$Tert == "With college degree"], na.rm = TRUE)  #college degree
    irr_ncd <- 1 - mean(df$repchoice[df$Tert == "Without college degree"], na.rm = TRUE)  #no college degree
    irr_50 <- 1 - mean(df$repchoice[df$Age == "50 years or older"], na.rm = TRUE)  #50 plus
    irr_n50 <- 1 - mean(df$repchoice[df$Age == "Below 50 years"], na.rm = TRUE)  #under 50
    irr_d <- 1 - mean(df$repchoice[df$TM == "No migration background"], na.rm = TRUE)  #dutch
    irr_tm <- 1 - mean(df$repchoice[df$TM == "Turkish/Moroccan background"], na.rm = TRUE)  #turk/moroc

    return(list(`With college degree` = irr_cd, `Without college degree` = irr_ncd, `50 years or older` = irr_50,
        `Below 50 years` = irr_n50, `No migration background` = irr_d, `Turkish/Moroccan background` = irr_tm))
}

# function to calculate swapping error
f_tau <- function(irr_list) {
    sapply(irr_list, function(irr) (1 - sqrt(1 - 2 * (1 - irr)))/2)
}

# define data-sets and (repeated) choice tasks
datasets <- list(df[!is.na(df$rep_org_choice), ], df[!is.na(df$rep_nbh_choice), ])
tasks <- list(c("first_org_choice", "rep_org_choice"), c("first_nbh_choice", "rep_nbh_choice"))
labels <- c("organizations", "neighborhoods")

# use a loop to store swapping errors per subgroup, disaggregated by setting (org/nbh)
swapping <- list()
for (i in 1:length(datasets)) {
    # select dataset
    dataset <- datasets[[i]]

    # check whether repeated choicetasks were answered the same by the respondent
    dataset$repchoice <- ifelse(dataset[[tasks[[i]][1]]] == dataset[[tasks[[i]][2]]], 1, 0)

    # calculate IRR (note that the options/alternatives were reversed! that's why we subtract the
    # IRR from 1)
    irr_list <- f_irr(dataset)

    # calculate tau
    tau_list <- f_tau(irr_list)

    # store results in the list
    swapping[[labels[i]]] <- list(irr = irr_list, tau = tau_list)
}

# store this, as we need it later:
swappingList <- swapping

print(swapping)

# necessary packages
library(parallel)
library(boot)

#set.seed(52235235) #seed for replication

fboot <- function(data, i) { 
  # sample respondent IDs with replacement
  ids <- sample(unique(data$id), length(unique(data$id)), replace = TRUE)
  
  # create a new dataset by iterating over each sampled id and binding the matching rows
  d <- map_dfr(ids, ~ data %>% filter(id == .x)) 

  # make two datasets, for both experiments
  datasets <- list(d[!is.na(d$rep_org_choice), ], d[!is.na(d$rep_nbh_choice), ])
  
  # get the repeated tasks in both experiments
  tasks <- list(c("first_org_choice", "rep_org_choice"), c("first_nbh_choice", "rep_nbh_choice"))
  
  # label the experiments
  labels <- c("organizations", "neighborhoods")
  
  swapping <- list()
  for (i in 1:length(datasets)) { # for experiment i
    # get the dataset 
    dataset <- datasets[[i]]
    
    # identify whether in the repeated choice task, a different (1) or similar (0) choice was made
    dataset$repchoice <- ifelse(dataset[[tasks[[i]][1]]] == dataset[[tasks[[i]][2]]], 1, 0)

    #calculate IRR and swapping error
    irr_list <- f_irr(dataset)
    tau_list <- f_tau(irr_list)
  
    swapping[[labels[i]]] <- list("irr" = irr_list, "tau" = tau_list)
    }
     
  # calculate conditional marginal means:
  x_vars <- cbind(c("Age", "TM", "Tert"), c("age", "eth", "educ"))
  y_vars <- c("enbh", "eciv")
    
  results <- data.frame()
  
  for(x in 1:nrow(x_vars)) {
    for(y in y_vars) {
      # formula  
      f <- as.formula( paste0(y, "_choice ~ ", y, "_comp_", x_vars[x,2]))
      by <- as.formula(paste0("~", x_vars[x,1]))
      cmms <- cregg::cj( d[!is.na(d[x_vars[x,1]]), ], #exclude missings on the subgrouping variable
                           formula = f, id = ~id, estimate = "mm", by = by)
      #bind results
      results <- rbind(results, cmms[,-12])
      }
  }
 
  #calculate corrected conditional marginal mean (for subgroups)
  results$estimate_corr <- NA
    
  for (i in 1:nrow(results)) { 
    #retrieve setting and subgroup
    setting <- ifelse(results$outcome[i] == "enbh_choice", "neighborhoods", "organizations")
    subgroup <- as.character(results$BY[i])

    tau <- swapping[[setting]]$tau[subgroup]
    
    #define bias corrected MM (eq. 11 of Clayton et al.)
    results$estimate_corr[i] <- (results$estimate[i] - tau) / (1 - 2*tau)
  }
   
  #save rownames
  #rnames <- paste0( results$outcome, "-", results$level, "-", results$BY)
  #save(rnames, file = "./data/rownames.Rda")
  
  return(results$estimate_corr)  
}
#fboot(df)#test

nCore <- detectCores() - 1 #detect number of cores (keep one for other tasks)
nIter <- 1000 #number of bootstrap iterations

#set up cluster
mycl <- makeCluster(rep("localhost", nCore)) 

#load necessary library on each cluster node
clusterEvalQ(mycl, {
  library(cregg)
  library(purrr)
  library(dplyr)
})

#copy data and custom functions to each node's workspace
clusterExport(mycl, varlist=c("df", "f_tau", "f_irr", "fboot")) 

#perform bootstrap (takes very long!!!)
system.time(test <- boot(df, fboot, R=nIter, parallel = "snow", ncpus = nCore, cl = mycl)) #about 2 hours on my pc

#save(test, file = "./results/boot.Rda")
stopCluster(mycl)

load("./results/boot.Rda")
load("./data/rownames.Rda")
nIter = 1000

# Some NA values may occur because, in certain bootstrap iterations, no individuals from minority groups are sampled due to their small sample size. replace these NA values with the mean of the corresponding bootstrapped CMM estimates, ignoring any NAs during the calculation.
for (i in 1:length(test$t0)) {
  if(is.na(test$t0[i])) {
    test$t0[i] <- mean(test$t[,i][is.finite(test$t[, i])], na.rm=TRUE)
  }
}

test1 <- summary(test)

#add rownames
rownames(test1) <- rnames
#str(test)
 
#first add the original bias corrected CMMs to the dataframe (instead of the bootstrap estimates):
#create a merge variable in the original CMM dataframe
mms$BY2 <- paste0(mms$outcome, "-", mms$level, "-", mms$BY)
#reorder
mms <- mms[match(rownames(test1), mms$BY2),]

#now calculate the corrected CMMS
for (i in 1:nrow(mms)) { # for each CMM
    #retrieve setting and subgroup
    setting <- ifelse(mms$outcome[i] == "enbh_choice", "neighborhoods", "organizations")
    subgroup <- as.character(mms$BY[i])
    
    #get tau from our saved list
    tau <- swappingList[[setting]]$tau[subgroup]
    #define bias corrected MM (eq. 11 of Clayton et al.)
    mms$estimate_corr[i] <- (mms$estimate[i] - tau) / (1 - 2*tau)
}

#cor(cbind(mms$estimate_corr, test1$original)) # bias is negligible

# i also calculate the difference in the bootstrapped CMM estimates for each composition attribute level, between (1) subgroups (e.g., college/no college degree); and (2) settigns (i.e., nbh vs org)

# create an empty list to store comparison pairs
comparisons <- list()

# loop through each element in rnames
for (i in rnames) {
  # extract relevant components
  components <- unlist(strsplit(i, "-"))
  demo_trait <- components[2] #e.g., age: half (50%)
  choice_type <- components[1] # e.g., enbh_choice-age
  group <- components[3] # e.g., 50 years or older

  # compare with other groups
  for (j in rnames) {
    if (i != j) { 
      other_components <- unlist(strsplit(j, "-"))
      other_choice_type <- other_components[1]
      other_demo_trait <- other_components[2]
      other_group <- other_components[3]
      # compare subgroups (within setting)
      if( group != other_group & demo_trait == other_demo_trait & choice_type == other_choice_type ) {
        comparisons[[length(comparisons) + 1]] <- data.frame(
          focal = i, other = j, comparison = "subgroup", stringsAsFactors = FALSE
        )
      }
      # compare settings (within subgroup)
        if( group == other_group & demo_trait == other_demo_trait & choice_type != other_choice_type ) {
        comparisons[[length(comparisons) + 1]] <- data.frame(
          focal = i, other = j, comparison = "setting", stringsAsFactors = FALSE
        )
      }
    }
  }
}
      
# convert list into dataframe
comparison_df <- do.call(rbind, comparisons)

#now, for each comparison-pair, get the bootstrapped (error-adjusted) CMMs, and calculate the mean and CI of the difference to assess whether the difference is significant (indicating interaction)
comparison_df$mean_dif <- NA
comparison_df$CI <- NA
comparison_df$sig <- NA

for (i in 1:nrow(comparison_df)) {
  #get bootstrapped CMMs 
  focal = comparison_df$focal[i] 
  other = comparison_df$other[i]
  
  #calculate difference between bootstrap estimates of groups (i.e., subgroups/settings)
  cmm_difs <- test$t[, which(rnames == focal)] - test$t[, which(rnames == other)]
  
  #exclude missings and infinite values
  cmm_difs <- cmm_difs[!is.na(cmm_difs) & is.finite(cmm_difs)]
  
  #get the mean difference
  mean_cmm_dif <- mean(cmm_difs)
  
  #standard error
  se_cmm_dif <- sd(cmm_difs) / sqrt(length(cmm_difs))
  
  #critical t value
  t_value <- qt(0.975, df = length(cmm_difs) - 1)
  
  #margin of error
  margin_error <- se_cmm_dif * t_value
  
  #95% CI
  lower_bound <- mean_cmm_dif - margin_error
  upper_bound <- mean_cmm_dif + margin_error
  
  #add to df
  comparison_df$mean_dif[i] <- mean_cmm_dif
  comparison_df$CI[i] <- paste0("[",round(lower_bound,3), ", ", round(upper_bound,3),"]")
  comparison_df$sig[i] <- ifelse((lower_bound <= 0) & (upper_bound >= 0), FALSE, TRUE)
  }

#calculate standard deviation of bootstrap statistics
for (i in 1:nrow(test1)) {
  if(test1$R[i] < nIter) {
    stat <- test$t[,i]
    test1$bootSE[i] <- sd(stat[!is.na(stat) & is.finite(stat)])
  }
}

test1$dif_subgr_CI <- test1$dif_sett_CI <- NA

for (i in 1:nrow(test1)) { 
  focal = rownames(test1)[i]
  test1$dif_subgr_CI[i] <- comparison_df$CI[comparison_df$focal == focal & comparison_df$comparison == "subgroup"]
  test1$dif_sett_CI[i] <- comparison_df$CI[comparison_df$focal == focal & comparison_df$comparison == "setting"]
}

fshowdf(test1, caption = "bootstrap results", footnote = "For each bootstrap iteration, we calculated the difference in estimated conditional marginal means for composition attribute levels (e.g., 25% college educated) *across* subgroups (e.g., with vs. without college degree) within the same setting (e.g., neighborhood). The 95% confidence interval (`dif_subgr_CI`) indicates whether the subgroup differences were significant. Additionally, we tested whether preferences for composition attribute levles *within* subgroups differed significantly between settings (`dif_sett_CI`).")

#add these statistics to plotting data
#create copy of mms
mms_corr <- mms
#use corrected estimates
mms_corr$estimate <- mms_corr$estimate_corr

#with bootstrapped standard errors
mms_corr$std.error <- test1$bootSE
#and confidence intervals
mms_corr$upper <- mms_corr$estimate + 1.96*mms_corr$std.error
mms_corr$lower <- mms_corr$estimate - 1.96*mms_corr$std.error

mms_corr$estimate[mms_corr$estimate<0] <- 0

p <- plot(mms_corr, vline=1/2, feature_headers = NA, size = 2, group = "Subgroup", xlab = "Marginal mean") +
  facet_grid2( Outcome ~ Feature, scales = "free", independent = "y", switch = "y") +
  geom_text( aes(label = sprintf("%0.2f (%0.2f)", estimate, std.error)), size = 3, colour = "black", vjust = -1, position = position_dodge(width = 0.9)) +
  scale_color_manual(values = cbPalette) +
  facetted_pos_scales(
    y = list(
      Feature == "Share of people with a college degree" ~ scale_y_discrete(
        labels = c("three quarters (75%)", "half (50%)", "a quarter (25%)")),
      Feature == "Share of people aged 50 years or older" ~ scale_y_discrete(
        labels = c("half (50%)", "just under half (40%)", "a quarter (25%)")),
      Feature == "Share of people with Turk/Moroc background" ~ scale_y_discrete(
        labels = c("a quarter (25%)", "minority (10%)", "none (0%)"))
      ),
    #freeing up the x scales does not work for some reason... so i set them manually:
    x = list(
      Feature == "Share of people with a college degree" ~ scale_x_continuous(
        limits = c(.35,.65), labels = scales::number_format(accuracy = 0.01)),
      Feature == "Share of people aged 50 years or older" ~ scale_x_continuous(
        limits = c(.40,.60), labels = scales::number_format(accuracy = 0.01)),
      Feature == "Share of people with Turk/Moroc background" ~ scale_x_continuous(
        limits = c(-0.09,1), labels = scales::number_format(accuracy = 0.01))
      )) +
  theme(legend.position = "none",
       # strip.text.x = element_markdown(),
        strip.placement = "outside")
 
p <- plot_grid(p, leg123, align = "h", nrow=2, rel_heights = c(7,1))

fshowdf(mms_corr[,-c(14, 15,16)], caption = "Measurement error bias corrected conditional marginal means", digits=3)

plot(p)

#ggsave("./figures/cmms_liss_corrected.png", p, width=13) 
#ggsave("./figures/figure02_subgroup_mm.png", p, width=13)

### Table for supplement
mms_corr2 <- mms_corr[, c("Subgroup", "Feature", "Outcome", "level", "estimate" )]
mms_corr2$level <- sub(".*: ", "", mms_corr2$level)

#to wide:
mms_corr_wide <- mms_corr2 %>%
  pivot_wider(
    names_from = Outcome,
    values_from = estimate
  )
names(mms_corr_wide)[4:5] <- c("nbh_estimate", "civ_estimate")

# calculate difference
mms_corr_wide$dif_estimate <- mms_corr_wide$nbh_estimate -  mms_corr_wide$civ_estimate
# calculate ratio
mms_corr_wide$ratio <- (mms_corr_wide$nbh_estimate - 0.5) / (mms_corr_wide$civ_estimate - 0.5)

#to 3 digits
mms_corr_wide_d <- mms_corr_wide %>%
  mutate(across(4:7, round, digits = 3))

#write.csv(mms_corr_wide_d, file = "tab.csv", row.names = FALSE)
#write.csv(comparison_df, file = "comp.csv", row.names=FALSE)

#create subgroups

df2$Tert <- NA_real_
df2$Tert[df2$tert_educ == 1] <- 1L
df2$Tert[df2$tert_educ == 0] <- 2L
df2$Tert <- factor(df2$Tert, 1:2, c("With college degree", "Without college degree"))

df2$Migr <- NA_real_
df2$Migr[df2$migration == 1] <- 1L
df2$Migr[df2$migration == 0] <- 2L
df2$Migr <- factor(df2$Migr, 1:2, c( "Migration background","No migration background"))

#outcome and conjoint features, including only composition features we are interested in;
#and the subgrouping variable

#conditional marginal means (composition dimension; disaggregated by respondent value for this dimension)

#1. club
cmmeduc <- cregg::cj(df2[!is.na(df2$Tert), ], choice ~ educ_club, id = ~id, estimate = "mm", by = ~Tert)
cmmeth <- cregg::cj(df2[!is.na(df2$Migr), ], choice ~ migrant_club, id ~id, estimate = "mm", by = ~Migr)

#2. team
tmmeduc <- cregg::cj(df2[!is.na(df2$Tert), ], choice ~ educ_team, id = ~id, estimate = "mm", by = ~Tert)
tmmeth <- cregg::cj(df2[!is.na(df2$Migr), ], choice ~ migrant_team, id ~id, estimate = "mm", by = ~Migr)

#rowbind the dataframes
#leave out last column, reflecting the subgrouping variable
mms <- rbind(cmmeduc[-12],tmmeduc[-12],
             cmmeth[-12],tmmeth[-12])

#rename labels
#1. of outcome (setting)
mms$Outcome <- factor(ifelse( grepl("club", mms$feature), "Club", "Team / training group"), levels = c("Club", "Team / training group"))

#2. of (demographic composition) features
mms$Feature <- factor(sub("_.+", "", mms$feature), levels = c("migrant", "educ"),
                      labels = c("Share of people with a migration background", "Share of people with a college degree"))

mms$Subgroup <- mms$BY

# nice color palette
cbPalette <- c("#F0E442","#0072B2", "#56B4E9","#D55E00")

p <- plot(mms, vline=1/2, feature_headers = NA, size = 2, group = "Subgroup") +
  facet_grid2( Outcome ~ Feature, scales = "free", independent = "y", switch = "y") +
  geom_text( aes(label = sprintf("%0.2f (%0.2f)", estimate, std.error)), size = 3, colour = "black", vjust = -1, position = position_dodge(width = 0.9)) +
  scale_color_manual(values = cbPalette) +
  facetted_pos_scales(
    y = list(
      Feature == "Share of people with a college degree" & Outcome == "Club" ~ scale_y_discrete(
        labels = c("three quarters (75%)", "half (50%)", "a quarter (25%)")),
      
      Feature == "Share of people with a college degree"  & Outcome == "Team / training group" ~ scale_y_discrete(
        labels = c("majority (90%)", "three quarters (75%)", "half (50%)", "a quarter (25%)", "minority (10%)")),
      
      Feature == "Share of people with a migration background" & Outcome == "Club" ~ scale_y_discrete(
        labels = c("a quarter (25%)", "minority (10%)", "a few (5%)")),
      
      Feature == "Share of people with a migration background" & Outcome == "Team / training group" ~ scale_y_discrete(
        labels = c("half (50%)", "a quarter (25%)", "minority (10%)", "a few (5%)", "none (0%)"))
    ),
    #freeing up the x scales does not work for some reason... so i set them manually:
    x = list(
      Feature == "Share of people with a college degree" ~ scale_x_continuous(
        limits = c(.42,.58), labels = scales::number_format(accuracy = 0.01)),
      Feature == "Share of people with a migration background" ~ scale_x_continuous(
        limits = c(.42,.58), labels = scales::number_format(accuracy = 0.01))
    )
  ) +
  theme(legend.position = "none",
        #strip.text.x = element_markdown(),
        strip.placement = "outside")

#split legend into multiple parts
split1 <- mms[mms$Subgroup %in% c("Migration background", "No migration background"), ]
split2 <- mms[mms$Subgroup %in% c("With college degree", "Without college degree"), ]

#reverse order factor.... bug in ggplot2 makes the order of the subgroups in the plot and the legend reversed
split1$Subgroup <- factor(split1$Subgroup, levels = rev(levels(split1$Subgroup)))
split2$Subgroup <- factor(split2$Subgroup, levels = rev(levels(split2$Subgroup)))

#plot legends
split1p <- ggplot(split1, aes(x = estimate, y = level, color=Subgroup)) + geom_point() + labs(color = "Respondent subgroups:") + theme_bw() + scale_color_manual(values = cbPalette[c(4,3)]) 
split2p <- ggplot(split2, aes(x = estimate, y = level, color=Subgroup)) + geom_point() + labs(color = "Respondent subgroups:") + theme_bw() + scale_color_manual(values = cbPalette[c(2,1)])

#get legends
leg1 <- get_legend(split1p)
leg2 <- get_legend(split2p)

#combine legends
leg12 <- plot_grid(leg1,leg2, nrow=1, align="h")

#create a blank plot for legend alignment
blank_p <- plot_spacer()

#combine all in one
p <- plot_grid(p, leg12, align = "h", nrow=2, rel_heights = c(9,1))

print(p)

#ggsave("./figures/cmms_trial.png", p, width=9) 

# create dataframe
swapping2 <- data.frame(x = c("With college degree", "Without college degree", "Migration background",
    "No migration background"), tau = NA)

# we retrieve LISS-based swapping errors
swapping2$tau <- ifelse(swapping2$x == "With college degree", swapping$organizations$tau[["With college degree"]],
    ifelse(swapping2$x == "Without college degree", swapping$organizations$tau[["Without college degree"]],
        ifelse(swapping2$x == "No migration background", swapping$organizations$tau[["No migration background"]],
            ifelse(swapping2$x == "Migration background", swapping$organizations$tau[["Turkish/Moroccan background"]],
                NA))))

fboot <- function(data, i) {
    # sample respondent IDs with replacement
    ids <- sample(unique(data$id), length(unique(data$id)), replace = TRUE)

    # create a new dataset by iterating over each sampled id and binding the matching rows
    d <- map_dfr(ids, ~data %>%
        filter(id == .x))

    # calculate conditional marginal means:
    x_vars <- cbind(c("Migr", "Tert"), c("migrant", "educ"))
    y_vars <- c("club", "team")

    results <- data.frame()

    for (x in 1:nrow(x_vars)) {
        for (y in y_vars) {

            # formula
            f <- as.formula(paste0("choice ~ ", x_vars[x, 2], "_", y))
            by <- as.formula(paste0("~", x_vars[x, 1]))

            cmms <- cregg::cj(d[!is.na(d[x_vars[x, 1]]), ], formula = f, id = ~id, estimate = "mm", by = by)

            # bind results
            results <- rbind(results, cmms[, -12])
        }
    }

    # calculate corrected conditional marginal mean (for subgroups)
    results$estimate_corr <- NA

    for (i in 1:nrow(results)) {
        # retrieve subgroup
        subgroup <- as.character(results$BY[i])
        tau <- swapping2$tau[swapping2$x == subgroup]

        # define bias corrected MM (eq. 11 of Clayton et al.)
        results$estimate_corr[i] <- (results$estimate[i] - tau)/(1 - 2 * tau)

    }
    # save rownames rnames3 <- paste0( results$outcome, '-', results$level, '-', results$BY)
    # save(rnames3, file = './data/rownames3.Rda')

    return(results$estimate_corr)
}
# fboot(df2)#test

nCore <- detectCores() - 1  #detect number of cores (keep one for other tasks)
nIter <- 1000  #number of bootstrap iterations

# set up cluster
mycl <- makeCluster(rep("localhost", nCore))

# load necessary library on each cluster node
clusterEvalQ(mycl, {
    library(cregg)
    library(purrr)
    library(dplyr)
})

# copy data and custom functions to each node's workspace
clusterExport(mycl, varlist = c("df2", "f_tau", "f_irr", "fboot", "swapping2"))

# perform bootstrap:
system.time(test <- boot(df2, fboot, R = nIter, parallel = "snow", ncpus = nCore, cl = mycl))  #about 30 min on my pc
# save(test, file = './results/boot5.Rda')
stopCluster(mycl)

load("./results/boot5.Rda")
load("./data/rownames3.Rda")
nIter = 1000

test1 <- summary(test)
#add rownames
rownames(test1) <- rnames3

#print(test1)

#first add the original bias corrected CMMs to the dataframe:
#create a merge variable in the original CMM dataframe
mms$BY2 <- paste0(mms$outcome, "-", mms$level, "-", mms$BY)
#reorder
mms <- mms[match(rownames(test1), mms$BY2),]

#now calculate the corrected CMMS
for (i in 1:nrow(mms)) { # for each CMM
    #retrieve setting and subgroup
    setting <- "organizations" #use LISS-based swapping error for organizations experiment
    subgroup <- as.character(mms$BY[i])
    subgroup <- ifelse(subgroup == "Migration background", "Turkish/Moroccan background", subgroup) # if migration background, use swapping error of TM
    #get tau from our saved list
    tau <- swappingList[[setting]]$tau[subgroup]
    #define bias corrected MM (eq. 11 of Clayton et al.)
    mms$estimate_corr[i] <- (mms$estimate[i] - tau) / (1 - 2*tau)
}

#add these statistics to plotting data

#create copy of mms
mms_corr <- mms

#use corrected estimates
mms_corr$estimate <- mms_corr$estimate_corr

#with bootstrapped standard errors
mms_corr$std.error <- test1$bootSE
#and confidence intervals
mms_corr$upper <- mms_corr$estimate + 1.96*mms_corr$std.error
mms_corr$lower <- mms_corr$estimate - 1.96*mms_corr$std.error

#plot
p <- plot(mms_corr, vline=1/2, feature_headers = NA, size = 2, group = "Subgroup", xlab = "Marginal mean") +
  facet_grid2( Outcome ~ Feature, scales = "free", independent = "y", switch = "y") +
  geom_text( aes(label = sprintf("%0.2f (%0.2f)", estimate, std.error)), size = 3, colour = "black", vjust = -1, position = position_dodge(width = 0.9)) +
  scale_color_manual(values = cbPalette) +
  facetted_pos_scales(
    y = list(
      Feature == "Share of people with a college degree" & Outcome == "Club" ~ scale_y_discrete(
        labels = c("three quarters (75%)", "half (50%)", "a quarter (25%)")),
      
      Feature == "Share of people with a college degree"  & Outcome == "Team / training group" ~ scale_y_discrete(
        labels = c("majority (90%)", "three quarters (75%)", "half (50%)", "a quarter (25%)", "minority (10%)")),
        
      Feature == "Share of people with a migration background" & Outcome == "Club" ~ scale_y_discrete(
        labels = c("a quarter (25%)", "minority (10%)", "a few (5%)")),
      
      Feature == "Share of people with a migration background" & Outcome == "Team / training group" ~ scale_y_discrete(
        labels = c("half (50%)", "a quarter (25%)", "minority (10%)", "a few (5%)", "none (0%)"))
      ),
    #freeing up the x scales does not work for some reason... so i set them manually:
    x = list(
      Feature == "Share of people with a college degree" ~ scale_x_continuous(
        limits = c(.38,.60), labels = scales::number_format(accuracy = 0.01)),
      Feature == "Share of people with a migration background" ~ scale_x_continuous(
        limits = c(.28,.70), labels = scales::number_format(accuracy = 0.01))
    )
    ) +
  theme(legend.position = "none",
        #strip.text.x = element_markdown(),
        strip.placement = "outside")

p <- plot_grid(p, leg12, align = "h", nrow=2, rel_heights = c(9,1))
fshowdf(mms_corr[,-c(15,16)], caption = "Measurement error bias corrected conditional marginal means")