setwd("X") ##set your own working directory

ls() ##shows active workspaces
rm(list=ls()) ##clears the active workspaces

library(irr) ##package for Fleiss kappa
library(lme4) ##package for glmm
library(ggplot2) ##package for plotting
library(plyr) ##package for plotting
library(lmerTest) ##package that shows significance for Gaussian GLMM

##FLEISS KAPPA
Set1Fleiss=read.table(file="Set1Fleiss.csv", header=T, stringsAsFactors=T, dec=".", sep=",")
str(Set1Fleiss)
kappam.fleiss(Set1Fleiss) ##calculate Fleiss kappa for all Set 1 videos

Set2Fleiss=read.table(file="Set2Fleiss.csv", header=T, stringsAsFactors=T, dec=".", sep=",")
str(Set2Fleiss)
kappam.fleiss(Set2Fleiss) ##calculate Fleiss kappa for all Set 2 videos

Set1MotherFleiss=read.table(file="Set1MotherFleiss.csv", header=T, stringsAsFactors=T, dec=".", sep=",")
str(Set1MotherFleiss)
kappam.fleiss(Set1MotherFleiss) ##calculate Fleiss kappa for Set 1 videos where the rater is a mother

Set1NonMotherFleiss=read.table(file="Set1NonMotherFleiss.csv", header=T, stringsAsFactors=T, dec=".", sep=",")
str(Set1NonMotherFleiss)
kappam.fleiss(Set1NonMotherFleiss) ##ccalculate Fleiss kappa for Set 1 videos where the rater is not a mother

Set2MotherFleiss=read.table(file="Set2MotherFleiss.csv", header=T, stringsAsFactors=T, dec=".", sep=",")
str(Set2MotherFleiss)
kappam.fleiss(Set2MotherFleiss) ##calculate Fleiss kappa for Set 1 videos where the rater is a mother

Set2NonMotherFleiss=read.table(file="Set2NonMotherFleiss.csv", header=T, stringsAsFactors=T, dec=".", sep=",")
str(Set2NonMotherFleiss)
kappam.fleiss(Set2NonMotherFleiss) ##ccalculate Fleiss kappa for Set 1 videos where the rater is not a mother


##PAIRED T-TEST FOR RATER CONFIDENCE
RaterConf=read.table(file="ConfidenceTTest.csv", header=T, stringsAsFactors=T, dec=".", sep=",")
str(RaterConf)
with(RaterConf, (t.test(High, Low, alternative='two.sided', 
	conf.level=.95, paired=TRUE)))
	

##GLMMS
CheckVShare=read.table(file="CheckVShare.csv", header=T, stringsAsFactors=T, dec=".", sep=",")
str(CheckVShare) ## this is used for all upcoming GLMMs unless otherwise stated

##GLMM FOR LOOK DURATION
LookTest=glmer(Infant.Look.Duration~Rater.Response+(1|Dyad.ID), 
	family=Gamma(link = "log"), data=CheckVShare)
summary(LookTest)

##GLMM FOR MUTUAL GAZE
MutGaze=glmer(Mutual.Gaze.Hurdle~Rater.Response+(1|Dyad.ID), 
	family=binomial(link = "logit"), data=CheckVShare)
summary(MutGaze)

##GLMM FOR INFANT FE
InfFE=glmer(Infant.FE~Rater.Response+(1|Dyad.ID), 
	family=binomial(link = "logit"), data=CheckVShare)
summary(InfFE)

##GLMM FOR INFANT GESTURE
InfGest=glmer(Infant.Gest~Rater.Response+(1|Dyad.ID), 
	family=binomial(link = "logit"), data=CheckVShare)
summary(InfGest)

##GLMM FOR INFANT VOC
InfVoc=glmer(Infant.Voc~Rater.Response+(1|Dyad.ID), 
	family=binomial(link = "logit"), data=CheckVShare)
summary(InfVoc)

##GLMM FOR MOTHER FE
MothFE=glmer(Mother.FE~Rater.Response+(1|Dyad.ID), 
	family=binomial(link = "logit"), data=CheckVShare)
summary(MothFE)

##GLMM FOR MOTHER GESTURE
MothGest=glmer(Mother.Gest~Rater.Response+(1|Dyad.ID), 
	family=binomial(link = "logit"), data=CheckVShare)
summary(MothGest)

##GLMM FOR MOTHER VOC
MothVoc=glmer(Mother.Voc~Rater.Response+(1|Dyad.ID), 
	family=binomial(link = "logit"), data=CheckVShare)
summary(MothVoc)

##GLMM FOR MOTHER COMMUNICATION FREQUENCY 
MothComm=lmer(Mother.Comm.Frequency~Rater.Response+(1|Dyad.ID), data=CheckVShare)
summary(MothComm)

##GLMM FOR INFANT COMMUNICATION FREQUENCY 
InfComm=lmer(Infant.Comm.Frequency~Rater.Response+(1|Dyad.ID), data=CheckVShare)
summary(InfComm)



##LOOK DURATION BARPLOT - figure 1
xdata=read.table(file="CheckVShare.csv", header=T, stringsAsFactors=T, dec=".", sep=",")
str(xdata)

lbar=ddply(xdata, c("Rater.Response"), summarise, 
	N=length(Infant.Look.Duration), mean =mean(Infant.Look.Duration), 
	sd=sd(Infant.Look.Duration), se=sd/sqrt(N))
bp=ggplot(lbar, aes(x=Rater.Response, y=mean, fill=Rater.Response))+
	ylab("Infant Look Duration (s)")+xlab("Type of Look")+
	geom_bar(position=position_dodge(), 
	stat="identity", width=.6)+geom_errorbar(aes(ymin=mean-se, 
	ymax=mean+se), width=.2, position=position_dodge(.9))
tiff("Fig1.tiff", units="in", width=3, height=2, res=300)
bp+guides(fill=FALSE)+ coord_cartesian(ylim = c(0, 4))+
	scale_y_continuous(expand = c(0, 0))+
	scale_fill_grey()+
	theme(axis.text.x=element_text(color="black", size=8),
	axis.title.x=element_text(size=8),
	axis.text.y=element_text(color="black", size=8),
	axis.title.y=element_text(size=8),
	panel.background=element_blank(),
	panel.border=element_blank(),
	axis.line=element_line(colour="black"))
dev.off()

##MOTHER COMMUNICATION FREQUENCY BARPLOT - figure 2
xdata=read.table(file="CheckVShare.csv", header=T, stringsAsFactors=T, dec=".", sep=",")
str(xdata)

mbar=ddply(xdata, c("Rater.Response"), summarise, 
	N=length(Mother.Comm.Frequency), mean =mean(Mother.Comm.Frequency), 
	sd=sd(Mother.Comm.Frequency), se=sd/sqrt(N))
bp=ggplot(mbar, aes(x=Rater.Response, y=mean, fill=Rater.Response))+
	ylab("Frequency of Mother Communication")+xlab("Type of Look")+
	geom_bar(position=position_dodge(), 
	stat="identity", width=.6)+geom_errorbar(aes(ymin=mean-se, 
	ymax=mean+se), width=.2, position=position_dodge(.9))
tiff("Fig2.tiff", units="in", width=3, height=2.5, res=300)
bp+guides(fill=FALSE)+ coord_cartesian(ylim = c(0, 4))+
	scale_y_continuous(expand = c(0, 0))+
	scale_fill_grey()+
	theme(axis.text.x=element_text(color="black", size=8), 
	axis.title.x=element_text(size=8),
	axis.text.y=element_text(color="black", size=8),
	axis.title.y=element_text(size=8),
	panel.background=element_blank(),
	panel.border=element_blank(),
	axis.line=element_line(colour="black"))
dev.off()


setwd("X") ##set your working directory
save.image("DetectingJAEvents.RData")
load("DetectingJAEvents.RData")