##############################################################################
######### MIXED MODEL ANALYSIS FOR DIVORCE RATE AT YEAR T ACCORDING ##########
#########        TO BREEDING SUCCESS AT YEAR T-1 AND YEAR           ##########
##############################################################################

library(lme4) # for multilevel models
library(tidyverse) # for data manipulation and plots
library(haven) #for reading sav data
library(sjstats) #for calculating intra-class correlation (ICC)
library(effects) #for plotting parameter effects
library(jtools) #for transformaing model summaries
library(ROCR) #for calculating area under the curve (AUC) statistics
library (multcomp)#for post-hoc analysis
library(emmeans)#for post-hoc analysis (BEST FOR NON-NORMAL DISTRIBUTION)
library(lattice)
library(EMAtools)
library(MuMIn)

setwd("...")
DivBrSt.1_Yrt_Partn <- read.csv("Pelletier2021_DivBrSt-1_yrt_partn_FORGLMM.csv", sep=";")

str(DivBrSt.1_Yrt_Partn) #where Divorce = 1 and retained = 0 (without lost) and BrS at t-1
attach(DivBrSt.1_Yrt_Partn)
names(DivBrSt.1_Yrt_Partn)
head(DivBrSt.1_Yrt_Partn)

DivBrSt.1_Yrt_Partn$ID <- as.factor(DivBrSt.1_Yrt_Partn$ID)
DivBrSt.1_Yrt_Partn$Divorce <- as.factor(DivBrSt.1_Yrt_Partn$Divorce)
DivBrSt.1_Yrt_Partn$Yrt <- as.factor(DivBrSt.1_Yrt_Partn$Yrt)
DivBrSt.1_Yrt_Partn$BrS <- as.factor(DivBrSt.1_Yrt_Partn$BrS)


#### 1. Data exploration

### 1.1 Inspect missing data
DivBrSt.1_Yrt_Partn %>%
  summarise_each(list(~sum(is.na(.)))) %>%
  gather()

### 1.2 Explore data
DivBrSt.1_Yrt_Partn %>%
  group_by(Yrt) %>%
  summarise(ID = length(unique(ID)))

DivBrSt.1_Yrt_Partn %>%
  group_by(BrS) %>%
  summarise(ID = length(BrS))


#### 2. GLMM

### 2.1  Test with dredge() of MuMin package for full model

#### 2.1.1 With Yrt as a random parameter and y = DivR
global.model.1 <- glmer(Divorce ~ BrS + (1|ID) + (1|Yrt),
                        data = DivBrSt.1_Yrt_Partn,
                        na.action = na.fail,
                        family = binomial(link = "logit"),
                        control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun=2e5)))
dredge(global.model.1, rank = "AICc")
summary(global.model.1)


#### 2.1.2 With Yrt as a fixed parameter and y = DivR
global.model.2 <- glmer(Divorce ~ BrS*Yrt + (1|ID),
                        data = DivBrSt.1_Yrt_Partn,
                        na.action = na.fail,
                        family = binomial(link = "logit"),
                        control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun=2e5)))
dredge(global.model.2, rank = "AICc")
summary(global.model.2)


#### 2.2 Models
mod1 <- glmer(Divorce ~ BrS*Yrt + (1|ID),
              data = DivBrSt.1_Yrt_Partn,
              na.action = na.fail,
              family = binomial(link = "logit"),
              control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun=2e5)))

mod2 <- glmer(Divorce ~ BrS + Yrt + (1|ID),
              data = DivBrSt.1_Yrt_Partn,
              na.action = na.fail,
              family = binomial(link = "logit"),
              control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun=2e5)))

mod3 <- glmer(Divorce ~ BrS + (1|ID),
              data = DivBrSt.1_Yrt_Partn,
              na.action = na.fail,
              family = binomial(link = "logit"),
              control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun=2e5)))

mod4 <- glmer(Divorce ~ Yrt + (1|ID),
              data = DivBrSt.1_Yrt_Partn,
              na.action = na.fail,
              family = binomial(link = "logit"),
              control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun=2e5)))

mod5 <- glmer(Divorce ~ BrS + (1|ID) + (1|Yrt),
              data = DivBrSt.1_Yrt_Partn,
              na.action = na.fail,
              family = binomial(link = "logit"),
              control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun=2e5)))



### 2.2 Model selection
model.sel(mod1, mod2, mod3, mod4, mod5,
          rank = AICc, extra = c(BIC, AICc))
anova(mod5, mod1, test ="Chisq")

### 2.3 Summary of the best model
summary(mod1)


### 2.4 Post-hoc pairwise comparison, emmeans and plots
emmeans_Div_BrS_Yrt <- emmeans(mod1, pairwise~BrS|Yrt, adjust="tukey")
emmeans_Div_BrS_Yrt

summ(mod1, exp = T)
plot(allEffects(mod1))
dotplot(ranef(mod1,condVar=TRUE),
        lattice.options=list(layout=c(1,2)))
qqmath(mod1)
plot(mod1,ID~resid(.))


### 2.5 AUR (area under the curve) for predicting Class with the model
Prob <- predict(mod1, type="response")
Pred <- prediction(Prob, as.vector(pull(DivBrSt.1_Yrt_Partn, Divorce)))
AUC <- performance(Pred, measure = "auc")
AUC <- AUC@y.values[[1]]
AUC       # Model discrimination is very good 82%