##############################################################################
############ MIXED MODEL ANALYSIS FOR BREEDING SUCCESS ACCORDING #############
############                TO PARTNERSHIP STATUS               ##############
##############################################################################

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)

setwd("...")
DivBrS <- read.csv2("Pelletier2021_DivBrS_individual.csv")
str(DivBrS)
attach(DivBrS)
names(DivBrS)
head(DivBrS)

DivBrS$ID <- as.factor(DivBrS$ID)


#### 1. Data exploration

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

### 1.2 Explore data
DivBrS %>%
  group_by(Sex) %>%
  summarise(BrS = sum(BrS))

DivBrS %>%
  group_by(PartnStatus) %>%
  summarise(BrS = sum(BrS))

DivBrS %>%
  group_by(Yr_PrtnSt) %>%
  summarise(BrS = sum(BrS))


#### 2. GLMM
### 2.1 Full model
Model_Multi_Full <- glmer(BrS ~ PartnStatus*Yr_PrtnSt + (1 + PartnStatus + Yr_PrtnSt|ID),
                          family = binomial(logit),
                          data = DivBrS,
                          control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun=2e5)))
summary(Model_Multi_Full)

### 2.2 Model with intercept only
Model_Multi_Interc <- glmer(BrS ~ 1 + (1|ID),
                          family = binomial(logit),
                          data = DivBrS,
                          control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun=2e5)))
summary(Model_Multi_Interc)

### 2.3 Model with years only
Model_Multi_years <- glmer(BrS ~ (1 + Yr_PrtnSt|ID),
                            family = binomial(logit),
                            data = DivBrS,
                            control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun=2e5)))
summary(Model_Multi_years)

### 2.4 Model with partnership status only
Model_Multi_partnership <- glmer(BrS ~ (1 + PartnStatus|ID),
                           family = binomial(logit),
                           data = DivBrS,
                           control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun=2e5)))
summary(Model_Multi_partnership)

### 2.5 Models comparisons
anova(Model_Multi_Interc, Model_Multi_years, Model_Multi_partnership,
      Model_Multi_Full, test ="Chisq")

### 2.6 Summary of the best model
summary(Model_Multi_Full)

### 2.7 Post-hoc pairwise comparison, emmeans and plots
emmeans(Model_Multi_Full, pairwise~Yr_PrtnSt|PartnStatus, adjust="tukey")
emmeans(Model_Multi_Full, pairwise~PartnStatus|Yr_PrtnSt, adjust="tukey")
summ(Model_Multi_Full, exp = T)
plot(allEffects(Model_Multi_Full))
dotplot(ranef(Model_Multi_Full,condVar=TRUE),
        lattice.options=list(layout=c(1,2)))
qqmath(Model_Multi_Full)
plot(Model_Multi_Full,ID~resid(.))

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