library(lme4)
library(glmmTMB)
library(car)


#####  PROXIMITY  ####
hist(x.data$Proximity)
x.data$pro.Proximity=(x.data$Proximity)/(x.data$Test_Duration)
x.data$tr.pro.Proximity=(x.data$pro.Proximity*(nrow(x.data) - 1) + 0.5)/nrow(x.data) 
#check which random slopes to include:
xx.fe.re=fe.re.tab(fe.model="tr.pro.Proximity~group*condition+sex+body_condition",re="(1|experimenter)", data=x.data)
xx.fe.re$summary[1:4] 
prox.data=xx.fe.re$data
prox.data$condition.happy=prox.data$condition.happy-mean(prox.data$condition.happy)
prox.data$condition.neutral=prox.data$condition.neutral-mean(prox.data$condition.neutral)
prox.data$sex.m=prox.data$sex.m-mean(prox.data$sex.m)
prox.data$body_condition.t=prox.data$body_condition.t-mean(prox.data$body_condition.t)
str(prox.data)
#model
fullP=glmmTMB(tr.pro.Proximity~ group*condition+sex+body_condition+
                (1|experimenter)+(0+condition.happy|experimenter)+(0+condition.neutral|experimenter)+
                (0+sex.m|experimenter)+(0+body_condition.t|experimenter),family=beta_family(link="logit"),data=prox.data)
summary(fullP)
#check for overdispersion:
overdisp.test(fullP) #ok
#absence of collinearity:
fullPc=lm(tr.pro.Proximity~ group+condition+sex+body_condition,data=x.data)
vif(fullPc) #ok
#check for BLUPs
ranef.diagn.plot(fullP)
#model stability
fullP.stab=glmmTMB.stab(model.res=fullP,para=T,data=prox.data)
table(fullP.stab$detailed$converged) 
round(fullP.stab$summary[,-1],3)
m.stab.plot(fullP.stab$summary[,-1])
#confident intervals
fullP.boot=boot.glmmTMB(model.res=fullP, data=prox.data,
                        excl.non.conv=F, nboots=1000, para=T, resol=100,
                        level=0.95, use=NULL, contr=NULL, n.cores=c("all-1", "all"))
fullP.boot$ci.estimates

#comparison red-full:testing the significance of the interaction
redP=glmmTMB(tr.pro.Proximity~ group+condition+sex+body_condition+
               (1|experimenter)+(0+condition.happy|experimenter)+(0+condition.neutral|experimenter)+
               (0+sex.m|experimenter)+(0+body_condition.t|experimenter),family=beta_family(link="logit"),data=prox.data)
anova(redP, fullP, test="Chisq") #the interaction is not significant
#comparison null-full: testing the effects of the predictors without the interaction
nullP=glmmTMB(tr.pro.Proximity~ 1+sex+body_condition+
                (1|experimenter)+(0+condition.happy|experimenter)+(0+condition.neutral|experimenter)+
                (0+sex.m|experimenter)+(0+body_condition.t|experimenter),family=beta_family(link="logit"),data=prox.data)
anova(nullP, fullP, test="Chisq")
#final model
summary(redP)

##### EATING AVAILABLE FOOD #####

#random slopes:
xx.fe.re=fe.re.tab(fe.model="Eat_all_Y1_N0 ~ group*condition+sex+body_condition",
                   re="(1|experimenter)", data=x.data)
xx.fe.re$summary[1:5] #
eat.data=xx.fe.re$data
str(eat.data)

eat.data$condition.happy=eat.data$condition.happy-mean(eat.data$condition.happy)
eat.data$condition.neutral=eat.data$condition.neutral-mean(eat.data$condition.neutral)
eat.data$sex.m=eat.data$sex.m-mean(eat.data$sex.m)
eat.data$body_condition.t=eat.data$body_condition.t-mean(eat.data$body_condition.t)

contr=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=100000))
fullEat.wac=glmer(Eat_all_Y1_N0~group*condition+sex+body_condition+
                    (1+condition.happy+condition.neutral+sex.m+body_condition.t|experimenter), 
                  data=eat.data, family=binomial, control=contr)
summary(fullEat.wac)$varcor
##all absolute correlations are estimated to be 1...
##try model without them:
fullEat.wnc=glmer(Eat_all_Y1_N0~group*condition+sex+body_condition+
                    (1+condition.happy+condition.neutral+sex.m+body_condition.t||experimenter), 
                  data=eat.data, family=binomial, control=contr)
logLik(fullEat.wac)
logLik(fullEat.wnc)
##not much of a difference, keep the model without the correlations
nullEat.wnc=glmer(Eat_all_Y1_N0~1+sex+body_condition+
                    (1+condition.happy+condition.neutral+sex.m+body_condition.t||experimenter), 
                  data=eat.data, family=binomial, control=contr)
ranef.diagn.plot(nullEat.wnc)
as.data.frame(anova(nullEat.wnc, fullEat.wnc))#the null-full comparison addresses the effect of group (in interaction with condition or as main effect), group or its interaction with condition has an effect on the response variable

redEat.wnc=glmer(Eat_all_Y1_N0~group+condition+sex+body_condition+
                   (1+condition.happy+condition.neutral+sex.m+body_condition.t||experimenter), 
                 data=eat.data, family=binomial, control=contr)
as.data.frame(anova(redEat.wnc, fullEat.wnc)) #the interaction is significant

summary(fullEat.wnc)

#model stability 
fullEat.wnc.stab=glmm.model.stab(model.res=fullEat.wnc, contr=contr, para=F,
                                 data=NULL)
table(fullEat.wnc.stab$detailed$warnings)
round(fullEat.wnc.stab$summary[, -1], 3)
is.re=grepl(x=rownames(fullEat.wnc.stab$summary), pattern="@")
m.stab.plot(fullEat.wnc.stab$summary[!is.re, -1])


#confident intervals (it takes long time)
plot.fullEat.wnc=glmer(Eat_all_Y1_N0~group*condition+sex.m+body_condition.t+(1+condition.happy+condition.neutral+sex.m+body_condition.t||experimenter), 
                       data=eat.data, family=binomial, control=contr)
plot.fullEat.wnc.boot=boot.glmm.pred(model.res=plot.fullEat.wnc,
                                     excl.warnings=F, use=c("group", "condition"),
                                     nboots=1000, para=F, resol=1000,
                                     level=0.95)
plot.fullEat.wnc.boot$ci.estimates


#####  TAIL WAGGING  ####
x.data$pro.Wagging=(x.data$Wagging)/(x.data$Test_Duration)
x.data$tr.pro.Wagging=(x.data$pro.Wagging*(nrow(x.data) - 1) + 0.5)/nrow(x.data) 
#random slopes
xx.fe.re=fe.re.tab(fe.model="tr.pro.Wagging~group*condition+sex+body_condition",re="(1|experimenter)", data=x.data)
xx.fe.re$summary[1:4] 
wag.data=xx.fe.re$data
wag.data$condition.happy=wag.data$condition.happy-mean(wag.data$condition.happy)
wag.data$condition.neutral=wag.data$condition.neutral-mean(wag.data$condition.neutral)
wag.data$sex.m=wag.data$sex.m-mean(wag.data$sex.m)
wag.data$body_condition.t=wag.data$body_condition.t-mean(wag.data$body_condition.t)
str(wag.data)

#model
fullW=glmmTMB(tr.pro.Wagging~ group*condition+sex+body_condition+
                (1+condition.happy+condition.neutral+sex.m+body_condition.t||experimenter),family=beta_family(link="logit"),data=wag.data)

#check for overdispersion
overdisp.test(fullW) 
#absence of collinearity (values should be <3)
fullWc=lm(tr.pro.Wagging~ group+condition+sex,data=x.data)
vif(fullWc) 
#check for BLUPs
ranef.diagn.plot(fullW)
#model stability 
fullW.stab=glmmTMB.stab(model.res=fullW,para=T,data=wag.data)
table(fullW.stab$detailed$converged) 
round(fullW.stab$summary[,-1],3)
m.stab.plot(fullW.stab$summary[,-1])
#confident intervals
fullW.boot=boot.glmmTMB(model.res=fullW, data=wag.data,
                       excl.non.conv=F, nboots=1000, para=T, resol=100,
                       level=0.95, use=NULL, contr=NULL, n.cores=c("all-1", "all"))
fullW.boot$ci.estimates

#comparison red-full:testing the significance of the interaction
redW=glmmTMB(tr.pro.Wagging~ group+condition+sex+body_condition+
               (1+condition.happy+condition.neutral+sex.m+body_condition.t||experimenter),family=beta_family(link="logit"),data=wag.data)
anova(redW, fullW, test="Chisq") #the interaction is not significant
#comparison full-null: testing the effects of the predictors without the interaction
nullW=glmmTMB(tr.pro.Wagging~ 1+sex+body_condition+
                (1+condition.happy+condition.neutral+sex.m+body_condition.t||experimenter),family=beta_family(link="logit"),data=wag.data)
anova(nullW, fullW, test="Chisq")
#final model
summary(redW)


#####  GAZE AVERSION  #####
#random slopes
xx.fe.re=fe.re.tab(fe.model="Avert_gaze ~ group * condition + sex+body_condition",
                   re="(1|experimenter)", data=x.data)
xx.fe.re$summary[1:5] 
gaze.data=xx.fe.re$data
gaze.data$condition.happy=gaze.data$condition.happy-mean(gaze.data$condition.happy)
gaze.data$condition.neutral=gaze.data$condition.neutral-mean(gaze.data$condition.neutral)
gaze.data$sex.m=gaze.data$sex.m-mean(gaze.data$sex.m)
gaze.data$body_condition.t=gaze.data$body_condition.t-mean(gaze.data$body_condition.t)

str(gaze.data)

#Including test duration in the dataframe
time=x.data$Test_Duration
gaze.data$Test_Duration=time

contr=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=100000))
fullGaze<-glmer.nb(Avert_gaze~group*condition+sex+body_condition+(1|experimenter)+
                     (0+condition.happy|experimenter)+(0+condition.neutral|experimenter)+
                     (0+sex.m|experimenter)+(0+body_condition.t|experimenter)+
                     offset(log(Test_Duration)),control=contr, family=poisson, data=gaze.data, na.action=na.omit)
overdisp.test(fullGaze)
ranef.diagn.plot(fullGaze) 
#collinearity
fullGazeColl<-lm(Avert_gaze~group+condition+sex+body_condition, data=x.data, na.action=na.omit)
vif(fullGazeColl) 
#model stability 
fullGaze.stab=glmm.model.stab(model.res=fullGaze, data=gaze.data)
table(fullGaze.stab$detailed$lme4.warnings)
table(fullGaze.stab$detailed$opt.warnings)
round(fullGaze.stab$summary[, -1], 3)
m.stab.plot(fullGaze.stab$summary[, -1])

#confident intervals #it takes long!
fullGaze.boot=boot.glmm.pred(model.res=fullGaze, excl.warnings=F,
                             nboots=1000, para=F, resol=1000, level=0.95, use=NULL)
fullGaze.boot$ci.estimates

#estimates
summary(fullGaze)
#comparison red-full:testing the significance of the interaction
redGaze<-glmer.nb(Avert_gaze~group+condition+sex+body_condition+(1|experimenter)+
                    (0+condition.happy|experimenter)+(0+condition.neutral|experimenter)+
                    (0+sex.m|experimenter)+(0+body_condition.t|experimenter)+
                    offset(log(Test_Duration)),control=contr, family=poisson, data=gaze.data, na.action=na.omit)
summary(redGaze)$varcor
anova(redGaze, fullGaze, test="Chisq")  #the interaction is not significant
#comparison full-null: testing the effects of the predictors without the interaction
nullGaze<-glmer.nb(Avert_gaze~1+sex+body_condition+(1|experimenter)+(0+condition.happy|experimenter)+
                     (0+condition.neutral|experimenter)+(0+sex.m|experimenter)+(0+body_condition.t|experimenter)+
                     offset(log(Test_Duration)),control=contr, family=poisson, data=gaze.data, na.action=na.omit)
anova(nullGaze, fullGaze, test="Chisq")#this is significant so I can remove the interaction from the full model and inspect the outcome
#final model
summary(redGaze)

##### LOOKING  #####
x.data$pro.LookE=(x.data$Looking)/(x.data$Test_Duration)
x.data$tr.pro.LookE=(x.data$pro.LookE*(nrow(x.data) - 1) + 0.5)/nrow(x.data) 
#random slopes
xx.fe.re=fe.re.tab(fe.model="tr.pro.LookE~group*condition+sex+body_condition",re="(1|experimenter)", data=x.data)
xx.fe.re$summary[1:5] 
look.data=xx.fe.re$data
look.data$condition.happy=look.data$condition.happy-mean(look.data$condition.happy)
look.data$condition.neutral=look.data$condition.neutral-mean(look.data$condition.neutral)
look.data$sex.m=look.data$sex.m-mean(look.data$sex.m)
look.data$body_condition.t=look.data$body_condition.t-mean(look.data$body_condition.t)
str(look.data)
#model
fullL=glmmTMB(tr.pro.LookE~ group*condition+sex+body_condition+
                (1+condition.happy+condition.neutral+sex.m+body_condition.t||experimenter),family=beta_family(link="logit"),data=look.data)
summary(fullL)
#check for overdispersion
overdisp.test(fullL) 
#absence of collinearity 
fullLc=lm(tr.pro.LookE~ group+condition+sex,data=x.data)
vif(fullLc) 
#check for BLUPs
ranef.diagn.plot(fullL)
#model stability 
fullL.stab=glmmTMB.stab(model.res=fullL,para=T,data=look.data)
table(fullL.stab$detailed$converged) 
round(fullL.stab$summary[,-1],3)
m.stab.plot(fullL.stab$summary[,-1])
#confident intervals
fullL.boot=boot.glmmTMB(model.res=fullL, data=look.data,
                        excl.non.conv=F, nboots=1000, para=T, resol=100,
                        level=0.95, use=NULL, contr=NULL, n.cores=c("all-1", "all"))
fullL.boot$ci.estimates
#comparison full-null: 
nullL=glmmTMB(tr.pro.LookE~ 1+sex+body_condition+
                (1+condition.happy+condition.neutral+sex.m+body_condition.t||experimenter),family=beta_family(link="logit"),data=look.data)
anova(nullL, fullL, test="Chisq")
#comparison red-full:testing the significance of the interaction
redL=glmmTMB(tr.pro.LookE~ group+condition+sex+body_condition+
               (1+condition.happy+condition.neutral+sex.m+body_condition.t||experimenter),family=beta_family(link="logit"),data=look.data)
anova(redL, fullL, test="Chisq") #the interaction is not significant

#final model
summary(redL)