Data analysis Code

---

Code

expe <- read_excel("DATA - 10683_2007_9171_MOESM2_ESM.xls")
colnames(expe) <- expe[4,]
expe <- as.data.frame(expe[-c(1:4),])
expe[is.na(expe)==FALSE] <- as.numeric(expe[is.na(expe)==FALSE])
#Original graph
#G1
forgr1.df <- select(expe, Period, Group, Treatment, mycont)
forgr1.df$Period <- as.double(forgr1.df$Period)
forgr1.df$mycont <- as.double(forgr1.df$mycont)
forgr1.df <- filter(forgr1.df, Period>5)
forgr1.df$Group <- as.character(forgr1.df$Group)
forgr1_res.df <- as.data.frame(aggregate(forgr1.df[,4], list(forgr1.df$Treatment), mean))
forgr1_res.df$x <- round(forgr1_res.df$x,digits = 2)
forgr1_res.df$Group.1 <- as.character(forgr1_res.df$Group.1)
EEDAFigure1 <- ggplot(forgr1_res.df, aes(x=Group.1,y=x,fill=Group.1),color="black")+
  geom_bar(stat="identity")+
  scale_fill_brewer(palette="Set2")+
  geom_text(label=forgr1_res.df$x,vjust=-.25)+
  labs(title = "Average contribution by treatment (data from periods 6–10)",
       x="Treatment",
       y="Average contribution",
       fill=NULL)+
  theme_light()+
  theme(legend.background = element_rect(linetype="solid",color="Black"))
#ggsave("EEDA-Figure1.png",EEDAFigure1)
#G2
forgr2.df <- select(expe, Period, Group, Treatment, mycont)
forgr2.df$Period <- as.double(forgr2.df$Period)
forgr2.df$mycont <- as.double(forgr2.df$mycont)
forgr2.df$Group <- as.character(forgr2.df$Group)
forgr2_res.df <- data.frame(aggregate(x=forgr2.df$mycont,
                                      by=list(forgr2.df$Period,forgr2.df$Treatment),
                                      FUN=mean))
forgr2_res.df$x <- round(forgr2_res.df$x,digits = 2)
forgr2_res.df$Group.1 <- as.integer(forgr2_res.df$Group.1)
forgr2_res.df$Group.2 <- as.character(forgr2_res.df$Group.2)
EEDAFigure2 <- ggplot(forgr2_res.df,aes(x=Group.1,y=x))+
  scale_x_continuous(breaks=c(1:10))+
  scale_y_continuous(breaks=seq(0,20,2))+
  geom_point(aes(shape=Group.2))+
  geom_line(aes(group =Group.2))+
  scale_shape_manual(values=c(1,0,15,5,8))+
  labs(title = "Average contribution over time across treatments",
       x="Period",
       y="Average contribution",
       shape=NULL)+
  theme_light()+
  theme(legend.position = c(0.07, 0.22),
        legend.background = element_rect(linetype="solid",color="Black"),
        plot.title = element_text(hjust = 0.5))
#ggsave("EEDA-Figure2.png",EEDAFigure2)
#G3
forgr3.df <- select(expe, Period, Group, Treatment, TotalProfit)
forgr3.df$Period <- as.double(forgr3.df$Period)
forgr3.df$TotalProfit <- as.double(forgr3.df$TotalProfit)
forgr3.df$Group <- as.character(forgr3.df$Group)
forgr3_res.df <- data.frame(aggregate(x=forgr3.df$TotalProfit,
                                      by=list(forgr3.df$Period,forgr3.df$Treatment),
                                      FUN=mean))
forgr3_res.df$x <- round(forgr3_res.df$x,digits = 2)
forgr3_res.df$Group.1 <- as.integer(forgr3_res.df$Group.1)
forgr3_res.df$Group.2 <- as.character(forgr3_res.df$Group.2)
Base <- forgr3_res.df$x[forgr3_res.df$Group.2==0]
Base <- c(rep(Base,4))
forgr3_res.df <- forgr3_res.df[-c(1:10),]
forgr3_res.df$rce <- (forgr3_res.df$x-Base-25)/Base
EEDAFigure3 <- ggplot(forgr3_res.df,aes(x=Group.1,y=rce))+
  geom_hline(aes(yintercept = 0),colour="gray28")+
  scale_x_continuous(breaks=c(1:10))+
  geom_point(aes(shape=Group.2))+
  geom_line(aes(group =Group.2))+
  scale_shape_manual(values=c(0,15,5,8))+
  labs(title = "CRE over time across treatments",
       x="Period",
       y="CRE",
       shape=NULL)+
  theme_bw()+
  theme(legend.position = c(0.9, 0.22),
        legend.background = element_rect(linetype="solid",color="Black"),
        plot.title = element_text(hjust = 0.5))
#ggsave("EEDA-Figure3.png",EEDAFigure3)

expe <- read_excel("DATA - 10683_2007_9171_MOESM2_ESM.xls")
colnames(expe) <- expe[4,]
expe <- as.data.frame(expe[-c(1:4),])
expe[is.na(expe)==FALSE] <- as.numeric(expe[is.na(expe)==FALSE])
#New
#Graph
#H1
#Y=delta Contribution (C t - C t-1)
#X=R.Punishment t-1
h1.df <- select(expe, Period, Group, Treatment, Subject, mycont, tot_red)
h1.df$Treatment <- as.integer(h1.df$Treatment)
h1.df <- filter(h1.df,Treatment>0)
h1.df$Period <- as.double(h1.df$Period)
h1.df$mycont <- as.double(h1.df$mycont)
h1.df$tot_red <- as.double(h1.df$tot_red)
h1.df$Group <- as.character(h1.df$Group)
h1.df$Subject <- as.character(h1.df$Subject)
allpopu1 <- unique(h1.df$Subject)
h1.df$delta <- NA
for(i in 1:length(allpopu1)){
  for(j in 2:10){
    this <- h1.df[(h1.df$Period==j & h1.df$Subject==allpopu1[i]),]$mycont
    last <- h1.df[(h1.df$Period==(j-1) & h1.df$Subject==allpopu1[i]),]$mycont
    h1.df[(h1.df$Period==j & h1.df$Subject==allpopu1[i]),]$delta <- (this - last)
  }
}
h1_res.df <- data.frame(lastpunish=h1.df[(h1.df$Period<10),]$tot_red,
                        delta=h1.df[(h1.df$Period>1),]$delta)
h1_res.df$count <- 1
h1_res2.df <- aggregate(count~lastpunish+delta,
                        data=h1_res.df,
                        sum)
EEDANewH1 <-  ggplot(h1_res2.df,aes(x=lastpunish,y=delta))+
  geom_point(aes(size=log(count),color=log(count)))+
  scale_color_viridis_c()+
  labs(title = "Last Punishment v.s. ΔCont",
       x="Last Punishment",
       y="ΔCont",
       size="ln(Count)",
       color="ln(Count)")+
  theme_bw()+
  theme(legend.background = element_rect(linetype="solid",color="Black"),
        plot.title = element_text(hjust = 0.5))
#ggsave("EEDA-NewH1.png",EEDANewH1)
#H2
#Y=Punishment t
#X=Ci-Cj(j!=i) t
#Xc=Pool t
h2.df <- select(expe,Group,Subject,Period,Treatment,mycont,totcont,`p[1]`,`p[2]`,`p[3]`,`p[4]`)
h2.df$Treatment <- as.integer(h2.df$Treatment)
h2.df <- filter(h2.df,Treatment>0)
h2.df$Group <- as.character(h2.df$Group)
h2.df$Subject <- as.integer(h2.df$Subject)
h2.df$Period <- as.double(h2.df$Period)
h2.df$mycont <- as.double(h2.df$mycont)
h2.df$totcont <- as.double(h2.df$totcont)
h2.df$`p[1]` <- as.double(h2.df$`p[1]`)
h2.df$`p[2]` <- as.double(h2.df$`p[2]`)
h2.df$`p[3]` <- as.double(h2.df$`p[3]`)
h2.df$`p[4]` <- as.double(h2.df$`p[4]`)
h2.df$d1 <- 0
h2.df$d2 <- 0
h2.df$d3 <- 0
h2.df$d4 <- 0
diffmat = function(x){
  D = matrix(as.numeric(NA), NROW(x), NROW(x))
  for (i in 1:NROW(x)){
    d = x[[i]] - x[-i]
    D[i,-i] = d
  }
  if (!all(is.na(diag(D)))){
    stop("Not all diagonal elements zero")
  }
  diag(D) = 0
  if (!is.null(names(x))) colnames(D) = rownames(D) = names(x)
  return(D)
}
k <- nrow(h2.df)/4
for( w in 1:k){
  j <- w*4
  h2.df[(j-3):j,(11:14)] <- diffmat(h2.df[(j-3):j,5])
}
h2.df$tp1 <- h2.df$Treatment*h2.df$`p[1]`
h2.df$tp2 <- h2.df$Treatment*h2.df$`p[2]`
h2.df$tp3 <- h2.df$Treatment*h2.df$`p[3]`
h2.df$tp4 <- h2.df$Treatment*h2.df$`p[4]`
h2_res.df <- c(NA,NA,NA,NA)
for(i in 1:nrow(h2.df)){
  place <- h2.df$Subject[i] %% 10
  for(j in 1:4){
    if(j!=place){
      p <- h2.df[i,(j+6)]
      tp <- h2.df[i,(j+14)]
      d <- h2.df[i,(j+10)]
      pool <- h2.df$totcont[i]
      h2_res.df <- rbind(h2_res.df, c(p,tp,d,pool))
    }
  }
}
h2_res.df <- as.data.frame(h2_res.df[-1,])
colnames(h2_res.df) <- c("punish","rpunish","delta","pool")
h2_res.df$count <- 1
h2_res2.df <- aggregate(count ~ punish + delta,
                        data=h2_res.df,
                        sum)
h2_res3.df <- filter(h2_res2.df,delta>=0)
EEDANewH2 <- ggplot(h2_res3.df,aes(x=delta,y=punish))+
  geom_point(aes(size=log(count),color=log(count)))+
  scale_color_viridis_c()+
  labs(title = "MORE Cont v.s. Punish",
       x="More Contribution",
       y="Punish",
       size="ln(Count)",
       color="ln(Count)",
       caption = "(Only when dif(mycont) is non-negative)")+
  theme_bw()+
  theme(legend.background = element_rect(linetype="solid",color="Black"),
        plot.title = element_text(hjust = 0.5))
#ggsave("EEDA-NewH2.png",EEDANewH2)
EEDANewH2.1 <- ggplot(h2_res2.df,aes(x=delta,y=punish))+
  geom_point(aes(size=log(count),color=log(count)))+
  scale_color_viridis_c()+
  labs(title = "Dif Cont v.s. Punish",
       x="DIfference in Contribution",
       y="Punish",
       size="ln(Count)",
       color="ln(Count)")+
  theme_bw()+
  theme(legend.background = element_rect(linetype="solid",color="Black"),
        plot.title = element_text(hjust = 0.5))
#ggsave("EEDA-NewH2.1.png",EEDANewH2.1)
#Regression
#H1
#Y=delta Contribution (C t - C t-1)
#X=R.Punishment t-1
resh1 <- lm(delta~lastpunish,h1_res.df)
resh1m1 <- data.frame(lastpunish=h1.df[(h1.df$Period<10 & h1.df$Treatment==1),]$tot_red,
                      delta=h1.df[(h1.df$Period>1 & h1.df$Treatment==1),]$delta)
resh1m2 <- data.frame(lastpunish=h1.df[(h1.df$Period<10 & h1.df$Treatment==2),]$tot_red,
                      delta=h1.df[(h1.df$Period>1 & h1.df$Treatment==2),]$delta)
resh1m3 <- data.frame(lastpunish=h1.df[(h1.df$Period<10 & h1.df$Treatment==3),]$tot_red,
                      delta=h1.df[(h1.df$Period>1 & h1.df$Treatment==3),]$delta)
resh1m4 <- data.frame(lastpunish=h1.df[(h1.df$Period<10 & h1.df$Treatment==4),]$tot_red,
                      delta=h1.df[(h1.df$Period>1 & h1.df$Treatment==4),]$delta)
resh11 <- lm(delta~lastpunish,resh1m1)
resh12 <- lm(delta~lastpunish,resh1m2)
resh13 <- lm(delta~lastpunish,resh1m3)
resh14 <- lm(delta~lastpunish,resh1m4)
tl <- c("All","'1'","'2'","'3'","'4'")
#stargazer(resh1,resh11,resh12,resh13,resh14, title="Hypothesis 1",column.labels=tl, align=TRUE, digits=2, out="EET1.htm")
#H2
#Y=Punishment t
#X=Ci-Cj(j!=i) t
#Xc=Pool t
h2_ols.df <- c(NA,NA,NA,NA,NA)
for(i in 1:nrow(h2.df)){
  place <- h2.df$Subject[i] %% 10
  for(j in 1:4){
    if(j!=place){
      tr <- h2.df$Treatment[i]
      p <- h2.df[i,(j+6)]
      tp <- h2.df[i,(j+14)]
      d <- h2.df[i,(j+10)]
      pool <- h2.df$totcont[i]
      h2_ols.df <- rbind(h2_ols.df, c(tr,p,tp,d,pool))
    }
  }
}
h2_ols.df <- as.data.frame(h2_ols.df)
h2_ols.df <- h2_ols.df[-1,]
h2_ols.df$Negative <- 0
colnames(h2_ols.df) <- c("Treatment","Punish","Rpunish","Delta","Pool","Negative")
h2_ols.df$Negative[h2_ols.df$Delta<0] <- 1
h2_ols.df1 <- filter(h2_ols.df,Treatment==1)
h2_ols.df2 <- filter(h2_ols.df,Treatment==2)
h2_ols.df3 <- filter(h2_ols.df,Treatment==3)
h2_ols.df4 <- filter(h2_ols.df,Treatment==4)
resh2 <- lm(Punish ~ Delta + Negative + Negative*Delta + Pool, h2_ols.df)
resh21<- lm(Punish ~ Delta + Negative + Negative*Delta + Pool, h2_ols.df1)
resh22<- lm(Punish ~ Delta + Negative + Negative*Delta + Pool, h2_ols.df2)
resh23<- lm(Punish ~ Delta + Negative + Negative*Delta + Pool, h2_ols.df3)
resh24<- lm(Punish ~ Delta + Negative + Negative*Delta + Pool, h2_ols.df4)
#stargazer(resh2,resh21,resh22,resh23,resh24, title="Hypothesis 2",column.labels=tl, align=TRUE, digits=2, out="EET2.htm")

Result

EEDAFigure1

EEDAFigure2

EEDAFigure3

Hypothesis 1 - Contribute more after more punishment

EEDANewH1

stargazer(resh1,resh11,resh12,resh13,resh14, title="Hypothesis 1",column.labels=tl, align=TRUE, digits=2,type = 'html',style = "qje",df = FALSE)

Hypothesis 1

	delta
	All	‘1’	‘2’	‘3’	‘4’
	(1)	(2)	(3)	(4)	(5)
lastpunish	0.44***	0.80***	0.36***	0.38***	0.39***
	(0.03)	(0.13)	(0.07)	(0.04)	(0.05)
Constant	-0.80***	-1.69***	-1.03***	-0.45**	-0.15
	(0.17)	(0.45)	(0.35)	(0.22)	(0.23)
N	864	216	216	216	216
R2	0.16	0.16	0.10	0.32	0.22
Adjusted R2	0.15	0.15	0.09	0.32	0.21
Residual Std. Error	4.44	6.12	4.56	2.94	3.18
F Statistic	158.18***	39.54***	23.16***	103.00***	58.89***
Notes:	***Significant at the 1 percent level.
	**Significant at the 5 percent level.
	*Significant at the 10 percent level.

Hypothesis 2 - Punish the free riders

EEDANewH2

EEDANewH2.1

stargazer(resh2,resh21,resh22,resh23,resh24, title="Hypothesis 2",column.labels=tl, align=TRUE, digits=2,type = 'html',style = "qje",df = FALSE)

Hypothesis 2

	Punish
	All	‘1’	‘2’	‘3’	‘4’
	(1)	(2)	(3)	(4)	(5)
Delta	0.14***	0.15***	0.10***	0.17***	0.12***
	(0.005)	(0.01)	(0.01)	(0.01)	(0.004)
Negative	0.10*	0.13	0.21*	-0.004	0.06
	(0.05)	(0.14)	(0.12)	(0.08)	(0.05)
Pool	0.003***	0.01***	0.004**	-0.003	0.003***
	(0.001)	(0.002)	(0.002)	(0.002)	(0.001)
Delta:Negative	-0.14***	-0.15***	-0.11***	-0.16***	-0.12***
	(0.01)	(0.02)	(0.01)	(0.02)	(0.01)
Constant	-0.15***	-0.25***	-0.16	0.24*	-0.19***
	(0.05)	(0.09)	(0.12)	(0.14)	(0.07)
N	2,880	720	720	720	720
R2	0.25	0.25	0.15	0.34	0.61
Adjusted R2	0.25	0.25	0.15	0.34	0.61
Residual Std. Error	0.89	1.24	0.97	0.72	0.32
F Statistic	237.74***	59.56***	31.64***	92.68***	284.74***
Notes:	***Significant at the 1 percent level.
	**Significant at the 5 percent level.
	*Significant at the 10 percent level.

---