#Load data ----
xdata=read.table(file="corrected_dataset.txt", header = T, sep="\t", stringsAsFactors=T)

test=xdata=read.table(file="teste.txt", header = T, sep="\t", stringsAsFactors=T)
tail(xdata)


#Data exploration ----
plot(table(table(xdata$Dog)))


plot(table(xdata$S_Choice))
#number of times subject chose stag or hare

table(xdata$S_Choice, xdata$Group)
#number of times subject chose stag or hare by group


table(xdata$S_Choice, xdata$P_Choice, xdata$Group)
#subject's choice according to partner's choice, per group




#Analysis - Response and Match ----
##
##
##WE ARE RUNNING TWO MODELS FOR CHOICE, ONE WITH CONDITION AND ONE WITH PARTNER CHOICE


#SUBJECT'S CHOICE WITH PARTNER'S CHOICE ----
##finding the random slopes ----
source("diagnostic_fcns.r")
xx.fe.re1=fe.re.tab(fe.model= "S_Choice ~ Group* P_Choice * Trial + Session + Age + Chosen_side + Distance", 
                    re="(1|Dog)", other.vars = "Strategy", data=xdata)

xx.fe.re1$summary[1:25]
##session is a random slope
##trial is a random slope 
##chosen_side is a random slope
##distance is a random slope
##P_Choice is a random slope
##P_Choice:Trial is a random slope

##z transforming ----

r.data=xx.fe.re1$data
str(r.data)
#z.transform the covariates
r.data$z.Trial=as.vector(scale(r.data$Trial))
r.data$z.Session=as.vector(scale(r.data$Session))
r.data$z.Age=as.vector(scale(r.data$Age))
#centre factors for random slopes part
r.data$P_Choice.P_Stag=r.data$P_Choice.P_Stag-mean(r.data$P_Choice.P_Stag)
r.data$Chosen_side.Right=r.data$Chosen_side.Right-mean(r.data$Chosen_side.Right)
r.data$Distance.Far=r.data$Distance.Far-mean(r.data$Distance.Far)

##building the models ----
library(lme4)
contr=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))

full_ch=glmer(S_Choice ~ Group * P_Choice * 
                z.Trial + z.Session + z.Age + Chosen_side + Distance +
                (1 + z.Trial*P_Choice.P_Stag + z.Session + Chosen_side.Right + Distance.Far|Dog),
              data=r.data, family=binomial, control=contr)

#it converges, but we get a singular warning
summary(full_ch)$varcor #No issue with extreme correlations

null_ch=glmer(S_Choice ~ 
                z.Trial + z.Session + z.Age + Chosen_side + Distance +
                (1 + z.Trial*P_Choice.P_Stag + z.Session + Chosen_side.Right + Distance.Far|Dog),
              data=r.data, family=binomial, control=contr)

as.data.frame(anova(null_ch, full_ch, method="Chisq"))

#F-N is significant, there is an influence of the predictors Group and/or P_Choice and/or their interaction
#So now we can test which of these is happening, by testing each individual predictor separately, starting with the highest order interaction

round(summary(full_rf)$coefficients, 3)

source("C:/Users/Juliana/OneDrive/Documentos/Comportamento_Caes/Doutorado/Stats/functions/drop1p.r")

##test individual predictors ----
drop1(full_ch)
full_ch.test1=drop1p(model.res = full_ch, para = T, n.cores = "all-1", contr = contr)$drop1.res
full_ch.test1

#the 3 way interaction is not significant

full_ch.2=glmer(S_Choice ~ Group * P_Choice + 
                  z.Trial + z.Session + z.Age + Chosen_side + Distance +
                  (1 + z.Trial*P_Choice.P_Stag + z.Session + Chosen_side.Right + Distance.Far|Dog),
                data=r.data, family=binomial, control=contr)

full_ch.2.test1=drop1p(model.res = full_ch.2, para = T, n.cores = "all-1", contr = contr)$drop1.res

full_ch.2.test1

round(summary(full_ch.2)$coefficients, 3)

#interaction between group and P_choice is not significant

full_ch.3=glmer(S_Choice ~ Group + P_Choice + 
                  z.Trial + z.Session + z.Age + Chosen_side + Distance +
                  (1 + z.Trial*P_Choice.P_Stag + z.Session + Chosen_side.Right + Distance.Far|Dog),
                data=r.data, family=binomial, control=contr)

full_ch.3.test1=drop1p(model.res = full_ch.3, para = T, n.cores = "all-1", contr = contr)$drop1.res

full_ch.3.test1

round(summary(full_ch.3)$coefficients, 3)


library(emmeans)


emmeans(full_ch.3, pairwise~Group*P_Choice, by="Group")
emmeans(full_ch.3, pairwise~Group*P_Choice, by="P_Choice")

library(dplyr)
rdata %>% 
  na.omit() %>% 
  group_by(Group, P_Choice) %>% 
  summarise(ave.S_ch=mean(as.numeric(S_Choice)-1))

summary(st.full.2)$coefficients


aa=emmeans(full_rf.2,~P_Choice|Group)
bb=contrast(aa, method="pairwise")
summary(bb)



#RUNNING DIAGNOSIS ON THE MODEL ----

source("diagnostic_fcns.r")
ranef.diagn.plot(full_ch.3)
##I don't really remember how to interpret that

diagnostics.plot(full_ch.3)
##model stability
source("glmm_stability.r")
m.stab=glmm.model.stab(model.res=full_ch.3, contr=contr)
round(m.stab$summary[, -1],4)
m.stab.plot(m.stab$summary[, -1])
is.re=grepl(x=rownames(m.stab$summary), pattern="@")
m.stab.plot(m.stab$summary[!is.re, -1]) #fixed effects
m.stab.plot(m.stab$summary[is.re, -1]) #random effects
hist(m.stab$detailed$X.Intercept.)
View(m.stab$detailed)

summary(st.full.2)$coefficients

boot.res$ci.estimates

##
##
##
#MATCH WITH STRATEGY ----
##Random slopes ----

source("diagnostic_fcns.r")
xx.fe.re3=fe.re.tab(fe.model= "Match ~ Group * Strategy * Trial + Session + Age + Chosen_side + Distance", 
                    re="(1|Dog)", other.vars = "S_Choice", data=xdata)

xx.fe.re3$summary[1:25]

##session is a random slope
##trial is a random slope 
##chosen_side is a random slope
##distance is a random slope
##Strategy_within_Dog 
##Strategy:Trial_within_Dog 

##preparing the data for the model

st.data=xx.fe.re3$data
str(st.data)
#z.transform the covariates
st.data$z.Trial=as.vector(scale(st.data$Trial))
st.data$z.Session=as.vector(scale(st.data$Session))
st.data$z.Age=as.vector(scale(st.data$Age))
#centre factors for random slopes part
st.data$Strategy.Random=st.data$Strategy.Random-mean(st.data$Strategy.Random)
st.data$Strategy.Stag=st.data$Strategy.Stag-mean(st.data$Strategy.Stag)
st.data$Chosen_side.Right=st.data$Chosen_side.Right-mean(st.data$Chosen_side.Right)
st.data$Distance.Far=st.data$Distance.Far-mean(st.data$Distance.Far)


##Fitting full model ----
library(lme4)
contr=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))
st.full=glmer(Match ~ Group * Strategy * z.Trial + 
             z.Session + z.Age + Chosen_side + Distance +
             (1 + z.Trial*(Strategy.Random+Strategy.Stag) + z.Session + Chosen_side.Right + Distance.Far|Dog), 
           data=st.data, family=binomial, control = contr)
summary(st.full)$varcor #No issue with extreme correlations


##Fitting null model ----
st.null=glmer(Match ~ z.Trial + 
                z.Session + z.Age + Chosen_side + Distance +
                (1 + z.Trial*(Strategy.Random+Strategy.Stag) + z.Session + Chosen_side.Right + Distance.Far|Dog), 
              data=st.data, family=binomial, control = contr)
##F-N comparison ----
as.data.frame(anova(st.null, st.full, method="Chisq"))
#F-N is significant, there is an influence of the predictors Group and/or P_Choice and/or their interaction
#So now we can test which of these is happening, by testing each individual predictor separately, starting with the highest order interaction

##Test individual predictors ----
source("C:/Users/Juliana/OneDrive/Documentos/Comportamento_Caes/Doutorado/Stats/functions/drop1p.r")


st.full.test1=drop1p(model.res = st.full, para = T, n.cores = "all-1", contr = contr)$drop1.res
st.full.test1

st.full.2=glmer(Match ~ Group * Strategy 
                + z.Trial + 
                  z.Session + z.Age + Chosen_side + Distance +
                  (1 + z.Trial*(Strategy.Random+Strategy.Stag) + z.Session + Chosen_side.Right + Distance.Far|Dog), 
                data=st.data, family=binomial, control = contr)

st.full.2.test1=drop1p(model.res = st.full.2, para = T, n.cores = "all-1", contr = contr)$drop1.res
st.full.2.test1

round(summary(st.full.2)$coefficients, 3)

library(emmeans)

emmeans(st.full.2,~Group*Strategy, type="response")
emmeans(st.full.2, pairwise~Group*Strategy)
emmeans(st.full.2, pairwise~Group*Strategy, by="Group")
emmeans(st.full.2, pairwise~Group*Strategy, by="Strategy")

summary(st.full.2)$coefficients

#relevel
st.data$Group <-relevel(st.data$Group, ref="Control")
st.data$Strategy <-relevel(st.data$Strategy, ref="Hare")

st.full.icpt=glmer(Match ~ Group * Strategy 
                + z.Trial + 
                  z.Session + z.Age + Chosen_side + Distance +
                  (1 |Dog), 
                data=st.data, family=binomial, control = contr)
summary(st.full.icpt)$coefficients
#Check if intercept is significant, if it's not, it's more or less at random


---
library(lme4)
  
source("diagnostic_fcns.r")
ranef.diagn.plot(st.full.2)
##I don't really remember how to interpret that


##model stability
source("glmm_stability.r")
m.stab=glmm.model.stab(model.res=st.full.2, contr=contr)
m.stab.plot(m.stab$summary[, -1])

##CI
##Confidence intervals
source("boot_glmm.r")
boot.res=boot.glmm.pred(model.res=st.full.2, excl.warnings=T,
                        nboots=1000, para=T)
boot.res$ci.estimates