---
title: "Money vs. time in interpersonal choice "
Author: "Alexander Pastukhov and Johannes Leder"
output: html_document
---


## Service

```{r global_options, include=FALSE}
library(arsenal)
library(beanplot)
library(correlation)
library(cowplot)
library(DiagrammeR)
library(DiagrammeRsvg)
library(glmmTMB)
library(gridExtra)
library(kableExtra)
library(effectsize)
library(lmerTest)
library(lsr)
library(patchwork)
library(pscl)
library(psych)
library(readxl)
library(rsvg)
library(skimr)
library(tidyverse)
```

```{r Loading fonts}
extrafont::loadfonts(quiet = TRUE)
```

```{r Custom plotting theme}
plot_theme <- function(){
  theme(text=element_text(family="Arial"), 
        axis.text.x = element_text(size = 8, colour = 'black'), 
        axis.text.y = element_text(size = 8, colour = 'black'), 
        axis.title.x = element_text(size = 10), 
        axis.title.y = element_text(size = 10), 
        panel.grid.minor.x = element_blank(), 
        panel.grid.minor.y =  element_line(size = 0.24), 
        axis.ticks = element_line(size = 0.24),
        plot.subtitle = element_text(size = 10))
}
```




```{r custom functions}
# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

```

## Getting Data

```{r}
data <- 
  fs::dir_ls(path="../Data", glob = "*.xlsx") %>% 
  map_df(~read_xlsx(.)) %>% 
  distinct(participant.mturk_worker_id,.keep_all = TRUE) %>% 
  mutate(pie_other = money_time_tradeoff.1.player.giveaway/money_time_tradeoff.1.player.salary,
         vpn = 1:nrow(.),
         worked_for = case_when(role_assignment.1.player.id_in_group==1 ~ "time",
                                role_assignment.1.player.id_in_group==2 ~ "money",
                                TRUE ~ "non"))%>%
  separate(timed_slider_task_copy.1.player.tasklog, remove = FALSE, into = c(NA,"blocks_completed", NA), sep=",") %>% 
  
  ## getting the number of completed blocks in the real effort task
  mutate(blocks_completed = as.numeric(stringr::str_replace(blocks_completed,'"completed":',"")),
         completed_work = timed_slider_task.1.player._nblocks == blocks_completed,
         completed_work = case_when(participant._current_page_name=="TaskPage"~FALSE,
                                    role_assignment.1.player._work_for_other=="0:00" ~ TRUE,
                                    TRUE ~ completed_work),
         quit_waiting = case_when(payment_info.1.player._payoff_dropout==2 ~ "yes",
                                  payment_info.1.player._payoff_dropout==0 ~ "no",
                                  participant._current_page_name=="RoleWaitPage"~"yes",
                                  is.na(payment_info.1.player._payoff_dropout)&
                                     participant._current_app_name=="payment_info"~"no"),
         salary_too_high = money_time_tradeoff.1.player.salary_offer>money_time_tradeoff.1.player.salary) %>% 
  filter (money_time_tradeoff.1.player.extrawork>-1)
```

Keeping only subject with responses in the central vars, and sorting out double ids due to data retrieval
```{r}
data %>% count(payment_info.1.player._payoff_dropout)
```

# Methods
## Sample

```{r}
table.sample <- arsenal::tableby(~ demographic_information.1.player.age + 
                                   demographic_information.1.player.gender + 
                                   demographic_information.1.player.education+ 
                                   demographic_information.1.player.monthly_income,
                                 numeric.stats=c("mean","sd","median","q1q3"),
                                 data = data)

mylabels <- list(demographic_information.1.player.gender = "Gender",
                 demographic_information.1.player.age = "Age, yrs",
                 demographic_information.1.player.education = "Education",
                 demographic_information.1.player.monthly_income = "Monthly Income")

summary(table.sample, labelTranslations = mylabels)
```

## Attrition during experiment

where can we loose participants after the decision? 
1. grouping page
a) quit waiting
b) assigned but left

Bonus payment for dropouts
```{r}
sum(data$payment_info.1.player._payoff_dropout==2,na.rm=TRUE)
data$payment_info.1.player._payoff_total
sum(data$payment_info.1.player.payoff==data$payment_info.1.player._payoff_total,na.rm = TRUE)
 
data %>% 
  filter(data$payment_info.1.player._payoff_dropout==2) %>% 
  select(contains("mturk")) %>% 
  write.csv("workerIDs_add_Bonus.csv")
```


### Who quit the real effort task?
```{r}
dropout_after_choice<-data  %>% count(participant._current_app_name) %>% pull(2) %>% .[2]
dropout_salary<-data %>% count(data$salary_too_high) %>% pull(2) %>% .[2]
dropout_waiting<-data %>% count(data$quit_waiting) %>% pull(2) %>% .[2] # get the second elemnt of the vector
dropout_real_effort<-data  %>% count(participant._current_app_name) %>% pull(2) %>% .[5:6] %>% sum()
dropout_followup<-data  %>% count(participant._current_app_name) %>% pull(2) %>% .[1]

```

creating a flowchart

```{r}


mermaid("
graph LR
  A[Decision n = 371]-->B[dropout n = 2]
  A-->C
  C{Salary n= 369}-->D[too high n = 138]
  C{Salary n= 369}-->E[Waiting Page n = 231]
  D[too high n = 138]-->I[Follow up n = 325]
  E[Waiting Page n = 231]-->F[not waiting n = 71]
  F[not waiting n = 71]--> L[dropout n = 38]
  F[not waiting n = 71]-->I[Follow up n = 326]
  E[Waiting Page n = 231]-->G[Real effort task n = 160]
  G[Real effort task n = 160]-->I[Follow up n = 326]
  G[Real effort task n = 160]-->H[dropout n = 5]
  I[Follow up n = 326]-->J[dropout n =1]
  I[Follow up n = 326]-->K[Payment n = 325]
")

```


```{r}
data  %>% count(participant._current_app_name)
```

The higher number of participants in the time condition is due the the default of time when arriving at the waiting page.


# Social Value Orientation

### Creating SVO scores

 
```{r}


get_svo_angle<-function (subject) {
  mean_to_self <- subject %>% select (contains("player.input_self")) %>% rowSums()/6
  mean_to_others <- subject %>% select (contains("player.input_other")) %>% rowSums()/6
  angle <- atan((mean_to_others - 50) / (mean_to_self - 50))*180/pi
  return(angle)
}

data$svo_angle<-data %>% get_svo_angle()

data<-data %>% mutate (svo_detailed_type = case_when(svo_angle < -12.04 ~ "competitive",
                                       svo_angle >= -12.04 & svo_angle < 22.45 ~ "individualist", 
                                       svo_angle >= 22.45 & svo_angle < 57.15 ~ "prosocial",
                                       svo_angle >= 57.15 ~ "altruistic"))


```

## Distribution of Social Value Orientation

```{r distribution of svo }
svo_angles <-
  data %>%
  group_by(vpn, svo_angle, svo_detailed_type) %>%
  drop_na(svo_angle) %>%
  summarise()

svo_type_counts <-
  svo_angles %>%
  ungroup() %>%
  mutate(svo_detailed_type = factor(svo_detailed_type, levels = c("competitive", "individualist", "prosocial", "altruistic"))) %>%
  group_by(svo_detailed_type) %>%
  summarise(count = n()) %>%
  ungroup() %>%
  tidyr::complete(svo_detailed_type, fill=list(count=0)) %>% 
  mutate(count = glue::glue("{svo_detailed_type}\nN={count}")) %>% 
  pull(count)

ggplot(svo_angles,
       aes(x = svo_angle)) +
  geom_histogram(aes(y=..count../sum(..count..))) + 
  geom_vline(xintercept = c(-12.04, 22.45, 57.15), color="darkgray") +
  # annotate("text", x=21.5, y=Inf,  hjust = 1, vjust = 1, family="Arial", color="#F8766D",
  #          label = glue::glue("← individualist\nN={sum(data$svo_type=='individualist', na.rm=TRUE)}"))+
  # annotate("text", x=23.5, y=Inf, hjust = 0, vjust = 1, family="Arial", color="#00BFC4",
  #          label = glue::glue("prosocial →\nN={sum(data$svo_type=='prosocial', na.rm=TRUE)}"))+
  theme_bw() + 
  theme(legend.position = "none", legend.title = element_blank()) + 
  scale_x_continuous(name="Angle from SVO slider measure", 
                     limits = c(-30, 70),
                     breaks = c(-12.04, 22.45, 57.15), 
                     labels = c("-12.04°", "22.45°", "57.15°"),
                     sec.axis = sec_axis(~ ., breaks = c(-25, 5, 40, 65), labels = svo_type_counts))+
  scale_y_continuous(name = "Proportion")+
  plot_theme()

ggsave("Figs/SVO_angle.png",units = "cm", dpi = 300, width = 10, height=10)
```

# Hexaco Scoring


```{r}
hexaco_items<-data %>% select(contains("hex")) %>% colnames()

list_items_hexaco <-
  read.csv2("Hexaco60Items.csv") %>% 
  pivot_longer(cols = Items1:Items3,values_to = "Item") %>% 
  select(-name) %>% 
  filter(Item != "") %>% 
  
  # creating a code for reversed items
  mutate(reverse = stringr::str_sub(Item,-1)) %>%
  
  # only keeping the R 
  mutate(reverse = case_when(reverse =="R" ~ reverse,
                             TRUE ~ "")) %>%
  
  # removing the R from item name
  mutate(Item = as.numeric(stringr::str_remove(Item,"R"))) %>%
  
  # Item names in the questionnaire extracted from data
  mutate(Item_code = paste0("hex_",Item))


items_reverse <- 
  list_items_hexaco %>% 
  filter (reverse == "R") %>% 
  select(Item_code) %>% 
  mutate(Item_code = paste0("questionnaires.1.player.",Item_code)) %>% 
  pull()

honesty_humility <-
  list_items_hexaco %>% 
  filter(Scale=="Honesty-Humility") %>% 
  pull(Item_code)

data<-data %>% 
  mutate_at(.vars=items_reverse, ~6-.) %>% 
  mutate(honesty_humility_score = reduce(select(., contains(honesty_humility)), `+`)/10)

data %>% 
  select(contains(honesty_humility)) %>% 
  psych::alpha()
```

```{r}
ggplot(data,
       aes(x = honesty_humility_score)) +
  geom_histogram(aes(y=..count../sum(..count..))) + 
  # geom_density()+
  theme_bw() + 
  theme(legend.position = "none", legend.title = element_blank()) + 
  scale_y_continuous(name = "Proportion")+
  plot_theme()
```

# Moral Identity Centrality

```{r}
miq_items<- 
  data %>% 
  select(contains("miq")) %>% 
  colnames()

data %>% 
  select(contains("miq")) %>% 
  psych::alpha()

data$miq_score <- 
  data %>% 
  select(contains("miq")) %>% 
  rowMeans()
```

```{r}
ggplot(data,
       aes(x = miq_score)) +
  geom_histogram(aes(y=..count../sum(..count..))) + 
  # geom_density()+
  theme_bw() + 
  theme(legend.position = "none", legend.title = element_blank()) + 
  scale_y_continuous(name = "Proportion")+
  plot_theme()
```

# Positive and negative affect 

```{r}
panas_items <- 
  data %>% 
  select (contains("PANAS ")) %>% 
  colnames()

list_items_panas<-read.csv2("PanasItems.csv") %>% 
  pivot_longer(cols = item1:item10,values_to = "Item") %>% 
  select (-name) %>% 
  
  # Item names in the questionnaire extracted from data
  mutate (Item_code = paste0("panas_measure.1.player.PANAS ",Item))

panas_positive <- 
  list_items_panas %>% 
  filter(subscale=="positive") %>% 
  pull(Item_code)

panas_negative <- 
  list_items_panas %>% 
  filter(subscale=="negative") %>% 
  pull(Item_code)


data %>% 
  select(panas_positive) %>% 
  psych::alpha()

data %>% 
  select(panas_negative) %>% 
  psych::alpha()

data$panas_positive <- 
  data %>% 
  select(panas_positive) %>% 
  rowSums()

data$panas_negative <- 
  data %>% 
  select(panas_negative) %>% 
  rowSums()
```

# Generosity index
* 'proportion of time OR money given to the other (Davis et al. 2015)'
* 'expressing time in money to have a comparison based on personal cost : The wage per second, which reflects the
opportunity cost of time, will be multiplied by the amount of time the participant decided he
or she would work for charity expressed in seconds to derive a common scale for both
resources.'

Problem: if the salary cannot be divided by 10 to have a even number, it rounded funny.

```{r generosity index}
data <- 
  data %>%
  mutate (generosity_index_money = round(case_when(
            money_time_tradeoff.1.player.salary > 0 ~
            money_time_tradeoff.1.player.giveaway /
            money_time_tradeoff.1.player.salary,
            TRUE ~ money_time_tradeoff.1.player.giveaway / 10),1),
          
    generosity_index_time = ((money_time_tradeoff.1.player.extrawork - 1) * 90) / 900,
    
    # salary per second multiplied with the time given in seconds
    generosity_index_time_as_money = (money_time_tradeoff.1.player.salary / 900 * ((money_time_tradeoff.1.player.extrawork - 1) * 90)) / money_time_tradeoff.1.player.salary,
    
    # salary per second multiplied with the time given in seconds
    generosity_as_money = money_time_tradeoff.1.player.giveaway,
    
    generosity_time_as_money = (money_time_tradeoff.1.player.salary / 900) * ((money_time_tradeoff.1.player.extrawork - 1) * 90),
    
    diff_time_money_index =  generosity_index_time-generosity_index_money,
    
    give_raw_money = case_when(money_time_tradeoff.1.player.salary > 0 ~ money_time_tradeoff.1.player.giveaway * 10/money_time_tradeoff.1.player.salary,
                               TRUE ~ money_time_tradeoff.1.player.giveaway * 10 / 10),
    
    give_raw_time = money_time_tradeoff.1.player.extrawork - 1
  )
```


```{r long data}
data.long <-
  data %>% 
  select(vpn,
         give_raw_time,
         give_raw_money,
         generosity_index_money,
         generosity_index_time,
         diff_time_money_index,
         money_time_tradeoff.1.player.salary,
         miq_score,
         honesty_humility_score,
         svo_angle,
         svo_detailed_type,
         panas_positive) %>% 
  pivot_longer(-c(miq_score, honesty_humility_score, svo_angle, svo_detailed_type, panas_positive, vpn, money_time_tradeoff.1.player.salary, diff_time_money_index),
               names_to= "measure", values_to = "score") %>%
  
  mutate(condition = stringr::str_sub(measure,-5,-1),
         condition = stringr::str_replace(condition,"_",""),
         svo_angle_centered = svo_angle-mean(svo_angle,na.rm = TRUE),
         honesty_humility_centered= honesty_humility_score-mean(honesty_humility_score,na.rm = TRUE),
         miq_centered = miq_score-mean(miq_score,na.rm = TRUE)) %>% 
  
  separate(measure,into = c(NA,"measure",NA), sep = "_") %>% 
  rename (req_salary = money_time_tradeoff.1.player.salary) %>% 
  pivot_wider(names_from = measure, values_from = score) %>% 
  mutate(raw = as.integer(raw),
         non_zero = index!=0,
         index_beta = abs(index-.000001),
         index_beta_trunc = case_when(index == 1 ~ index-.000001,
                                      TRUE ~ index)) %>% 
  drop_na(miq_score)
```


# Analysis

## Correlation between all variables
```{r correlations}

##https://easystats.github.io/correlation/
correlations <-
  data %>% select(
    generosity_index_money,
    give_raw_money,
    generosity_index_time,
    demographic_information.1.player.monthly_income,
    give_raw_time,
    generosity_time_as_money,
    generosity_index_time_as_money,
    diff_time_money_index,
    miq_score,
    honesty_humility_score,
    svo_angle)

kable(summary(correlation(correlations)),digits = 2) %>% save_kable(file = "correlationTable.html", self_contained = T)
```

## Plotting stuff
```{r}
ggplot(data, aes(x = demographic_information.1.player.monthly_income, y = give_raw_money)) + 
  geom_point() + 
  scale_x_log10() + 
  geom_density_2d()

ggplot(data, aes(x = demographic_information.1.player.monthly_income, y = give_raw_time)) + 
  geom_point() + 
  scale_x_log10() + 
  geom_density_2d()

ggplot(data, aes(x = demographic_information.1.player.monthly_income, y = diff_time_money_index)) + 
  geom_point() + 
  scale_x_log10() + 
  geom_density_2d()
```

```{r}
prop_index <-
  data.long %>% 
  group_by(condition,index) %>% 
  summarise(n(),nrow(.)/2,prop = n()/(nrow(.)/2))

prop_index_opp_cost <-
  data %>% 
  group_by(generosity_index_time_as_money) %>% 
  summarise(n(),nrow(.)/2,prop = n()/(nrow(.)/2))
```

```{r}
ggplot(data,aes(x = generosity_index_time-generosity_index_money)) + 
  geom_histogram(bins = 50) + 
  facet_grid(~miq_score>3.5)
```

```{r}
ggplot(prop_index_opp_cost,aes(x = as.factor(generosity_index_time_as_money),y=prop)) + 
  geom_bar(stat = 'identity') +
  scale_y_continuous("Proportion") +
  scale_x_discrete("Generosity time as money")
```


```{r}
data %>% 
  filter(money_time_tradeoff.1.player.extrawork>-1) %>% 
  ggplot(aes(x = money_time_tradeoff.1.player.salary)) + 
  geom_histogram(aes(y=..count../sum(..count..))) + 
  scale_x_log10("Log (req. Salary)") +
  scale_y_continuous("Proportion") 

ggsave("Figs/salary.png", units = "cm", dpi = 300, width = 10,height=10)
```




```{r}
ggplot(prop_index %>% filter(condition == "money"),aes(x = as.factor(index),y=prop)) + 
  geom_bar(stat = 'identity') +
  scale_y_continuous("Proportion") +
  scale_x_discrete("Generosity for money")
ggsave("Figs/given_money.png",units = "cm", dpi = 300, width = 10,height=10)
```


```{r}
plot_pie_money <- 
  ggplot(prop_index %>% filter(condition == "money"),aes(x = as.factor(index),y=prop)) + 
  geom_bar(stat = 'identity') +
  scale_y_continuous("Proportion",limits = c(0,.5)) +
  scale_x_discrete("Generosity for money")
```


```{r}
plot_pie_time <- 
  ggplot(prop_index %>% filter(condition=="time"),aes(x = as.factor(index),y=prop)) + 
  geom_bar(stat = 'identity') +
  scale_y_continuous("Proportion",,limits = c(0,.5)) +
  scale_x_discrete("Generosity for time")
```

```{r}
ggplot(prop_index %>% filter(condition=="time"),aes(x = as.factor(index),y=prop)) + 
  geom_bar(stat = 'identity') +
  scale_y_continuous("Proportion") +
  scale_x_discrete("Generosity for time")

ggsave("Figs/given_time.png",units = "cm", dpi = 300, width = 10, height=10)
```

```{r}
ggplot(data,aes(x =generosity_index_money, y = generosity_index_time))+ 
  geom_point(position = position_jitter())+
  coord_equal()+
  geom_density2d()+
  geom_abline(slope = 1)+
  scale_x_continuous("Generosity for money") +
  scale_y_continuous("Generosity for time")+
  geom_smooth(method = "lm")

ggsave("Figs/given_time_vs_money.png", units = "cm", dpi = 300, width = 10,height=10)
```

```{r}
ggplot(data,aes(x =generosity_index_time_as_money, y = generosity_index_money))+ 
  geom_point(position = position_jitter())+
  geom_smooth(method = "lm", formula = y ~ x)+ 
   coord_equal()+
  geom_density2d()+
  geom_abline(slope = 1)
```


```{r generosity time vs. opportunity cost}
ggplot(data %>% filter(money_time_tradeoff.1.player.salary<50),aes(x =money_time_tradeoff.1.player.salary, y = generosity_index_time))+ 
  geom_point(position = position_jitter())+
 # coord_equal()+
#  geom_density2d()+
 # geom_abline(slope = 1)+
  geom_smooth(method = "lm")
 # scale_x_continuous("Generosity for money") +
#  scale_y_continuous("Generosity for time")
```


```{r}
count_per_choice <-
  data %>%
  group_by(generosity_index_money, generosity_index_time) %>%
  summarise(Count= n(),
            MedianSVO = median(svo_angle)) %>%
  ungroup() %>%
  mutate(Freq = 100 * Count / sum(Count))


plot_time_vs_money <-
  ggplot(data, aes(x=generosity_index_money, y = generosity_index_time))+ 
  # geom_point(position = position_jitter(height = 0.01, width=0.01), alpha = .4, aes(color=svo_angle, fill=svo_angle))+
  geom_abline(slope = 1, color="white", size=2)+
  geom_density2d(show.legend = FALSE) +
  geom_point(data=count_per_choice, aes(size=Freq)) +
  coord_equal()+
  scale_x_continuous("Generosity for money", breaks=seq(0, 1, 0.1)) +
  scale_y_continuous("Generosity for time", breaks=seq(0, 1, 0.1)) + 
  theme(legend.position = c(0.9, 0.2)) +
  labs(size = "Proportion") +
  guides(size = guide_legend(reverse=TRUE)) +
  plot_theme()

plot_time_vs_money
```

```{r}
# layout<- matrix(c(1,1,1,2,2,2,3,3,3,3), nrow = 1, byrow = TRUE)
# multiplot(plot_pie_money,plot_pie_time,plot_time_vs_money, layout=layout)
```


```{r}
layout <- cbind(rbind(c(1,1),c(2,2)),rbind(c(3,3),c(3,3)))
png("Figs/all_index_plots.png",units = "cm", res = 300, width = 15,height=10)
grid.arrange(plot_pie_money,plot_pie_time,plot_time_vs_money, layout_matrix = layout) 
dev.off()
```

```{r}
left_col <- plot_grid(plot_pie_money,plot_pie_time,ncol = 1, labels = c("A","B"))
png("Figs/all_index_plots.png",units = "cm", res = 300, width = 15,height=10)
plot_grid(left_col,plot_time_vs_money,labels = c("","C"))
dev.off()
```


```{r}
ggplot(data %>% filter(!is.na(svo_detailed_type)),aes(x =generosity_index_money, y = generosity_index_time))+ 
  geom_point(position = position_jitter())+
  coord_equal()+
  geom_density2d()+
  geom_abline(slope = 1)+
  scale_x_continuous("Generosity for money") +
  scale_y_continuous("Generosity for time")+
  facet_wrap(~svo_detailed_type)

ggsave("Figs/svo_type_given_time_vs_money.png",units = "cm", dpi = 300, width = 15,height=10)
```


```{r}
ggplot(data %>% filter(!is.na(svo_detailed_type)),aes(x =generosity_index_money, y = generosity_index_time))+ 
  # geom_density_2d_filled(contour_var = "ndensity") +
  geom_density_2d()+
  geom_point(position = position_jitter(), color = "white")+
  coord_equal()+
  geom_density2d()+
  geom_abline(slope = 1)+
  scale_x_continuous("Generosity for money") +
  scale_y_continuous("Generosity for time")+
  facet_wrap(~svo_detailed_type)
```



```{r}
ggplot(data,aes(x =svo_angle, y = generosity_index_time))+ 
  # geom_density_2d_filled(contour_var = "ndensity") +
  geom_density_2d()+
  geom_point(position = position_jitter(), color = "white")+
  scale_x_continuous("SVO angle") +
  scale_y_continuous("Generosity for time")
```

```{r}
ggplot(data,aes(x =svo_angle, y = generosity_index_money))+ 
  # geom_density_2d_filled(contour_var = "ndensity") +
  geom_density_2d()+
  geom_point(position = position_jitter(), color = "white")+
  scale_x_continuous("SVO angle") +
  scale_y_continuous("Generosity for money")
```


```{r}
ggplot(data,aes(x =generosity_index_money, y = generosity_index_time, color = svo_angle))+ 
  geom_rug(aes(color=svo_angle))+
  geom_point(position = position_jitter())+
  coord_equal()+
  geom_density2d()+
  geom_abline(slope = 1)+
  scale_x_continuous("Generosity for money") +
  scale_y_continuous("Generosity for time")
ggsave("Figs/svo_angle_given_time_vs_money.png",units = "cm", dpi=300, width = 15,height=10)
```

## doing some models

# H1 t.test

```{r long data}

data.long %>% 
  group_by(condition) %>% 
  summarise(mean = mean(index),
            sd = sd(index))

t.test(data = data.long, index ~ condition, paired = TRUE)
#t.test(data = data.long, raw_score ~ condition, paired = TRUE)


cohensD(index ~ condition,
       data   = data.long,
       method = "paired")
```

We find a medium effect size for the difference of generosity between time and money, d = .45, t(370) = -8.704, p < .001. Participants were more generous when giving time than money.

Controlling for the effect of opportunity costs.

```{r}
m1 <- lmer (index ~ condition + (1|vpn), data = data.long)
summary(m1)


m2.me <- lmer (index ~ condition+req_salary + (1|vpn), data = data.long)
summary(m2.me)
m2.ie <- lmer (index ~ condition*req_salary + (1|vpn), data = data.long)
summary(m2.ie)

anova(m1,m2.me)
sjPlot::tab_model(m1,m2.me,m2.ie)
```

```{r}
data.long$index_predicted <- predict(m1,data.long)


ggplot(data.long,aes(x =condition, y = index_predicted)) + 
  geom_violin()

ggplot(data.long,aes(x =condition, y = index)) + 
  geom_violin()

```


The effect is not influenced by the opportunity cost of time.

# H3 What about personality?

## how many individuals become more prosocial?
```{r}
data %>% summarise(sum(diff_time_money_index > 0)/n())
```


### main effects only

```{r}
m3.hh_alone <- lmer (index ~ condition+honesty_humility_centered  + (1|vpn), data = data.long)

m3 <- lmer (index ~ condition+svo_angle_centered+honesty_humility_centered +miq_centered+req_salary + (1|vpn), data = data.long)
summary(m3.hh_alone)
```

### adding interaction and comparing models


```{r}
m3.svo <- lmer (index ~ condition*svo_angle_centered+honesty_humility_centered +miq_centered+req_salary+ (1|vpn), data = data.long)
m3.hh <- lmer (index ~ condition*honesty_humility_centered +svo_angle_centered+miq_centered+req_salary+ (1|vpn), data = data.long)
m3.miq <- lmer (index ~ condition*miq_centered+honesty_humility_centered +svo_angle_centered+req_salary+(1|vpn), data = data.long)


anova(m3,m3.svo)
anova(m3,m3.hh)
anova(m3,m3.miq)
```

Model results as table
```{r}
sjPlot::tab_model(m3,m3.svo,m3.hh,m3.miq)
```

### why is hh negative?
```{r}
ggplot(data.long %>% filter(!is.na(index)),aes(x = honesty_humility_centered, y = index)) + 
  geom_point()+
  geom_smooth(method = "loess")+
  facet_grid(~condition)
```
Beccause it is not linear.

```{r}
ggplot(data.long %>% filter(!is.na(index)),aes(x = svo_angle_centered, y = index)) + 
  geom_point()+
  geom_smooth(method = "loess")+
  facet_grid(~condition)
```

```{r}
ggplot(data.long %>% filter(!is.na(index)),aes(x = miq_centered, y = index)) + 
  geom_point()+
  geom_smooth(method = "loess")+
  facet_grid(~condition)
```

```{r}
ggplot(data.long %>% filter(!is.na(index)),aes(x = miq_centered, y = diff_time_money_index)) + 
  geom_point()+
  geom_smooth(method = "lm")

ggplot(data.long %>% filter(!is.na(index)),aes(x = svo_angle, y = diff_time_money_index)) + 
  geom_point()+
  geom_smooth(method = "lm")

ggplot(data.long %>% filter(!is.na(index)),aes(x = honesty_humility_score, y = diff_time_money_index)) + 
  geom_point()+
  geom_smooth(method = "lm")
```

```{r}
svo_plot_index <- 
  ggplot(data.long,aes(x =svo_angle, y = index))+ 
  # geom_density_2d_filled(contour_var = "ndensity") +
  geom_density_2d()+
  geom_point(position = position_jitter(), color = "white")+
  scale_x_continuous("SVO angle") +
  scale_y_continuous("Generosity index")+
  geom_smooth(method = "lm", color = "white")+
  facet_grid(~condition)+
  theme(legend.position = "none")

miq_plot_index <-
  ggplot(data.long,aes(x =miq_score, y = index))+ 
  # geom_density_2d_filled(contour_var = "ndensity") +
  geom_density_2d()+
  geom_point(position = position_jitter(), color = "white")+
  scale_x_continuous("Moral identity centrality") +
  geom_smooth(method = "lm", color = "white")+
  scale_y_continuous("Generosity index")+
  facet_grid(~condition) +
  theme(legend.position = "none")

hh_plot_index <- 
  ggplot(data.long,aes(x =honesty_humility_score, y = index))+ 
  # geom_density_2d_filled(contour_var = "ndensity") +
  geom_density_2d()+
  geom_point(position = position_jitter(), color = "white")+
  scale_x_continuous("Honesty-humility") +
  scale_y_continuous("Generosity index")+
  geom_smooth(method = "lm", color = "white")+
  facet_grid(~condition)+
  theme(legend.position = "none")

legend<-get_legend(hh_plot_index + 
                     guides(color = guide_legend(nrow = 1)) +
                     theme(legend.position = "bottom"))

prow <- plot_grid(
  svo_plot_index,
  miq_plot_index,
  hh_plot_index,
  align = 'vh',
  labels = c("A", "B", "C"),
  hjust = -1,
  #vjust = -.3,
  nrow = 3
)

png("Figs/personality_time_vs_money.png",units = "cm", res = 300, width = 30,height=45)
plot_grid(prow,legend, ncol = 1, rel_heights = c(1, .1))
dev.off()
cowplot::plot_grid(svo_plot_index,miq_plot_index,hh_plot_index, nrow = 3)

```


```{r}
data.long$index_predict <- predict(m3,data.long,allow.new.levels = TRUE)
pred.long <- 
  data.long %>% 
  pivot_longer(cols = c("index","index_predict"))

ggplot(pred.long,aes(x =svo_angle_centered, y = value))+ 
  # geom_density_2d_filled(contour_var = "ndensity") +
  geom_density_2d()+
  geom_point(position = position_jitter(), color = "white")+
  scale_x_continuous("SVO angle") +
  scale_y_continuous("Proportion of ressource given to alter")+
  facet_wrap(~condition*name)

```
Our model obviously misses the effect of SVO on giving time.


```{r}
m4 <- lmer (index ~ condition*honesty_humility_score+svo_angle+req_salary + (1|vpn), data = data.long)
summary(m4)
```

```{r}
m5 <- lmer (index ~ condition+miq_score+honesty_humility_score+svo_angle+req_salary + (1|vpn), data = data.long)
summary(m5)
```

## Hurdle model
We assume a two step decision process: give something yes vs. no and a second decision about how much to give (see Engel, 2011).  
Once for the raw count data, which is a truncated poisson distribution.
```{r}
m_hurdle <- glmmTMB(raw ~ condition*svo_angle_centered+(1|vpn), data = data.long,
                    zi=~.,
                    family=truncated_poisson)
summary(m_hurdle)
```

For the generosity index, the main variable of interest

```{r}
m_hurdle.me <- glmmTMB(index_beta_trunc ~ condition+svo_angle_centered + miq_centered+honesty_humility_centered +(1|vpn), data = data.long,
                    zi=~.,
                    family=beta_family())

m_hurdle.svo <- glmmTMB(index_beta_trunc ~ condition*svo_angle_centered +miq_centered + honesty_humility_centered+(1|vpn), data = data.long,
                    zi=~.,
                    family=beta_family())

m_hurdle.miq <- glmmTMB(index_beta_trunc ~ condition*miq_centered+svo_angle_centered + honesty_humility_centered +(1|vpn), data = data.long,
                    zi=~.,
                    family=beta_family())

m_hurdle.hh <- glmmTMB(index_beta_trunc ~ condition* honesty_humility_centered +miq_centered+svo_angle_centered +(1|vpn), data = data.long,
                    zi=~.,
                    family=beta_family())


sjPlot::tab_model(m_hurdle.me,m_hurdle.svo,m_hurdle.miq,m_hurdle.hh,
                  dv.labels = c("Main Effects","Condition*SVO","Condition*Miq","Condition*Honesty-humility"))
anova(m_hurdle.me,m_hurdle.svo)
anova(m_hurdle.me,m_hurdle.miq)
anova(m_hurdle.me,m_hurdle.hh)

summary(m_hurdle.me)

```

# predictions from the main effects hurdle model

-- this seems strange...i have multiple points, because it is not unique for values of the other predictors, so there are multiple values.

```{r}
sims <- simulate(m_hurdle.me, seed = 1, nsim = 100)
simdatlist <- lapply(sims, function(index_beta_trunc){
  cbind(index_beta_trunc, data.long[,c("vpn",'condition','svo_angle_centered','miq_centered','honesty_humility_centered')])
                         })
simdatsums <- lapply(simdatlist, function(x){
  plyr::ddply(x, ~ condition + svo_angle_centered + miq_centered, summarize,
        absence=mean(index_beta_trunc==0),
        mu=mean(index_beta_trunc))})

ssd <- do.call(rbind, simdatsums)

real <- plyr::ddply(data.long, ~condition + svo_angle_centered + miq_centered, summarize,
            absence=mean(index_beta_trunc==0),
            mu=mean(index_beta_trunc))

ggplot(ssd,aes(x=absence, color=condition))+
  geom_density(adjust=4) +
  geom_point(data=real,aes(x=absence, y=1, color=condition), size=2)+
  facet_wrap(~miq_centered) +
  xlab("Probability that no generosity is observed")+
  ylab(NULL)
```

```{r}
ggplot(ssd,aes(x=mu, color=condition))+
  geom_density(adjust=4)+
  geom_point(data=real,aes(x=mu, y=.5, color=condition), size=2)+
  xlab("Generosity")+ylab(NULL)
```


the regressions separated by hand, important the :

```{r}
m1 <- glmer(non_zero ~ condition*svo_angle_centered+(1|vpn), data = data.long, family = binomial(link = logit))
m2 <- glmmTMB(index ~ condition*svo_angle_centered+(1|vpn), data = data.long %>% filter(non_zero==TRUE), beta_family())

summary(m1)
summary(m2)
```

### Predictions of the hurdle model

```{r hurdle model predictions}
hurdle_predictions <- 
  data.long %>%
  na.omit() %>%
  mutate(Proportion = round(index, 1))

# predictions from both steps of the model
hurdle_predictions$non_zero_pred <- predict(m1, hurdle_predictions, 
                                            allow.new.levels = TRUE, type = "response")
hurdle_predictions$score_predict <- predict(m2, hurdle_predictions, allow.new.levels = TRUE, type = "response")

hurdle_predictions <-
  hurdle_predictions %>%
  mutate(Prediction = round(score_predict, 1),
         Prediction = ifelse(non_zero_pred < 0.5, 0, Prediction))

original <-
  hurdle_predictions %>%
  group_by(condition) %>%
  mutate(Ntotal = n()) %>%
  ungroup() %>%
  
  # computing number of participants who gave out each proportion of resource
  group_by(condition, Proportion, Ntotal) %>%
  dplyr::summarise(Nvp = n()) %>%
  mutate(Pvp = Nvp/Ntotal,
         Lvp = binom::binom.bayes(Nvp, Ntotal)$lower,
         Uvp = binom::binom.bayes(Nvp, Ntotal)$upper) %>%
  rename(Condition = condition)



# compute difference between observed proportions
compare_proportions_model <- brms::brm(Nvp|trials(Ntotal) ~ Condition, 
                                       data = single_proportion,
                                       family = "binomial",
                                       chains=0)

compare_proportions <- function (single_proportion){
  compare_proportions_fit <- update(compare_proportions_model, newdata=single_proportion, refresh=0)
  
  bayesian <- 
    brms::posterior_summary(compare_proportions_fit, probs=c( (1-0.89)/2, 1-(1-0.89)/2), pars="b_Conditiontime") %>%
    data.frame() %>%
    mutate(CI = sprintf('%.2f\n%.2f', exp(Q94.5), exp(Q5.5))) %>%
    select(CI)
  
  
  # same via Fisher exact test
  single_wide <-
    single_proportion %>%
    ungroup() %>%
    select(Condition, Proportion, Ntotal, Nvp) %>%
    pivot_wider(id_cols=c("Proportion", "Ntotal"), names_from=Condition, values_from=Nvp)
  ft <- fisher.test(matrix(c(single_wide$time[1], single_wide$Ntotal[1]-single_wide$time[1], 
                             single_wide$money[1], single_wide$Ntotal[1]-single_wide$money[1]), 
                           ncol=2))
  bayesian %>%
    mutate(FT = sprintf("%.2f\n%.2f", ft$conf.int[2], ft$conf.int[1]), 
           NoOverlap = ft$conf.int[1] > 1 | ft$conf.int[2] < 1)
}

proportion_difference <-
  original %>%
  group_by(Proportion) %>%
  do(compare_proportions(.)) %>%
  mutate(Face = ifelse(NoOverlap, "bold", "plain"))


predicted <-
  hurdle_predictions %>%
  group_by(condition) %>%
  mutate(Ntotal = n()) %>%
  ungroup() %>%
  
  # computing number of participants who gave out each proportion of resource
  group_by(condition, Prediction, Ntotal) %>%
  dplyr::summarise(Nvp = n()) %>%
  mutate(Pvp = Nvp/Ntotal,
         Lvp = binom::binom.bayes(Nvp, Ntotal)$lower,
         Uvp = binom::binom.bayes(Nvp, Ntotal)$upper) %>%
  rename(Condition = condition)
  

condition_plot <-
  ggplot(original, aes(x=Proportion, y=Pvp, color=Condition)) + 
  geom_ribbon(data=predicted, aes(x=Prediction, ymin=Lvp, ymax=Uvp, fill=Condition), color=NA, alpha=0.25) + 
  geom_errorbar(aes(ymin=Lvp, ymax=Uvp)) + 
  geom_point() +
  geom_line(data=predicted, aes(x=Prediction)) + 
  theme(legend.position = c(0.8, 0.7), legend.title = element_blank()) +
  scale_y_continuous(name="Proportion of participants", breaks = seq(0, 0.5, 0.1)) +
  scale_x_continuous(name="Generosity", breaks=seq(0, 1, 0.1),
                     sec.axis = dup_axis(name="Fisher exact test (odds)", labels=proportion_difference$CI)) +
  theme(axis.text.x.top = element_text(face = proportion_difference$Face),
        legend.background=element_blank()) +
  plot_theme()

condition_plot
```



```{r}
count_per_choice <-
  data %>%
  group_by(generosity_index_money, generosity_index_time) %>%
  summarise(Count= n(),
            MedianSVO = median(svo_angle)) %>%
  ungroup() %>%
  mutate(Freq = 100 * Count / sum(Count))


plot_time_vs_money <-
  ggplot(data, aes(x=generosity_index_money, y = generosity_index_time))+ 
  # geom_point(position = position_jitter(height = 0.01, width=0.01), alpha = .4, aes(color=svo_angle, fill=svo_angle))+
  geom_abline(slope = 1, color="white", size=2)+
  geom_density2d(show.legend = FALSE) +
  geom_point(data=count_per_choice, aes(size=Freq)) +
  coord_equal()+
  scale_x_continuous("Generosity for money", breaks=seq(0, 1, 0.1)) +
  scale_y_continuous("Generosity for time", breaks=seq(0, 1, 0.1)) + 
  theme(legend.position = c(0.9, 0.2)) +
  labs(size = "Proportion") +
  guides(size = guide_legend(reverse=TRUE)) +
  plot_theme()

plot_time_vs_money
```


```{r}
generosity_difference <- 
  data %>%
  mutate(Time_Money = round(generosity_index_time-generosity_index_money, 1)) %>%
  drop_na(Time_Money) %>%
  group_by(Time_Money) %>%
  summarise(Count = n()) %>%
  ungroup() %>%
  mutate(Prop = Count / sum(Count),
         LoCI = binom::binom.bayes(Count, sum(Count))$lower,
         HiCI = binom::binom.bayes(Count, sum(Count))$upper)

matched <- 
  generosity_difference %>%
  filter(Time_Money == 0) %>%
  pull()

approximately_matched <-
  generosity_difference %>%
  filter(abs(Time_Money) <= 0.1) %>%
  summarise(TotalProp = sum(Prop)) %>%
  pull(TotalProp)

matched_summary <- sprintf('Matched generosity: %.1f%%, exact: %.1f%%', 100*approximately_matched, 100*matched)
  

approx_matched <- 
  generosity_difference %>%
  filter(abs(Time_Money) <= 0.1)

generosity_diff_plot <-
  ggplot(generosity_difference, aes(x=Time_Money, y=Prop)) +
  # geom_ribbon(aes(ymin=LoCI, ymax=HiCI), fill="black", alpha=0.25)+
  geom_ribbon(data=approx_matched, aes(ymin=0, ymax=Prop), fill="black", alpha=0.25)+
  geom_line() + 
  # geom_vline(xintercept = c(-0.1, 0.1)) +
  labs(subtitle = matched_summary) +
  scale_x_continuous(name="Generosity Time - Generosity Money") +
  scale_y_continuous(name="Proportion of participants") +
  plot_theme()

generosity_diff_plot
```

```{r warning=FALSE}
all_plots <- (condition_plot / generosity_diff_plot) | plot_time_vs_money

all_plots + plot_annotation(tag_levels = 'A') + plot_layout(widths = c(2, 2))

ggsave("Figs/Results.png", units = "cm", dpi=300, width=18, height=12)
```



## Other plotting


```{r}
data.long$non_zero_pred <- 
  predict(m1,data.long,allow.new.levels = TRUE,type = "response")
pred.long <-
  data.long %>% 
  mutate(non_zero_pred = case_when(non_zero_pred>.5~1,
                                   TRUE ~ 0)) %>% 
  pivot_longer(cols = c("non_zero","non_zero_pred"))


ggplot(pred.long ,aes(x =svo_angle, y = value, color = condition))+ 
  stat_smooth(method="glm",
                method.args = list(family="binomial"), se=TRUE,
                fullrange=TRUE)+
  geom_point( color = "black")+
  scale_x_continuous("SVO angle") +
  scale_y_discrete("Proportion of ressource given to alter")+
  facet_wrap(~name)
```


```{r}
data.long$score_predict <- predict(m2,data.long,allow.new.levels = TRUE,type = "response")
pred.long<- 
  data.long %>% 
  pivot_longer(cols = c("index","score_predict"))

ggplot(pred.long,aes(x =svo_angle, y = value))+ 
  # geom_density_2d_filled(contour_var = "ndensity") +
  geom_density_2d(color="white")+
  # geom_point(position = position_jitter(), color = "white")+
  scale_x_continuous("SVO angle") +
  scale_y_continuous("Proportion of ressource given to alter")+
  facet_grid(condition~name)
```


```{r}
# getting it all in one prediciton
coef(m1)
```



## Sasha: binomial mixture of individualists/prosocialists/altruists.


### Binomial mixture model with exponential penalty for non-canonical give-aways

```{stan output.var="exponential_3mixture_model"}
data{
  int<lower=1> rowN;

  real<lower=0,upper=1> Proportion[rowN];
  int<lower=0, upper=1> condition[rowN];
  int<lower=0, upper=1> svo_type[rowN];
  int<lower=0> Nvp[rowN];
  int<lower=0> Ntotal[rowN];
}
parameters{
  real indF_a;
  real indF_bC;
  real indF_bT;

  real proF_a;
  real proF_bC;
  real proF_bT;

  real<lower=0> indKappa;
  real<lower=0> proKappa;
  real<lower=0> altKappa;
}
transformed parameters{
  real<lower=0,upper=1> fProAlt;
  real<lower=0,upper=1> indF[rowN];
  real<lower=0,upper=1> proF[rowN];
  real<lower=0,upper=1> altF[rowN];
  
  real thetaRaw[rowN];
  real thetaSum[2, 2];
  real<lower=0,upper=1> theta[rowN];
  
  // initialize counters to 0
  for(iC in 1:2){
    for(iT in 1:2){
      thetaSum[iC, iT] = 0;
    }
  }
  
  for(i in 1:rowN){
    // computing fraction of individualists of everybody
    indF[i] = inv_logit(indF_a + indF_bC * condition[i] + indF_bT * svo_type[i]);
    
    // computing fraction of prosocial out of prosocial + altruists group
    fProAlt = inv_logit(proF_a + proF_bC * condition[i] + proF_bT * svo_type[i]);
    proF[i] = (1 - indF[i]) * fProAlt;
    altF[i] = (1 - indF[i]) * (1 - fProAlt);
    
    thetaRaw[i] = indF[i] * exp(-fabs(0-Proportion[i]) / indKappa) + 
                  proF[i] * exp(-fabs(0.5-Proportion[i]) / proKappa) + 
                  altF[i] * exp(-fabs(1-Proportion[i]) / altKappa);

    thetaSum[condition[i]+1, svo_type[i]+1] += thetaRaw[i];
  }
  
  // normalizing theta to sum up to 1 for each condition and resource type
  for(i in 1:rowN){
    theta[i] = thetaRaw[i] / thetaSum[condition[i]+1, svo_type[i]+1];
  }
}
model{
  for (i in 1:rowN){
    target += binomial_lpmf(Nvp[i] | Ntotal[i], theta[i]);
  }

  target += gamma_lpdf(indKappa | 2, 0.05);
  target += gamma_lpdf(proKappa | 2, 0.05);
  target += gamma_lpdf(altKappa | 2, 0.05);

  target += normal_lpdf(indF_a | 0, 10);
  target += normal_lpdf(indF_bC | 0, 1);
  target += normal_lpdf(indF_bT | 0, 1);

  target += normal_lpdf(proF_a | 0, 10);
  target += normal_lpdf(proF_bC | 0, 1);
  target += normal_lpdf(proF_bT | 0, 1);
}
generated quantities {
  vector[rowN] log_lik;
  for (i in 1:rowN) {
    log_lik[i] = binomial_lpmf(Nvp[i] | Ntotal[i], theta[i]);
  }
}
```

### Preprocessing data for Stan

```{r data-for-stan}
tibble4stan <-
  data.long %>%
  select(vpn, svo_angle, condition, raw, index) %>%
  
  # rename index into proportion (makes more sense to me)
  rename(Ngiven = raw,
         Proportion = index) %>%
  
  # NO NA!
  na.omit() %>%
  
  # split participants into individualists and prosocialists
  mutate(svo_type = ifelse(svo_angle < 22.45, "individualist", "prosocial")) %>%
  
  # make sure proportion is within 0..1
  mutate(Proportion = ifelse(Proportion<0, 0, Proportion),
         Proportion = ifelse(Proportion>1, 1, Proportion)) %>%

  
  # turn condition into binary column IsTime
  mutate(IsTime = as.integer(condition == "time")) %>%
  
  # computing total number of participants 
  group_by(condition, svo_type) %>%
  mutate(Ntotal = n()) %>%
  ungroup() %>%
  
  # computing number of participants who gave out each proportion of resource
  group_by(condition, svo_type, Proportion, Ntotal) %>%
  dplyr::summarise(Nvp = n()) %>%
  mutate(Pvp = Nvp/Ntotal,
         Lvp = binom::binom.bayes(Nvp, Ntotal)$lower,
         Uvp = binom::binom.bayes(Nvp, Ntotal)$upper) %>%
  ungroup() %>%
  mutate(IsTime = as.integer(condition == "time"),
         IsPro = as.integer(svo_type == "prosocial"))

# list with data for STAN model
data4stan <- list(
  rowN = nrow(tibble4stan),
  Proportion = tibble4stan$Proportion,
  condition = tibble4stan$IsTime,
  svo_type = tibble4stan$IsPro,
  Nvp = tibble4stan$Nvp,
  Ntotal = tibble4stan$Ntotal
)
```

### Fitting stan model
```{r}
exp_3mix_fit <- rstan::sampling(exponential_3mixture_model, data=data4stan, chains=4, cores=future::availableCores())
rstan::check_hmc_diagnostics(exp_3mix_fit)
```

### Plotting predictions

Real vs. predicted proportions
```{r}
samples <- rstan::extract(exp_3mix_fit)

# lower 97% CI for predicted proportion
tibble4stan$thetaL <- 
  samples$theta %>% data.frame() %>%
  summarize_all(quantile, probs=(1-0.97) / 2) %>%
  as.matrix() %>%
  t()

# upper 97% CI for predicted proportion
tibble4stan$thetaU <- 
  samples$theta %>% data.frame() %>%
  summarize_all(quantile, probs=1-(1-0.97) / 2) %>%
  as.matrix() %>%
  t()

# mean predicted proportion
tibble4stan$theta <- colMeans(samples$theta)

# original versus predicted data
ggplot(tibble4stan, aes(x=Proportion, y=Pvp, color=condition, fill=condition)) + 
  geom_ribbon(aes(ymin=thetaL, ymax=thetaU), alpha=0.25, color=NA) +
  geom_line(aes(y=theta)) +
  geom_errorbar(aes( ymin=Lvp, ymax=Uvp)) +
  geom_point(shape=21, fill="white") + 
  facet_grid(.~svo_type) +
  theme(legend.position = c(0.9, 0.8))
```
Estimated proportions of each type for each resource
```{r}
compute_proportions <- function(samples, isTime, isPro){
  indF <- boot::inv.logit(samples$indF_a + samples$indF_bC * isTime + samples$indF_bT * isPro)
  fProAlt <- boot::inv.logit(samples$proF_a + samples$proF_bC * isTime + samples$proF_bT * isPro);

  bind_rows(
    tibble(GiveAway = "0%",
           Fraction= indF),
    tibble(GiveAway = "50%",
           Fraction= (1 - indF) * fProAlt),
    tibble(GiveAway = "100%",
           Fraction= (1 - indF) * (1 - fProAlt)))
}

giveaway_type_proportions <- 
  expand.grid(Resource = c("Money", "Time"), SVO_type = c("individualist", "prosocial")) %>%
  group_by(Resource, SVO_type) %>%
  do(compute_proportions(samples, as.integer(.$Resource[1] == "Time"), as.integer(.$SVO_type[1]=="prosocial"))) %>%
  ungroup() %>%
  mutate(GiveAway = factor(GiveAway, levels = c("0%", "50%", "100%")))


ggplot(giveaway_type_proportions, aes(x=interaction(SVO_type, GiveAway), y=Fraction, color=SVO_type)) + 
  geom_violin() +
  facet_grid(.~Resource) + 
  # theme(legend.position = c(0.9, 0.8)) +
  ylab("Proportion of participants in a give-away group") 


giveaway_type_proportions_avg <-
  giveaway_type_proportions %>%
  group_by(SVO_type, Resource, GiveAway) %>%
  summarise(Favg = mean(Fraction),
            Flo = quantile(Fraction, (1-0.97)/2),
            Fhi = quantile(Fraction, 1-(1-0.97)/2))
  

ggplot(giveaway_type_proportions_avg, aes(x=GiveAway, y=Favg, color=SVO_type, fill=SVO_type, group=SVO_type)) + 
  geom_ribbon(aes(ymin=Flo, ymax=Fhi), alpha=0.5, color=NA) +
  geom_line() +
  facet_grid(.~Resource) + 
  # theme(legend.position = c(0.9, 0.8)) +
  ylab("Proportion of participants in a give-away group") 

```
Posterior plots for each proportion

```{r}
fraction_posterior <-
  as.tibble(samples[c('indF_bC', 'indF_bT', 'proF_bC', 'proF_bT')]) %>%
  pivot_longer(cols=everything(), names_to="Beta", values_to = "Value") %>%
  mutate(Type = ifelse(stringr::str_starts(Beta, "ind"), "Change in 0% vs. 50/100%", "Change in 50% vs 100%"),
         Term = ifelse(stringr::str_ends(Beta, "C"), "Time (Money baseline)", "prosocial (individualist baseline)"))

fraction_below_zero <-
  fraction_posterior %>%
  group_by(Type, Term) %>%
  summarise(Below = mean(Value<0)) %>%
  mutate(Label = ifelse(Below > 0.5, 
                        glue::glue("β<0: {round(100 * Below, 1)}"),
                        glue::glue("β>0: {round(100 * (1-Below), 1)}")))

ggplot(fraction_posterior, aes(x=Value)) + 
  geom_histogram(, bins=100) + 
  geom_text(data=fraction_below_zero, aes(label=Label, x=Inf, y=Inf), hjust="right", vjust="top") +
  geom_vline(xintercept = 0, color="white") +
  facet_grid(Type ~ Term)
```

# H2 Affect differs depending on who you work for

```{r}
ggplot(data.long, aes(y = panas_positive, x = index, color = condition)) + geom_smooth()
```


```{r}
beanplot(panas_positive~worked_for, data = data %>% filter(worked_for!="non"))
beanplot(panas_negative~worked_for, data = data %>% filter(worked_for!="non"))

```
```{r}
data %>% filter(worked_for!="non") %>% 
  group_by(worked_for) %>% 
  summarize(mean = mean(panas_positive,na.rm = TRUE),
            sd = sd(panas_positive,na.rm = TRUE))

data %>% filter(worked_for!="non") %>% 
  drop_na(panas_positive) %>% 
  group_by(worked_for) %>% summarise(n())
```


```{r}

t.test(panas_positive~worked_for, data = data %>% filter(worked_for!="non"))
cohensD(panas_positive~worked_for,
       data   =data %>% filter(worked_for!="non"))

```
```{r}

data %>% filter(worked_for!="non") %>% 
  group_by(worked_for) %>% 
  select(panas_positive,generosity_index_time,generosity_index_money,svo_angle,miq_score,honesty_humility_score) %>% correlation()

typeof(data$panas_positive)
```

