#Statistical Analyses on Gesture-Like Properties
#Looking at binary dataset (1 if property present; 0 if absent)
#Analyses conducted by Brittany Florkiewicz



################################################
###STEP 1: Import datasets. 
library("readxl")

AvgScores <- read_excel("C:/Users/Florkie/Google Drive/UCLA 2016/Dissertation/Dissertation Publications/Chapter 3 & 4 Paper/Submission 1/Florkiewicz and Campbell Data.xlsx", 
                        sheet = "Average Composite Gesture Score") #Import average composite gesture scores. 

GestComp <- read_excel("C:/Users/Florkie/Google Drive/UCLA 2016/Dissertation/Dissertation Publications/Chapter 3 & 4 Paper/Submission 1/Florkiewicz and Campbell Data.xlsx", 
                        sheet = "Raw Data") #Import raw composite gesture score data.

CTIScores <- read_excel("C:/Users/Florkie/Google Drive/UCLA 2016/Dissertation/Dissertation Publications/Chapter 3 & 4 Paper/Submission 1/Florkiewicz and Campbell Data.xlsx", 
                       sheet = "CTIScores") #Import CTI scores for signal types. 

CTIScores_PDWTest <- read_excel("C:/Users/Florkie/Google Drive/UCLA 2016/Dissertation/Dissertation Publications/Chapter 3 & 4 Paper/Submission 1/Florkiewicz and Campbell Data.xlsx", 
                        sheet = "CTIScores_PDWTest") #Import CTI scores for signal types which were also considered in Pollick and de Waal 2007. 
#Rename columns to make identification easier. 
names(GestComp)<-c("Video", "StartTime", "Modality", "SignalType", "Elaboration", "MechIneff", "Persistence", "I_InteractionOutcome", "ReceiverAttn", "RecipientID", "ResponseWait_E", "ResponseWait_P", "ResponseWait_O", "SignalerID", "F_InteractionOutcome", "PresumedGoal","GeneralizedBehavioralContext")
################################################



################################################
###STEP 2: Run a Mann-Whitney U Test on Average Composite Gesture Scores
wilcox.test(AvgScores$ACS ~ AvgScores$GT) #GT is whether signal type is a bodily gesture or not. 
################################################



################################################
###STEP 3: Run binomial GLMM's to compare the proportion of signals which contain each property. 
library("lme4")

#Run the models for facial signals (FE) versus bodily gestures (GE). 
RAM<-glmer(ReceiverAttn ~ Modality + (1|SignalerID), data=GestComp, family="binomial")
RID<-glmer(RecipientID ~ Modality + (1|SignalerID), data=GestComp, family="binomial")
RWEM<-glmer(ResponseWait_E ~ Modality + (1|SignalerID), data=GestComp, family="binomial")
RWPM<-glmer(ResponseWait_P ~ Modality + (1|SignalerID), data=GestComp, family="binomial")
RWOM<-glmer(ResponseWait_O ~ Modality + (1|SignalerID), data=GestComp, family="binomial")
PM<-glmer(Persistence ~ Modality + (1|SignalerID), data=GestComp, family="binomial")
EM<-glmer(Elaboration ~ Modality + (1|SignalerID), data=GestComp, family="binomial")
IIO<-glmer(I_InteractionOutcome ~ Modality + (1|SignalerID), data=GestComp, family="binomial")
FIO<-glmer(F_InteractionOutcome ~ Modality + (1|SignalerID), data=GestComp, family="binomial")
PG<-glmer(PresumedGoal ~ Modality + (1|SignalerID), data=GestComp, family="binomial")

#Run null models for comparison. 
RAM2<-glmer(ReceiverAttn ~ (1|SignalerID), data=GestComp, family="binomial")
RID2<-glmer(RecipientID ~ (1|SignalerID), data=GestComp, family="binomial")
RWEM2<-glmer(ResponseWait_E ~ (1|SignalerID), data=GestComp, family="binomial")
RWPM2<-glmer(ResponseWait_P ~ (1|SignalerID), data=GestComp, family="binomial")
RWOM2<-glmer(ResponseWait_O ~ (1|SignalerID), data=GestComp, family="binomial")
PM2<-glmer(Persistence ~ (1|SignalerID), data=GestComp, family="binomial")
EM2<-glmer(Elaboration ~ (1|SignalerID), data=GestComp, family="binomial")
IIO2<-glmer(I_InteractionOutcome ~ (1|SignalerID), data=GestComp, family="binomial")
FIO2<-glmer(F_InteractionOutcome ~ (1|SignalerID), data=GestComp, family="binomial")
PG2<-glmer(PresumedGoal ~ (1|SignalerID), data=GestComp, family="binomial")

#Compare null models to models with modality to determine whether modality has a significant effect on proportions observed. 
anova(RAM,RAM2, test="Chisq") #Significant
anova(RID,RID2, test="Chisq") #Significant
anova(RWEM,RWEM2, test="Chisq") #Significant
anova(RWPM,RWPM2, test="Chisq") #Significant
anova(RWOM,RWOM2, test="Chisq") #NOT Significant
anova(PM,PM2, test="Chisq") #Significant
anova(EM,EM2, test="Chisq") #Significant
anova(IIO,IIO2, test="Chisq") #Significant
anova(FIO,FIO2, test="Chisq") #NOT Significant
anova(PG,PG2, test="Chisq") #Significant

#Get summaries of regular models.  
summary(RAM)
summary(RID)
summary(RWEM) 
summary(RWPM) 
summary(RWOM) 
summary(PM)
summary(EM)
summary(IIO)
summary(FIO) 
summary(PG)

#Get the odds ratio for each of the models.
library("emmeans")
data.frame(confint(pairs(emmeans(RAM, ~ Modality,type="response"))))
data.frame(confint(pairs(emmeans(RID, ~ Modality,type="response"))))
data.frame(confint(pairs(emmeans(RWEM, ~ Modality,type="response"))))
data.frame(confint(pairs(emmeans(RWPM, ~ Modality,type="response"))))
data.frame(confint(pairs(emmeans(RWOM, ~ Modality,type="response"))))
data.frame(confint(pairs(emmeans(PM, ~ Modality,type="response"))))
data.frame(confint(pairs(emmeans(EM, ~ Modality,type="response"))))
data.frame(confint(pairs(emmeans(IIO, ~ Modality,type="response"))))
data.frame(confint(pairs(emmeans(FIO, ~ Modality,type="response"))))
data.frame(confint(pairs(emmeans(PG, ~ Modality,type="response"))))

#Look at conditional effects for FE vs. GE.
library("effects")

eff1<- effect("Modality", RAM)
eff2<- effect("Modality", RID)
eff3<- effect("Modality", RWEM)
eff4<- effect("Modality", RWPM)
eff5<- effect("Modality", RWOM)
eff6<- effect("Modality", PM)
eff7<- effect("Modality", EM)
eff8<- effect("Modality", IIO)
eff9<- effect("Modality", FIO)
eff10<- effect("Modality", PG)

plot(eff1) #GE
plot(eff2) #GE
plot(eff3) #GE
plot(eff4) #FE
plot(eff5) #None.
plot(eff6) #FE
plot(eff7) #FE
plot(eff8) #GE
plot(eff9) #None. 
plot(eff10) #GE 
################################################



################################################
###STEP 4: Calculate CTI scores (based on Pollick and de Waal 2007) and run a Mann-Whitney U test.
wilcox.test(CTIScores$CTI ~ CTIScores$MODE) #In mode, 1 is gesture and 0 is facial signal. 

#Check to see how this compares to the subset of data presented in Pollick and de Waal 2007 in figure 3 (those select signal types)
wilcox.test(CTIScores_PDWTest$CTI ~ CTIScores_PDWTest$MODE) #In mode, 1 is gesture and 0 is facial signal. 
################################################



################################################
###STEP 5: Run binomial tests to determine if the number of signals observed having each gesture-like property is significant. 

#Facial signals (n=1090)
binom.test(1090, 1090, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Mechanically ineffective 
binom.test(535, 1090, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Receiver Attention
binom.test(827, 1090, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Recipient ID
binom.test(200, 1090, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Response Waiting - End
binom.test(678, 1090, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Response Waiting - Persist
binom.test(702, 1090, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Response Waiting - Overall
binom.test(975, 1090, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Persistence
binom.test(1024, 1090, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Elaboration
binom.test(778, 1090, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Immediate IO
binom.test(961, 1090, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Final IO
binom.test(772, 1090, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Presumed goal

#Bodily gestures (n=2356)
binom.test(2331, 2356, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Mechanically ineffective 
binom.test(1389, 2356, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Receiver Attention
binom.test(2261, 2356, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Recipient ID
binom.test(594, 2356, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Response Waiting - End
binom.test(1413, 2356, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Response Waiting - Persist
binom.test(1692, 2356, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Response Waiting - Overall
binom.test(1667, 2356, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Persistence
binom.test(2053, 2356, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Elaboration
binom.test(1902, 2356, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Immediate IO
binom.test(2107, 2356, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Final IO
binom.test(2151, 2356, p=0.5, alternative=c("two.sided"), conf.level=0.95) #Presumed goal
################################################



################################################
###STEP 6: Run a Mann-Whitney U for composite gesture scores and CTI scores to compare facial expressions with vocalizations to facial expressions without.

FEV_vs_FENV <- read_excel("C:/Users/Florkie/Google Drive/UCLA 2016/Dissertation/Dissertation Publications/Chapter 3 & 4 Paper/Submission 1/Florkiewicz and Campbell Data.xlsx", 
                        sheet = "FE_V vs FE_NV") #Import average CTI scores and CGS scores for facial expressions.

wilcox.test(FEV_vs_FENV$`CTI Scores` ~ FEV_vs_FENV$`With Vocalizaiton? (1=yes, 0=no)`) #CTI SCORES 
wilcox.test(FEV_vs_FENV$`CGS Scores` ~ FEV_vs_FENV$`With Vocalizaiton? (1=yes, 0=no)`) #CGS SCORES 
################################################