##### Data analysis for eyes and littering paper #####
# Set working directory 
#setwd("H:/Human Research/Robinson Greenlees Eyes paper/Analysis")
#setwd("D:/Dropbox/Robinson Greenlees Eyes paper/Analysis")
setwd("C:/Users/Daniel/Dropbox/Robinson Greenlees Eyes paper/Analysis")

#### Packages ######
require(ggplot2)
require(lme4)
require(gridExtra)
require(metafor)

#### Experiment one ####
### Get data ###
d=read.csv("Experiment one data.csv")
colnames(d)
# There were a few observations (12) from the old library and the students union
# but these were few and only in one treatment, so we exclude these
xtabs(~d$Location)
d=d[d$Location=="hershal"|d$Location=="merz court",]
d$Location=droplevels(d$Location)
# Code the actions verbally
d$Action=as.factor(d$Action)
levels(d$Action)=c("Threw on ground", "Kept on person", "Put elsewhere", "Left without removing")
xtabs(~d$Action)
xtabs(~d$Action+d$Condition)

# Plot of the basic frequencies
f1=ggplot(d, aes(x=Action, fill=Condition)) +
  geom_bar(position="dodge") +
  theme_bw() +
  ylab("Observations")
print(f1)

# Make dichotomous littering variable
# NB 'left without removing' is excluded
d$Litter=as.numeric(d$Action=="Threw on ground")
d$Litter[d$Action=="Left without removing"]=NA
xtabs(~d$Litter)/284

# Dichotomise number of people in vicinity
d$Someone=(d$Total.people>0)

# Plot of treatment and people in vicinity effect
f2=ggplot(d, aes(x=Someone, y=Litter, fill=Condition)) +
  stat_summary(fun.y=mean, geom="bar", position="dodge") +
  ylab("Proportion of people littering") +
  xlab("Someone else in vicinity") +
  coord_cartesian(ylim=c(0, 0.5)) +
  scale_fill_manual(values=c("black", "grey50")) +
  theme_bw()+
  geom_text(x=2.25, y=0.45, label="A")
print(f2)

# Statistical model
m1=glm(data=d, Litter~Condition*Someone, family=binomial(link=logit))
summary(m1)
m1null=glm(data=d, Litter~1, family=binomial(link=logit))
anova(m1, m1null, test="Chisq")

# Coefficient for eyes effect expressed as odds ratio
exp(m1$coefficients[2])
# Coefficient for someone in vicinity effect expressed as odds ratio
exp(m1$coefficients[3])

# Model with random effect of location
m2=glmer(data=d, Litter~Condition*Someone + (1|Location), family=binomial)
summary(m2)

##### Experiment 2 ####
s=read.csv("Experiment two data.csv")
levels(s$Action)=c("Put elsewhere", "Kept on person", "Left without removing", "Put in bin", "Threw on ground")
s$Condition2[s$Condition=="Big Eyes"]="Large eyes"
s$Condition2[s$Condition=="Control"]="Control"
s$Condition2[s$Condition=="Small eyes"]="Small eyes"
s$Condition=as.factor(as.character(s$Condition2))


# Plot of basic frequencies
fs1=ggplot(s, aes(x=Action, fill=Condition)) + 
  geom_bar(position="dodge") +
  theme_bw() + 
  ylab("Observations")
print(fs1)

# Make dichotomous littering variable
s$Litter=as.numeric(s$Action=="Threw on ground")
s$Litter[s$Action=="Left without removing"]=NA

# Make dichotomous someone in vicinity variable
s$Someone=(s$Total.people>0)

# Plot of treatment effects
fs2=ggplot(s, aes(x=Someone, y=Litter, fill=Condition)) +
  stat_summary(fun.y=mean, geom="bar", position="dodge") + 
  ylab("Proportion of people littering") +
  xlab("Someone else in vicinity") +
  coord_cartesian(ylim=c(0, 0.5)) +
  scale_fill_manual(values=c("black", "grey50", "grey80")) +
  theme_bw() +
  geom_text(x=2.25, y=0.45, label="B")
print(fs2)

# Statistical model
ms1=glm(data=s, Litter~Condition*Someone, family=binomial)
summary(ms1)
msnull=glm(data=s, Litter~1, family=binomial)
anova(ms1, msnull, test="Chisq")

# Delve into the interactions by running model with just (Someone==TRUE) and (Someone==FALSE)
ms1.2=glm(data=s, subset=(Someone==TRUE), Litter~Condition, family=binomial)
summary(ms1.2)
ms1.3=glm(data=s, subset=(Someone==FALSE), Litter~Condition, family=binomial)
summary(ms1.3)

# Now with random effect of location
ms2=glmer(data=s, Litter~Condition*Someone +(1|Location), family=binomial)
summary(ms2)

### Make figures ###
##### Figure 3 - plot of experimental effects ####
pdf("figure3.pdf", width=8, height=4)
grid.arrange(f2, fs2, ncol=2)
dev.off()

a=300
png("figure3.png", res=a, width=8*a, height=4*a)
grid.arrange(f2, fs2, ncol=2, widths=c(1,1.3))
dev.off()

#### Figure 4 meta-analysis plot ####

png("figure 4.png", res=300, width=2400, height=900)
forest(x=c(-1.13, -1.07, -1.51, -1.38, -3.12), sei=c(0.55, 0.47, 0.54, 0.57, 1.05), slab=c("Exp. 1: Eyes vs. control", "Exp. 2: Large eyes vs. control", "Exp. 2: Small eyes vs. control", "Exp. 1: Someone vs. no-one", "Exp. 2: Someone vs. no-one"), xlab="Odds ratio for littering", transf=exp, ref=1)
dev.off()

pdf("figure 4.pdf", width=8, height=3)
forest(x=c(-1.13, -1.07, -1.51, -1.38, -3.12), sei=c(0.55, 0.47, 0.54, 0.57, 1.05), slab=c("Exp. 1: Eyes vs. control", "Exp. 2: Large eyes vs. control", "Exp. 2: Small eyes vs. control", "Exp. 1: Someone vs. no-one", "Exp. 2: Someone vs. no-one"), xlab="Odds ratio for littering", transf=exp, ref=1)
dev.off()