## Simulation code for Ache & Hütter (submitted). A consistency account of advice taking: Integrating the effects of advice consensus and advice distance

# Packages ----------------------------------------------------------------
rm(list = ls())

library(distr)
library(tidyverse)

# Simulation in main text (mixture distribution) --------------------------

# Create empty data frame to store results
dat <- NULL

# Set parameters for prior knowledge
meanA <- 0
sdA <- 1

# Set seed to make reproducible
set.seed(20201013)

# Loop through combinations of prior knowledge (l), advice distance (k), consensus (j), and sample size (i)
# Prior knowledge
for (l in c(4, 6, 8)){
  nA <- l
  
  # Distance
  for (k in seq(0.25,1.75,0.25)){
    meanB <- meanA + (k)
    
    # Consensus
    for (j in seq(0.5,1.5,0.5)){
      sdB <- j*sdA
      
      # Sample size
      for (i in 1:20){
        weightB <- i/(i+nA)
        
        myMix <- UnivarMixingDistribution(Norm(mean=meanA, sd=sdA), 
                                          Norm(mean=meanB, sd=sdB),
                                          mixCoeff=c((1-weightB), weightB))
        
        ## ... and then a function for sampling random variates from it
        rmyMix <- r(myMix)
        
        x <- rmyMix(1e6)
        meanX <- mean(x)
        sdX <- sd(x)
        madX <- mad(x)
        dat <- rbind(dat, c(meanA, sdA, nA, meanB, k, i, sdB, meanX, sdX, madX))
      }
      
    }
  }
}

# Restructure data frame and assign column names
dat2 <- dat
dat <- as.data.frame(dat2)
names(dat) <- c("meanA", "sdA", "nA", "meanB", "distance", "sample", "sdB", "meanX", "sdX", "madX")

write.table(dat, file = paste0("simA.csv"), row.names = FALSE)

# Figure 1 and 1A ---------------------------------------------------------
rm(list = ls())

dat <- read.table(file = paste0("simA.csv"), header = TRUE)

# Set graphic parameters
pd <- position_dodge(0.1)
line <- 3
point <- 14
base <- 20
key <- 1.5

# Figure 1

samp.labs <- c("N Advice = 3", "N Advice = 6", "N Advice = 9")
names(samp.labs) <- c("3", "6", "9")

dat %>%
  filter(nA %in% c(4), sample %in% c(3, 6, 9), sdB %in% c(0.5, 1, 1.5)) %>%
  ggplot(aes(x = distance, y = sdX, linetype = as.factor(sdB))) + 
  geom_abline(intercept = 1, slope = 0, colour = "black", size = 1) +
  geom_line(size = line, colour = "#A51E37") +
  scale_x_continuous(name = "Distance of advice",
                     limits = c(0.25, 1.75), breaks = c(0:7*0.25)) +
  scale_y_continuous(name = "post-advice SD",
                     limits = c(0.0, 2.2), breaks = c(0:25*0.1)) +
  scale_linetype_manual(name="Consensus of advice",
                        breaks=c("0.5", "1", "1.5"),
                        labels=c("0.5", "1", "1.5"),
                        values=c("dotted","solid", "dashed")) +
  theme_bw(base_size = base) +
  theme(legend.key.size = unit(key,"cm")) +
  theme(axis.text = element_text(colour = "black")) +
  theme(legend.justification = c(0,0), legend.position = c(0,0), legend.direction = "horizontal", legend.background = element_rect(colour = "black", size = 1)) +
  facet_wrap(~sample, labeller = labeller(sample = samp.labs)) +
  theme(
    strip.text.x = element_text(
      size = base))

# Figure 1A

sd <- dat %>%
  select(meanA:sdX) %>%
  rename(consistency = sdX) %>%
  mutate(parameter = "sd")

mad <- dat %>%
  select(meanA:meanX,madX) %>%
  rename(consistency = madX) %>%
  mutate(parameter = "mad")

cons <- rbind(sd, mad)

cons$parameter <- factor(cons$parameter, levels = c("sd", "mad"), labels = c("SD", "MAD"))

samp.labs <- c("N Advice = 3", "N Advice = 6", "N Advice = 9")
names(samp.labs) <- c("3", "6", "9")

cons %>%
  filter(nA %in% c(4), sample %in% c(3, 6, 9), sdB %in% c(0.5, 1, 1.5)) %>%
  ggplot(aes(x = distance, y = consistency, linetype = as.factor(sdB), colour = parameter)) + 
  geom_abline(intercept = 1, slope = 0, colour = "black", size = 1) +
  geom_line(size = line) +
  scale_x_continuous(name = "Distance of advice",
                     limits = c(0.25, 1.75), breaks = c(0:7*0.25)) +
  scale_y_continuous(name = "post-advice consistency",
                     limits = c(0.0, 2.2), breaks = c(0:25*0.1)) +
  scale_colour_manual(name="Parameter",
                      # breaks=c("sd", "mad"),
                      # labels=c("SD", "MAD"),
                      values=c("#A51E37", "#32414B")) +
  scale_linetype_manual(name="Consensus of advice",
                        breaks=c("0.5", "1", "1.5"),
                        labels=c("0.5", "1", "1.5"),
                        values=c("dotted","solid", "dashed")) +
  theme_bw(base_size = base) +
  theme(legend.key.size = unit(key,"cm")) +
  theme(axis.text = element_text(colour = "black")) +
  theme(legend.justification = c(0,0), legend.position = c(0,0), legend.direction = "horizontal", legend.background = element_rect(colour = "black", size = 1)) +
  facet_wrap(~sample, labeller = labeller(sample = samp.labs)) +
  theme(
    strip.text.x = element_text(
      size = base))

# Simulation in S1 (pieces of information) --------------------------------
rm(list = ls())

# Create empty data frame to store results
dat <- NULL

# Set parameters for prior knowledge
A <- c(-1, 0, 1)
meanA <- mean(A)
sdA <- sd(A)
madA <- mad(A)
nA <- 3

# Set seed to make reproducible
set.seed(20201012)

# Loop through combinations of advice distance (k), consensus (j), and sample size (i)
for (m in 1:10){
  ## Loop
  for (l in 1:1000){
    # Prior knowledge
    # Distance
    for (k in seq(0.25,1.75,0.25)){
      meanB <- meanA + k
      
      # Consensus
      for (j in c(0.5, 1, 1.5)){
        consensus <- j
        sdB <- j*sdA
        advice <- rnorm(9, mean = meanB, sd = sdB)
        
        # Sample size
        for (i in c(1, 3, 6, 9)){
          A <- c(A, advice[1:i])
          
          meanX <- mean(A)
          sdX <- sd(A)
          madX <- mad(A)
          dat <- rbind(dat, c(meanA, sdA, madA, nA, meanB, k, i, consensus, sdB, meanX, sdX, madX))
        }
        # Reset prior
        A <- c(-1, 0, 1)
        
      }
    }
  }
  
  # Restructure data frame and assign column names
  dat2 <- dat
  dat <- as.data.frame(dat2)
  names(dat) <- c("meanA", "sdA", "madA", "nA", "meanB", "distance", "sample", "consensus", "sdB", "meanX", "sdX", "madX")
  # Save results to save computing time
  write.table(dat, file = paste0(m, "simB.csv"), row.names = FALSE)
  dat <- NULL
}

# Get the files names for the separate simulation results
files <- list.files(pattern="*simB.csv")
# Add separate results to single file
myfiles <- do.call(rbind, lapply(files, function(x) read.table(x, header = TRUE, stringsAsFactors = FALSE)))

dat <- myfiles

write.table(dat, file = paste0("complete_simB.csv"), row.names = FALSE)

# Figure S1 ---------------------------------------------------------------
rm(list = ls())

dat <- read.table(file = paste0("complete_simB.csv"), header = TRUE)

sd <- dat %>%
  select(meanA:sdX) %>%
  rename(consistency = sdX) %>%
  mutate(parameter = "sd")

mad <- dat %>%
  select(meanA:meanX,madX) %>%
  rename(consistency = madX) %>%
  mutate(parameter = "mad")

cons <- rbind(sd, mad)

cons$parameter <- factor(cons$parameter, levels = c("sd", "mad"), labels = c("SD", "MAD"))

# Set graphic parameters
pd <- position_dodge(0.1)
line <- 3
point <- 14
base <- 20
key <- 1.5

samp.labs <- c("N Advice = 3", "N Advice = 6", "N Advice = 9")
names(samp.labs) <- c("3", "6", "9")

cons %>%
  filter(nA %in% c(3), sample %in% c(3, 6, 9), consensus %in% c(0.5, 1, 1.5)) %>%
  ggplot(aes(x = distance, y = consistency, linetype = as.factor(consensus), colour = parameter)) + 
  geom_abline(intercept = 1, slope = 0, colour = "black", size = 1) +
  geom_smooth(size = line, method = "lm", se = FALSE) +
  scale_x_continuous(name = "Distance of advice",
                     limits = c(0.25, 1.75), breaks = c(0:7*0.25)) +
  scale_y_continuous(name = "post-advice consistency",
                     limits = c(0.0, 7.0), breaks = c(0:25*0.1)) +
  scale_linetype_manual(name="Consensus of advice",
                        breaks=c("0.5", "1", "1.5"),
                        labels=c("0.5", "1", "1.5"),
                        values=c("dotted","solid", "dashed")) +
  scale_colour_manual(name="Parameter",
                      # breaks=c("sd", "mad"),
                      # labels=c("SD", "MAD"),
                      values=c("#A51E37", "#32414B")) +
  coord_cartesian(ylim = c(0, 2.2)) +
  theme_bw(base_size = base) +
  theme(legend.key.size = unit(key,"cm")) +
  theme(axis.text = element_text(colour = "black")) +
  theme(legend.justification = c(0,0), legend.position = c(0,0), legend.direction = "horizontal", legend.background = element_rect(colour = "black", size = 1)) +
  facet_wrap(~sample, labeller = labeller(sample = samp.labs)) +
  theme(
    strip.text.x = element_text(
      size = base)) +
  guides(linetype = guide_legend(override.aes= list(colour = "black"), order = 2), 
         colour = guide_legend(order = 1))