Analysis Plan

# ---------- Mr. A vs Mr. B2 Scenario ---------- #

# Coding legend: response: 0 (Happier/More upset); 1 (Less Happy/Less upset) 2
# (No difference)

# Create dataframe scenario 1
data_s1 <- data_short %>%
    select(subject, contains("A", ignore.case = FALSE)) %>%
    # Recode response
mutate(`gain-gain VS gain` = case_when(A1 == "A" ~ 0, A1 == "B" ~ 1, T ~ 2), `loss-loss VS loss` = case_when(A2 ==
    "A" ~ 0, A2 == "B" ~ 1, T ~ 2), `gain-loss VS gain` = case_when(A3 == "A" ~ 0,
    A3 == "B" ~ 1, T ~ 2), `loss-gain VS loss` = case_when(A4 == "A" ~ 0, A4 == "B" ~
    1, T ~ 2)) %>%
    select(-contains("A", ignore.case = FALSE)) %>%
    melt(id.var = "subject", variable.name = "scenario", value.name = "response") %>%
    mutate(scenario_group = ifelse(scenario == "gain-gain VS gain" | scenario ==
        "gain-loss VS gain", "gain", "loss"))

data_s1 %>%
    ggplot(aes(response)) + geom_bar(aes(y = (..count..)/sum(..count..)), fill = jcolors::jcolors(palette = "pal6")[1]) +
    mytheme() + theme(legend.position = "none", axis.title.y = element_blank(), axis.text.y = element_blank()) +
    scale_fill_brewer(palette = "Set1") + scale_y_continuous(labels = scales::percent,
    guide = "prism_offset") + scale_x_continuous(guide = "prism_offset", breaks = 0:2,
    labels = c("A", "B", "Same")) + labs(y = NULL, x = NULL) + facet_wrap(~scenario)

functions {
  /* 
   * Args: 
   *   y:  response value 0: A, 1: B, 2: No difference
   *   ds: different (0|1) - same (2) probability
   *   cA: conditional (different) A probability
   * Returns:  
   *   a scalar to be added to the log posterior 
   */ 
   real mental_accounting_lpdf(real y,  real ds, real cA) {
     if (y == 0) { // Stay (Explore)
       return bernoulli_lpmf(1 | ds) + bernoulli_lpmf(1 | cA); 
     } else if (y == 1) { // Switch (Explore)
       return bernoulli_lpmf(1 | ds) + bernoulli_lpmf(0 | cA);
     } else { // Select
       return bernoulli_lpmf(0 | ds);
     }
   }
}

data {
  int<lower=1> N;  // total number of observations
  vector[N] Y;  // response variable
  int<lower=1> K_ds;  // number of population-level effects
  matrix[N, K_ds] X_ds;  // population-level design matrix
  int<lower=1> K_cA;  // number of population-level effects
  matrix[N, K_cA] X_cA;  // population-level design matrix
  
  int prior_only;  // should the likelihood be ignored?
}
parameters {
  vector[K_ds] b_ds;  // population-level effects
  vector[K_cA] b_cA;  // population-level effects
}

transformed parameters {
}

model {
  // likelihood including constants
  if (!prior_only) {
    // initialize different - same linear predictor term
    vector[N] ds = X_ds * b_ds;
    // initialize linear predictor term
    vector[N] cA = X_cA * b_cA;
    
    for (n in 1:N) {
      // apply the inverse link function
      ds[n] = inv_logit(ds[n]);
    }
    for (n in 1:N) {
      // apply the inverse link function
      cA[n] = inv_logit(cA[n]);
    }
    for (n in 1:N) {
      target += mental_accounting_lpdf(Y[n] | ds[n], cA[n]);
    }
  }
  // priors including constants
  target += normal_lpdf(b_ds | 0, 5);
  target += normal_lpdf(b_cA | 0, 5);
}

# ---- Fit Model ---- #
library(rstan)
library(tidybayes)
options(mc.cores = 4)
rstan_options(auto_write = TRUE)

# ------------- ~ scenario ------------- #
# Create Stan Data
X <- model.matrix(~ 0 + scenario, data_s1 %>% filter)
 
#Stan Data
stan_data <- list(
    N        = nrow(data_s1),
    Y        = data_s1$response,
    X_ds     = X,
    X_cA     = X,
    K_ds     = ncol(X),
    K_cA     = ncol(X),
    
    N_sbj    = n_sbj,
    J_sbj    = data_s1$subject,
    
    prior_only = 0
  )

# Fit the model
fit_S1 <- rstan::sampling(mixture_scenario_1, 
                          iter = 2000, 
                          cores = 4, 
                          data = stan_data, 
                          save_warmup = FALSE)

draws <- extract(fit_S1)
colnames(draws$b_ds) <- words_replace(colnames(X), "scenario")
colnames(draws$b_cA) <- words_replace(colnames(X), "scenario")
df <- as.data.frame(draws$b_cA)

rbind(data.frame(post = with(df, `gain-gain VS gain` - `gain-loss VS gain`), effect = "gain-gain VS gain\n-\ngain-loss VS gain"),
    data.frame(post = with(df, `loss-loss VS loss` - `loss-gain VS loss`), effect = "loss-loss VS loss\n-\nloss-gain VS loss")) %>%
    ggplot(aes(x = post, y = effect)) + stat_halfeye(fill = jcolors::jcolors(palette = "pal8")[6],
    color = jcolors::jcolors(palette = "pal8")[12]) + mytheme() + theme(axis.text.y = NULL) +
    geom_vline(xintercept = 0, linetype = 2) + labs(y = NULL, x = "Posterior Probability")

# --------- Coding legend: --------- #
# -- Dependent Variable -- #
# response: 0
#           5 
#           10 
#           Other
# -- Explanatory Variables -- #
# cost: 0, 5, 10 (defined as p in Thaler 1985) 
#      [cost influence the "fair price" (p_star), which is operationalized
#       with the price asked to a friend]
# market_value : 5, 10
# buyer: friend (proxy for a fair price), stranger


# Create dataframe scenario 2
df <- data_short_MA %>% 
  select(contains('B', ignore.case =  F))

col_names <- names(df)
cost <- rep(c(0, 5, 10), each=2)
market_value <- rep(c(5, 10), 3)
col_idx <- seq(1,ncol(df),2)
df1 <- map_df(seq_along(col_idx), function(i){
  
  data.frame( response = c(df[,col_names[col_idx[i]]], df[,col_names[col_idx[i]+1]]),
              cost = cost[i],
              market_value = market_value[i],
              buyer = c(rep('Friend', n_sbj), rep('Stranger', n_sbj)) )
  
}) 

data_s2 <- df1 %>% 
  mutate(response=ifelse(response%in%unique(cost), response, 'Other'),
         response = factor(response, levels = c('0', '5', '10', 'Other')))

data_s2%>% 
  group_by(cost, market_value, buyer) %>% 
  count(response, .drop = F) %>% 
  ggplot(aes(response, n, fill=buyer)) +
  geom_bar(stat="identity", position=position_dodge(), color="black") + 
  mytheme() + 
  scale_fill_jcolors(palette = 'pal6') +
  labs(y=NULL, x=NULL) +
  facet_grid(cost~market_value)

data_s2 <- df1 %>%
    mutate(response = ifelse(response >= market_value, 1, 0))

data_s2 %>%
    group_by(buyer) %>%
    summarise(m = mean(response)) %>%
    ggplot(aes(buyer, m, fill = buyer)) + geom_bar(stat = "identity", position = position_dodge(),
    color = "black") + scale_fill_jcolors(palette = "pal6") + labs(y = "P(Price ≥ Market Value)",
    x = "Buyer") + mytheme() + theme(legend.position = "none") + scale_y_continuous(limits = 0:1)

bernoulli_scenario_2 <- rstan::stan_model("bernoulli_scenario_4.stan")

# Create Stan Data
X <- model.matrix(~ buyer, data_s2)
  
#Stan Data
stan_data <- list(
  N       = nrow(data_s2),
  Y       = data_s2$response,
  K       = ncol(X),
  X       = X,
  
  prior_only = 0
)

# Fit the model
fit_S2 <- rstan::sampling(bernoulli_scenario_2, 
                          iter = 2000, 
                          cores = 4, 
                          data = stan_data, 
                          save_warmup = FALSE)  

extract(fit_S2)$b %>%
    melt(value.name = "post") %>%
    mutate(Var2 = ifelse(Var2 == 1, "Intercept", "")) %>%
    filter(Var2 != "Intercept") %>%
    ggplot(aes(x = post, y = Var2)) + stat_halfeye(fill = jcolors::jcolors(palette = "pal8")[6],
    color = jcolors::jcolors(palette = "pal8")[12]) + mytheme() + geom_vline(xintercept = 0,
    linetype = 2) + labs(y = NULL, x = "Posterior Probability", title = "Stranger - Friend")

# --------- Coding legend: --------- # -- Dependent Variable -- # response --
# Explanatory Variables -- # store: Resort hotel vs Grocery store

# Create dataframe scenario 3
data_s3 <- data_short_MA %>%
    select(subject, contains("C", ignore.case = FALSE)) %>%
    melt(id.var = "subject", variable.name = "store", value.name = "response") %>%
    mutate(store = ifelse(store == "C1_1", "Resort Hotel", "Grocery Store"), response = as.numeric(response))

data_s3 %>%
    ggplot(aes(response, fill = store, color = store)) + geom_density(alpha = 0.4,
    adjust = 1.5, size = 2) + mytheme() + theme(legend.position = c(0.7, 0.8)) +
    scale_fill_jcolors(palette = "pal6") + scale_color_jcolors(palette = "pal6") +
    # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') +
scale_x_continuous(breaks = seq(0, 40, 10), limits = c(0, 40)) + labs(y = NULL, x = NULL)

data {
  int<lower=1> N;  // total number of observations
  vector[N] Y;  // response variable
  int<lower=1> K;  // number of population-level effects
  matrix[N, K] X;  // population-level design matrix
  
  int prior_only;  // should the likelihood be ignored?
}

parameters {
  vector[K] b;  // population-level effects
  real<lower=0> shape;  // shape parameter
}

transformed parameters {
}

model {
  // likelihood including constants
  if (!prior_only) {
    // initialize linear predictor term
    vector[N] mu = X * b;
    for (n in 1:N) {
      // apply the inverse link function
      mu[n] = shape * exp(-(mu[n]));
    }
    target += gamma_lpdf(Y | shape, mu);
  }
  // priors including constants
  target += normal_lpdf(b | 0, 5);
  target += gamma_lpdf(shape | 0.01, 0.01);
}
 

# Create Stan Data
X <- model.matrix(~ store, data_s3) 

#Stan Data
stan_data <- list(
  N = nrow(data_s3),
  Y = data_s3$response,
  X = X,
  K = ncol(X),
  
  N_sbj    = n_sbj,
  J_sbj    = data_s3$subject,
  
  prior_only = 0
)

# Fit the model
fit_S3 <- rstan::sampling(gamma_scenario_3, 
                          iter = 2000, 
                          cores = 4, 
                          data = stan_data, 
                          save_warmup = FALSE)

draws <- extract(fit_S3)

draws$b %>%
    melt(value.name = "post") %>%
    mutate(Var2 = ifelse(Var2 == 1, "Intercept", "")) %>%
    filter(Var2 != "Intercept") %>%
    ggplot(aes(x = post, y = Var2)) + stat_halfeye(fill = jcolors::jcolors(palette = "pal8")[6],
    color = jcolors::jcolors(palette = "pal8")[12]) + mytheme() + geom_vline(xintercept = 0,
    linetype = 2) + labs(y = NULL, x = "Posterior Probability", title = "Resort Hotel - Grocery Store")

# --------- Coding legend: --------- #
  # -- Dependent Variable -- #
  # response 
  # -- Explanatory Variables -- #
  # store: Resort hotel vs Grocery store
  
  # Create dataframe scenario 3
  data_s4 <- data_short_MA %>% 
    select(subject, contains('D', ignore.case = FALSE)) %>% 
    rename(low_price2 = D3, low_price1 = D2_1,
           high_price2 = D4, high_price1 = D2_2) %>% 
    melt(id.var='subject',
         variable.name='price',
         value.name='response') %>% 
    mutate(price=ifelse(str_detect(price, 'low'), 'low', 'high'),
           response=ifelse(response=='No', 0, 1))
  
  data_s4 %>% 
    filter(response==1) %>% 
    group_by(price) %>%
    count(response, .drop = F) %>% 
    ungroup() %>% 
    ggplot(aes(price, n, fill=price)) +
    geom_bar(stat = "identity", 
             # aes(y = (..count..)/sum(..count..)),
             position = position_dodge()) +
    mytheme() + theme(legend.position = 'none') +
    scale_fill_jcolors(palette = 'pal6') +
    scale_color_jcolors(palette = 'pal6') +
    # scale_y_continuous(labels=scales::percent, guide = "prism_offset") +
    labs(y='Number of "Yes"', x='Price') 

bernoulli_scenario_4 <- rstan::stan_model("bernoulli_scenario_4.stan")

# Create Stan Data
X <- model.matrix(~ price, data_s4) 

#Stan Data
stan_data <- list(
  N = nrow(data_s4),
  Y = data_s4$response,
  X = X,
  K = ncol(X),
  
  N_sbj    = n_sbj,
  J_sbj    = data_s4$subject,
  
  prior_only = 0
)

# Fit the model
fit_S4 <- rstan::sampling(bernoulli_scenario_4, 
                          iter = 2000, 
                          cores = 4, 
                          save_warmup = FALSE,
                          data = stan_data, 
                          save_warmup = FALSE)

extract(fit_S4)$b %>%
    melt(value.name = "post") %>%
    mutate(Var2 = ifelse(Var2 == 1, "Intercept", "")) %>%
    filter(Var2 != "Intercept") %>%
    ggplot(aes(x = post, y = Var2)) + stat_halfeye(fill = jcolors::jcolors(palette = "pal8")[6],
    color = jcolors::jcolors(palette = "pal8")[12]) + mytheme() + geom_vline(xintercept = 0,
    linetype = 2) + labs(y = NULL, x = "Posterior Probability", title = "Low - High")

# --------- Coding legend: --------- #
# -- Dependent Variable -- #
# response 
# -- Explanatory Variables -- #
# loss: ticket vs cash

# Create dataframe scenario 3
data_s5 <- data_short_MA %>% 
  select(subject, contains('E', ignore.case = FALSE)) %>%
  rename(ticket = E1, 
         cash = E2) %>% 
  melt(id.var='subject',
       variable.name='loss',
       value.name='response') %>% 
  mutate(response=ifelse(response=='No', 0, 1))

data_s5 %>% 
  filter(response==1) %>% 
  group_by(loss) %>%
  count(response, .drop = F) %>% 
  ungroup() %>% 
  ggplot(aes(loss, n, fill=loss)) +
  geom_bar(stat = "identity", 
           # aes(y = (..count..)/sum(..count..)),
           position = position_dodge()) +
  mytheme() + theme(legend.position = 'none') +
  scale_fill_jcolors(palette = 'pal6') +
  scale_color_jcolors(palette = 'pal6') +
  # scale_y_continuous(labels=scales::percent, guide = "prism_offset") +
  labs(y='Number of "Yes"', x='Loss') 

# Create Stan Data
X <- model.matrix(~ loss, data_s5) 

#Stan Data
stan_data <- list(
  N = nrow(data_s5),
  Y = data_s5$response,
  X = X,
  K = ncol(X),
  
  N_sbj    = n_sbj,
  J_sbj    = data_s3$subject,
  
  prior_only = 0
)

# Fit the model
fit_S5 <- rstan::sampling(bernoulli_scenario_4, 
                          iter = 2000, 
                          cores = 4, 
                          save_warmup = FALSE,
                          data = stan_data, 
                          save_warmup = FALSE)

extract(fit_S5)$b %>%
    melt(value.name = 'post') %>% 
    mutate(Var2=ifelse(Var2==1, 'Intercept', '')) %>% 
    filter(Var2!='Intercept') %>%
    ggplot(aes(x = post, y = Var2)) + 
    stat_halfeye(fill = jcolors::jcolors(palette = "pal8")[6], 
                 color = jcolors::jcolors(palette = "pal8")[12]) + 
    mytheme() + geom_vline(xintercept = 0, linetype = 2) + 
    labs(y = NULL, x = "Posterior Probability", title='Ticket - Cash') 

answers = c("I feel like I wasted $20", "I feel like I wasted something but no specific amount or measure comes to mind",
    "I feel like I wasted nothing, since my visit had already been paid for")
# Create dataframe scenario 3
data_s6 <- data_short_MA %>%
    select(subject, contains("F", ignore.case = FALSE)) %>%
    mutate(F1 = case_when(F1 == answers[1] ~ 1, F1 == answers[2] ~ 0, F1 == answers[3] ~
        0), F2 = case_when(F2 == answers[1] ~ 1, F2 == answers[2] ~ 0, F2 == answers[3] ~
        0)) %>%
    rename(`Per-session` = F1, Yearly = F2) %>%
    melt(id.var = "subject", variable.name = "frame", value.name = "response")

data_s6 %>%
    filter(response == 1) %>%
    group_by(frame) %>%
    count(response, .drop = F) %>%
    ungroup() %>%
    ggplot(aes(frame, n, fill = frame)) + geom_bar(stat = "identity", position = position_dodge()) +
    mytheme() + theme(legend.position = "none") + scale_fill_jcolors(palette = "pal6") +
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + scale_color_jcolors(palette
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + =
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + "pal6")
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + +
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + #
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + scale_y_continuous(labels=scales::percent,
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + guide
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + =
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + 'prism_offset')
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + +
labs(y = "Number of \"Wasted $20\"", x = "Frame")

bernoulli_scenario_6 <- rstan::stan_model("bernoulli_scenario_6.stan")

# Create Stan Data
X <- model.matrix(~ frame, data_s6 %>% mutate(frame=ifelse(frame=='Per-session', 1, 0 )))

#Stan Data
stan_data <- list(
  N = nrow(data_s6),
  Y = data_s6$response,
  X = X,
  K = ncol(X),
  
  N_sbj    = n_sbj,
  J_sbj    = data_s6$subject,
  
  prior_only = 0
)

# Fit the model
fit_S6 <- rstan::sampling(bernoulli_scenario_6, 
                          iter = 2000, 
                          cores = 4, 
                          pars    = c('b'),
                          save_warmup = FALSE,
                          data = stan_data, 
                          save_warmup = FALSE)

extract(fit_S6)$b %>%
    as.data.frame() %>%
    mutate(post = V2, type = "") %>%
    ggplot(aes(x = post, y = type)) + stat_halfeye(fill = jcolors::jcolors(palette = "pal8")[6],
    color = jcolors::jcolors(palette = "pal8")[12]) + mytheme() + geom_vline(xintercept = 0,
    linetype = 2) + labs(y = NULL, x = "Posterior Probability", title = "Per-session - Yearly")

# --------- Coding legend: --------- # -- Dependent Variable -- # response --
# Explanatory Variables -- # coupon: purchased vs free

answers = c("Pay your friend $35 for the coupon.", "Pay some, but not the full amount for the coupon (for example, half the price).",
    "Consider it a gift and not pay for the coupon.")
# Create dataframe scenario 3
data_s7 <- data_short_MA %>%
    select(subject, contains("G", ignore.case = FALSE)) %>%
    mutate(G1 = case_when(G1 == answers[1] ~ 1, G1 == answers[2] ~ 0, G1 == answers[3] ~
        0), G2 = case_when(G2 == answers[1] ~ 1, G2 == answers[2] ~ 0, G2 == answers[3] ~
        0)) %>%
    rename(purchased = G1, free = G2) %>%
    melt(id.var = "subject", variable.name = "coupon", value.name = "response")

data_s7 %>%
    filter(response == 1) %>%
    group_by(coupon) %>%
    count(response, .drop = F) %>%
    ungroup() %>%
    ggplot(aes(coupon, n, fill = coupon)) + geom_bar(stat = "identity", position = position_dodge()) +
    mytheme() + theme(legend.position = "none") + scale_fill_jcolors(palette = "pal6") +
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + scale_color_jcolors(palette
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + =
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + "pal6")
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + +
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + #
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + scale_y_continuous(labels=scales::percent,
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + guide
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + =
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + 'prism_offset')
    scale_color_jcolors(palette = "pal6") + # scale_y_continuous(labels=scales::percent, guide = 'prism_offset') + +
labs(y = "Pay your friend $35 for the coupon", x = "")

# Create Stan Data
  X <- model.matrix(~ coupon, data_s7)
  
  #Stan Data
  stan_data <- list(
    N = nrow(data_s7),
    Y = data_s7$response,
    X = X,
    K = ncol(X),
    
    N_sbj    = n_sbj,
    J_sbj    = data_s7$subject,
    
    prior_only = 0
  )
  
  # Fit the model
  fit_S7 <- rstan::sampling(bernoulli_scenario_6, 
                            iter = 2000, 
                            cores = 4, 
                            pars    = c('b'),
                            save_warmup = FALSE,
                            data = stan_data, 
                            save_warmup = FALSE)

extract(fit_S7)$b %>%
    as.data.frame() %>%
    mutate(post = V2, type = "") %>%
    ggplot(aes(x = post, y = type)) + stat_halfeye(fill = jcolors::jcolors(palette = "pal8")[6],
    color = jcolors::jcolors(palette = "pal8")[12]) + mytheme() + geom_vline(xintercept = 0,
    linetype = 2) + labs(y = NULL, x = "Posterior Probability", title = "Per-session - Yearly")