library(foreign)
library(ggplot2)
library(Rmisc)
library(pastecs)
library(scales)
library(lsr)
library(Hmisc)
#library(Rcmdr)
library(extrafont)
library(reshape)
library (lawstat)
library(tibble)
library(pwr)
library(lme4)
library(lmerTest)
library(sme)

##############################################################################
# Code for paper "How selfish is a thirsty man? A pilot study on comparing sharing behavior with primary and secondary rewards" (Kause, Vitouch & Gl??ck, 2017).
# Study conducted 2013 in the lab of the Group for Adaptive Behavior and Cognition, Max Planck Institute for Human Development, Berlin
##############################################################################

# The dataset excludes 1 person who participated by accident twice and 2 who had technical problems with the online questionnaire.

# Path to file needs to be replaced and adapted to user of code and data
Results_DG <- read.csv("~/Documents/EcGames/EcGames_MPIB/Submissions/Summer_2017/Leeds University Data Repository/WaterMoney_Data.csv", sep=";", dec=",", header= T)
d<-Results_DG[1:82,1:36]
d[d == "-77" ] = NA
d=as.data.frame(d);

## Recode variables from absolute into percentage shares
d$prop_dictator_water_p<-(d$prop_dictator_water/d$water_overall) # kept by proposer
d$prop_dictator_money_p<-(d$prop_dictator_money/d$amount_overall) 
d$respo_water_p<-(d$respo_water/d$water_overall) # received from proposer
d$respo_money_p<-(d$respo_money/d$amount_overall)
d$exp_respo_water_p<-(d$exp_respo_water/d$water_overall)
d$exp_respo_money_p<-(d$exp_respo_money/d$amount_overall)

# Condition
## recode into verbal labels
for (i in 1:nrow(d))
{if (d$Condition[i]=="2")  {d$condition_verb[i]<-"baseline"}}

for (i in 1:nrow(d))
{if (d$Condition[i]=="1")  {d$condition_verb[i]<-"earned"}}

# Offers in %
d$respo_water_pe<-d$respo_water_p*100
d$respo_money_pe<-d$respo_money_p*100

write.csv(d, file="WaterMoney_Data_processed.csv")

## Split dataset according to condition
d_earn<-subset(d, d$Condition=="1") # earned
d_base<-subset(d, d$Condition=="2") # windfall

## Split dataset according to roles
 d_proposer = subset(d, role =="1")
 d_responder = subset(d, role =="2")

## Description of sample
### age; both overall & per condition
stat.desc(d$age);
stat.desc(d_base$age);
stat.desc(d_earn$age);

## Demographics and sharing behavior; separated by roles

### Linear model, money offers by proposers
lm_proposers.m = lm(prop_dictator_money_p~Condition+Reihenfolge+sex+age+education+sport.week+smoker+drink.time+puls_post+watt_18, data = d_proposer)
## proposers: water offers
lm_proposers.w = lm(prop_dictator_water_p~Condition+Reihenfolge+sex+age+education+sport.week+smoker+drink.time+puls_post+watt_18, data = d_proposer)

# Linear model, expectations by proposers
## responders: money offers
lm_responders.m = lm(exp_respo_money_p~Condition+Reihenfolge+sex+age+education+sport.week+smoker+drink.time+puls_post+watt_18, data = d_responder)
# responders: water offers
lm_responders.w = lm(exp_respo_water_p~Condition+Reihenfolge+sex+age+education+sport.week+smoker+drink.time+puls_post+watt_18, data = d_responder)

## Rise of thirst during bike session
options(digits=2)
thirst<-stat.desc(d[,22:26]);

## How do participants share both goods and what do they expect?
### Overall descriptive statistics
desc_overall<-stat.desc(d[,37:42]) # descriptives

### Distribution of offers
table(d$respo_water_p); 
table(d$respo_money_p);

## What do receivers expect, depending on condition?
summarySE(d,measurevar="exp_respo_money_p", groupvars="Condition", na.rm=T)
summarySE(d,measurevar="exp_respo_water_p", groupvars="Condition", na.rm=T)

## Overall T-tests: Height of offers and expectations
### Levene: Equality of variances
#### Offers
boem<-c(d$respo_money_p, d$respo_water_p)
boeme<-as.factor(c(rep(1, length(d$respo_money_p)), rep(2, length(d$respo_water_p))))
levene.test(boem,boeme);

#### Expectations
eboem<-c(d$exp_respo_money_p, d$exp_respo_water_p)
eboeme<-as.factor(c(rep(1, length(d$respo_money_p)), rep(2, length(d$respo_water_p))))
levene.test(eboem,eboeme); # for similarity of distributions

## T-Test (student), pooled
options(scipen=999)
d_resp=subset(d, role=="2") # only responders
t.test(d_resp$respo_money_p, d_resp$respo_water_p, var.equal=T, paired=T); # Actual offers
# t.test(d_resp$exp_respo_money_p, d_resp$exp_respo_water_p, var.equal=T, paired=T); # Expectations, don't use t#test as levene indicates inequality of variances
wilcox.test(d_resp$exp_respo_money_p, d_resp$exp_respo_water_p, paired=TRUE, alternative="less") # use Wilcoxon Test instead

### Cohen's d for t-test, pooled
cohensD(d_resp$respo_water_p, d_resp$respo_money_p, method="paired"); # Actual offers
cohensD(d_resp$exp_respo_water_p, d_resp$exp_respo_money_p, method="paired"); # Expectations

# PROPOSERS
### Long dataset for linear model analysis
da=subset(d, role=="2")
datalong=melt(da, id=c("Participant.nr.", "Condition","thirst_23", "Reihenfolge","sex","age","education","sport.week","smoker","drink.time","puls_post","watt_18"), measure=c("respo_water_p", "respo_money_p"), variable.name="good") # melts data into format that allows comparing values of participants across goods

### ANOVA of good shared x condition (earned versus baseline); actual offers
x=aov(value~Condition*variable+Error(Participant.nr./variable), data=datalong) # ANOVA with value as DV, Condition as between and variable(=shared good) as within variable

### Linear model of offers by proposers, depending on good shared
lm_proposers.simple = lmer(value~Condition* variable+(1 | Participant.nr.), data = datalong)
lm_proposers.simple2 = lmer(value~Condition+(1 | Participant.nr.), data = datalong)
lm_proposers.long = lmer(value~Condition*variable+Reihenfolge+sex+age+education+sport.week+smoker+drink.time+puls_post+thirst_23+(1 | Participant.nr.), data = datalong)
lm_proposers.long.gender = lmer(value~Condition*variable+Reihenfolge+variable*sex+age+education+sport.week+smoker+drink.time+puls_post+thirst_23+(1 | Participant.nr.), data = datalong)

AICc(lm_proposers.long)

# RESPONDERS
datalongex=melt(d_resp, id=c("Participant.nr.", "Condition","thirst_23", "Reihenfolge","sex","age","education","sport.week","smoker","drink.time","puls_post","watt_18"), measure=c("exp_respo_water_p", "exp_respo_money_p"), variable.name="good") # melts data into format that allows comparing values of participants across goods

### ANOVA of good shared x condition (earned versus baseline); expectations
ex=aov(value~Condition*variable+Error(Participant.nr./variable), data=datalongex) # ANOVA with value as DV, Condition as between and variable(=shared good) as within variable

## Linear model predicting expecation by responders
lm_responders.simple = lmer(value~Condition* variable+(1 | Participant.nr.), data = datalongex) # Simple model
lm_responders.long = lmer(value~Condition*variable+Reihenfolge+sex+age+education+sport.week+smoker+drink.time+puls_post+thirst_23+(1 | Participant.nr.), data = datalongex) # Full model
lm_responders.long.gender = lmer(value~Condition*variable+Reihenfolge+variable*sex+age+education+sport.week+smoker+drink.time+puls_post+thirst_23+(1 | Participant.nr.), data = datalongex) # Full model

## Do proposers make different offers depending on the condition (earning versus baseline)?
### T-test, Welch for testing differences between conditions 
mt_overall<-t.test(d_earn$respo_money_p, d_base$respo_money_p,var.equal=T) # money: earned vs. baseline
wt_overall<-t.test(d_earn$respo_water_p, d_base$respo_water_p,var.equal=T) # water: earned vs. baseline

# Cohen's d - mt_overall
a_var<-d_base$respo_money_p
b_var<-d_earn$respo_money_p
cohensD(a_var,b_var); 

## Cohen's d - for wt_overall 
a_var<-d_base$respo_water_p
b_var<-d_earn$respo_water_p
cohensD(a_var,b_var);

## Visual inspection of results
### WATER: Do water offers differ between experimental conditions?
condition_w<-group.CI(respo_water_p~Condition, d, ci=0.95)
condition_w$Condition <- factor(condition_w$Condition)
ggplot(condition_w, aes(x=Condition, colour=Condition, y=respo_water_p.mean, title='Water: offer by condition'))+
  geom_errorbar(aes(ymin=respo_water_p.lower, ymax=respo_water_p.upper), width=0.1)+
  geom_point(size=5) + scale_colour_brewer(palette="Set1", type="qual", name = " ") +
  scale_y_continuous(limits = c(0, 1))+
  ylab("Percentage of water offered") +
  theme_bw(base_size=20)

### MONEY: Do offers differ between experimental conditions?
condition_m<-group.CI(respo_money_p~Condition, d, ci=0.95)
condition_m$Condition <- factor(condition_m$Condition)
ggplot(condition_m, aes(x=Condition, colour=Condition, y=respo_money_p.mean, title='Money: offer by condition'))+
  geom_errorbar(aes(ymin=respo_money_p.lower, ymax=respo_money_p.upper), width=0.1)+
  geom_point(size=5) + scale_colour_brewer(palette="Set1", type="qual", name = " ") +
  scale_y_continuous(limits = c(0, 1))+
  ylab("Percentage of money offered") +
  theme_bw(base_size=20)

######################################
## Do participants behave differently depending on amount shared?

### ANOVA, on water offer by amount shared
#### Levene
?LeveneTest
anova.water1=aov(respo_water_p~water_overall, data = d)
anova.water2=oneway.test(respo_water_p~water_overall, data = d) # Welch's F, adjusted for unequal variances
summary.lm(anova.water1) # provides R2 as effect size (see Field, p. 454)
pwr.anova.test(k=3,n=41,f=0.03977928,sig.level=.82) # Power. Input: Effect size F was calculated with G*power (omnibus ANOVA) = 0.03977928; based on means in summarySE(d, measurevar="respo_water_p", groupvar="water_overall", na.rm=T)

### ANOVA, on money offer by amount shared 
anova.money1=aov(respo_money_p~amount_overall, data = d)
anova.money2=oneway.test(respo_money_p~amount_overall, data = d) # Welch's F, adjusted for unequal variances
summary.lm(anova.money1)
# Effect size F, calculated with G*power (omnibus ANOVA) = 0.03977928; based on means in summarySE(d, measurevar="respo_money_p", groupvar="water_overall", na.rm=T)
pwr.anova.test(k=3,n=41,f=3.03,sig.level=.01)

## Visual inspection
# Water offers, separated by amount
water<-group.CI(respo_water_p~water_overall, d, ci=0.95)
water$water_overall <- factor(water$water_overall)
ggplot(water, aes(x=water_overall, colour=water_overall, y=respo_water_p.mean, title="Water: Offer by amount"))+
  geom_errorbar(aes(ymin=respo_water_p.lower, ymax=respo_water_p.upper), width=0.1)+
  geom_point(size=5) + scale_colour_brewer(palette="Set2", type="qual", name = " ") +
  scale_y_continuous(limits = c(0, 1))+
  theme_bw(base_size=20) +
  xlab("Amount shared")+
  ylab("Mean (+CI)") 

# Money offers, separated by amount
money<-group.CI(respo_money_p~amount_overall, d, ci=0.95)
money$amount_overall <- factor(money$amount_overall)
ggplot(money, aes(x=amount_overall, colour=amount_overall, y=respo_money_p.mean, title='Money: Offer by amount'))+
  geom_errorbar(aes(ymin=respo_money_p.lower, ymax=respo_money_p.upper), width=0.1)+
  geom_point(size=5) + scale_colour_brewer(palette="Set2", type="qual", name = " ") +
  scale_y_continuous(limits = c(0, 1))+
  theme_bw(base_size=20) +
  xlab("Amount shared")+
  ylab("Mean (+CI)")

### Did order effects occur in water/ money offer?
summarySE(d_proposer, measurevar="prop_dictator_money_p", groupvar="Reihenfolge", na.rm=T)
t.test(d_proposer$prop_dictator_money_p[which(d$Reihenfolge==2)], d_proposer$prop_dictator_money_p[which(d$Reihenfolge==1)], var.equal=T); # Money
cohensD(d_proposer$prop_dictator_money_p[which(d$Reihenfolge==2)], d_proposer$prop_dictator_money_p[which(d$Reihenfolge==1)])

summarySE(d_proposer, measurevar="prop_dictator_water_p", groupvar="Reihenfolge", na.rm=T)
t.test(d_proposer$prop_dictator_water_p[which(d$Reihenfolge==2)], d_proposer$prop_dictator_water_p[which(d$Reihenfolge==1)], var.equal=T); # Water
cohensD(d_proposer$prop_dictator_water_p[which(d$Reihenfolge==2)], d_proposer$prop_dictator_water_p[which(d$Reihenfolge==1)])

####### Additional analysis, not reported in paper
########
# factorize: Condition
d$Condition_f<-factor(d$Condition)
levels(d$Condition_f) <- c("earned","windfall")
########

# Comparison of water windfall (MPIB) with money windfall DG from Forsythe et al.
### Values of DG Forsythe et al. (1994), taken from Appendix B
DGFors10=c(5,3,1,0,2,3,3,0,2,3,5,5,3,1,5,1,3,3,0,0,0,1,5,2) # raw values
DGFors10p=(DGFors10/10)*100 # recalculation into percent of overall amount which was 10 dollars
d_resp_wind=subset(d_resp, Condition=="2")
t.test(DGFors10p, (d_resp_wind$respo_water_p*100))
cohensD(DGFors10p, (d_resp_wind$respo_water_p*100))

######################################
## Graphs like in paper
# MPIB study with 2x2 design: Overall distribution of water offers
respo_water_g <- qplot(respo_water_pe, data=d, geom="histogram", binwidth=2, xlim=c(0,100), ylim=c(0,20), main="Frequency of Water Offers (in %)",facets=.~ d$Condition_f,  na.rm=T) 
respo_water_g <- respo_water_g + theme_bw(base_size=15) + labs(x="% offered", y="frequency", fill="") +theme(axis.text=element_text(size=9)) +theme(axis.title=element_text(size=10))+theme(plot.title=element_text(size=12))#+theme(text=element_text(family="Helvetica57-Condensed"))
respo_water_g

## Overall distribution of monetary offers
respo_money_g <- qplot(respo_money_pe, data=d, geom="histogram", binwidth=2, xlim=c(0,100),ylim=c(0,11), main="Frequency of Monetary Offers (in %)", facets=. ~ Condition_f, na.rm=T)
respo_money_g <- respo_money_g + theme_bw(base_size=15) + labs(x="% offered", y="frequency", fill="") +theme(axis.text=element_text(size=9))  +theme(axis.title=element_text(size=10))+theme(plot.title=element_text(size=12))#+theme(text=element_text(family="Helvetica57-Condensed"))
respo_money_g