---
title: "Appendix 2"
author: "Supplementary material of: Sex-specific fitness consequences of mate change in Scopoli's shearwaters (Calonectris diomedea)"
date: "2022-12-17"
output:
  html_document:
    toc: yes
    toc_depth: 3
    toc_float: yes
    collapsed: yes
    smooth_scroll: yes
    theme: readable
    highlight: kate
    df_print: paged
    code_folding: show
  '': default
bibliography: [biblio.bib]
csl: animal-behaviour.csl
link-citations: true
---

Loading the required libraries

```{r, warning=FALSE, message=FALSE}
library(glmmTMB)
library(MuMIn)
library(effects)
library(DHARMa)
library(performance)
library(bootComb)
library(ggplot2)
```

We have previously placed in the working directory the source files `skip.txt` and `breed.txt`.

# Hypothesis 1 - *Skipping behavior and mate change*

We test the *hypothesis 1*, that skipping behaviour depends on mate change. Note that `FIDELITY` in this study differs from `FIDELITY` in the subsequent breeding success analyses. In this analysis, an individual in state *changer* is when it switches partners in season t. However, it will be again in state *keeper* in the season t+1 if it maintains the same partner of season t. In the other breeding success analyses, an individual is in the *changer* state when it has changed its partner and this state does not revert back to *keeper* ever.

Loading and preparing the dataset.

```{r, warning=FALSE}
# Loading the data
datiskip<- read.table("skip.txt",header=T,sep="\t")

# Transforming some variables to categorical variables:
datiskip$ID<- as.factor(datiskip$ID)
datiskip$SEX<- as.factor(datiskip$SEX)
datiskip$NEST<- as.factor(datiskip$NEST)
datiskip$FIDELITY<- as.factor(datiskip$FIDELITY)# 0, same mate (first mate); 1, change mate
datiskip$SKIP<- as.factor(datiskip$SKIP)# 0, bred in previous season; 1, skipped previous season
datiskip$yearsfirst<- as.numeric(datiskip$yearsfirst)

# Removing the NAs:
datiskip0<- datiskip[,c("ID","SEX","YEAR","NEST","FIDELITY","SKIP", "yearsfirst"),]
datiskip0<- datiskip0[complete.cases(datiskip0),]
```

### Females

Preparing the dataset of females. We scale (mean=0,sd=1) the `yearsfirst` variable (reproductive age, i.e. years elapsed since the first capture) so that an individual of average age has value = 0.

```{r}
datiskip0f <- datiskip0[datiskip0$SEX=="2",]
datiskip0f$yearsfirst<- as.numeric(scale(datiskip0f$yearsfirst))

```

Model selection. We compare different crossed random intercepts structures for the model with `SKIP` (0, no skip; 1, skip) as response variable, `FIDELITY` (0, mate retention; 1, mate change) as focal predictor and `yearsfirst` as control predictor.

```{r}
fSkipFidAN<- glmmTMB( SKIP ~ FIDELITY + yearsfirst + (1|YEAR) + (1|NEST) , family = binomial, data = datiskip0f)
fSkipFidA<- glmmTMB( SKIP ~ FIDELITY + yearsfirst + (1|YEAR), family = binomial, data = datiskip0f)
fSkipFidI<- glmmTMB( SKIP ~ FIDELITY + yearsfirst + (1|ID) , family = binomial, data = datiskip0f)
fSkipFidNI <- glmmTMB( SKIP ~ FIDELITY + yearsfirst + (1|ID) + (1|NEST) , family = binomial, data = datiskip0f)
fSkipFidANI<- glmmTMB( SKIP ~ FIDELITY + yearsfirst + (1|YEAR) + (1|NEST) + (1|ID) , family = binomial, data = datiskip0f)
fSkipFidN<- glmmTMB( SKIP ~ FIDELITY + yearsfirst + (1|NEST) , family = binomial, data = datiskip0f)
fSkipFidAI<- glmmTMB( SKIP ~ FIDELITY + yearsfirst + (1|YEAR) + (1|ID) , family = binomial, data = datiskip0f)
fSkipFid<- glmmTMB( SKIP ~ FIDELITY + yearsfirst , family = binomial, data = datiskip0f)
```

We check for overdispersion, zero inflation and outliers for all models to identify the models to exclude from model selection.

```{r}
plot(simulateResiduals(fSkipFidAN))
plot(simulateResiduals(fSkipFidA))
plot(simulateResiduals(fSkipFidI))
plot(simulateResiduals(fSkipFidNI))
plot(simulateResiduals(fSkipFidANI))
plot(simulateResiduals(fSkipFidN))
plot(simulateResiduals(fSkipFidAI))
plot(simulateResiduals(fSkipFid))
```

All models meet the correct criteria in order to be considered for the model selection.

```{r}

model.sel(fSkipFidAN,fSkipFidA,fSkipFidI,fSkipFidNI,fSkipFidANI,fSkipFidN,fSkipFidAI,fSkipFid,rank = AICc)

```

The best model structure, according to model selection, is without any random intercept.

Now, we test whether the effect of mate change on the probability of skipping is significant.

```{r}
fSkipFidnull<- glmmTMB(SKIP ~ yearsfirst, family = binomial, data = datiskip0f)
anova(fSkipFid,fSkipFidnull)
fSkipFidnoyearsfirst<- glmmTMB(SKIP ~ FIDELITY, family = binomial, data = datiskip0f)
anova(fSkipFid,fSkipFidnoyearsfirst)
```

**The effect of mate change on the probability of skipping is positive (`r round(summary(fSkipFid)$coefficients$cond[2,1],2)`) and significant (P = `r round(anova(fSkipFid,fSkipFidnull)$Pr[2],3)`)**.

We plot the partial effects of the predictors on the probability of skipping and visualize the model's summary.

```{r}
plot(allEffects(fSkipFid))
summary(fSkipFid)
```

Now, we build the dataset to predict the probabilities of skipping breeding of females of average reproductive age (`yearsfirst=0`):

```{r}
newdatadatiskip0f<- data.frame(FIDELITY=c(0,1),yearsfirst=rep(0,2))
predsf<- predict(fSkipFid,newdatadatiskip0f, type="link", se.fit = T)
critval <- 1.96 ## approx 95% CI

fupr <- predsf$fit + (critval * predsf$se.fit)
flwr <- predsf$fit - (critval * predsf$se.fit)
ffit <- predsf$fit

ffit2 <- plogis(ffit)
fupr2 <- plogis(fupr)
flwr2 <- plogis(flwr)
```

The 95%CI estimate of the **probability of skipping for females who have NOT changed mate is `r round(ffit2[1],2)` (`r round(flwr2[1],2)` - `r round(fupr2[1],2)`)**. The 95%CI estimate of the **probability of skipping for females who have changed mate is `r round(ffit2[2],2)` (`r round(flwr2[2],2)` - `r round(fupr2[2],2)`)**.

### Males

Preparing the dataset of males. We scale (mean=0,sd=1) the `yearsfirst` variable (reproductive age, i.e. years elapsed since the first capture) so that an individual of average age has value = 0.

```{r}
datiskip0m <- datiskip0[datiskip0$SEX=="1",]
datiskip0m$yearsfirst<- as.numeric(scale(datiskip0m$yearsfirst))
```

Model selection. We compare different crossed random intercepts structures for the model with `SKIP` (0, no skip; 1, skip) as response variable, `FIDELITY` (0, mate retention; 1, mate change) as focal predictor and `yearsfirst` as control predictor.

```{r}
mSkipFidAN<- glmmTMB( SKIP ~ FIDELITY + yearsfirst + (1|YEAR) + (1|NEST) , family = binomial, data = datiskip0m)
mSkipFidA<- glmmTMB( SKIP ~ FIDELITY + yearsfirst + (1|YEAR), family = binomial, data = datiskip0m)
mSkipFidI<- glmmTMB( SKIP ~ FIDELITY + yearsfirst + (1|ID) , family = binomial, data = datiskip0m)
mSkipFidNI <- glmmTMB( SKIP ~ FIDELITY + yearsfirst + (1|ID) + (1|NEST) , family = binomial, data = datiskip0m)
mSkipFidANI<- glmmTMB( SKIP ~ FIDELITY + yearsfirst + (1|YEAR) + (1|NEST) + (1|ID) , family = binomial, data = datiskip0m)
mSkipFidN<- glmmTMB( SKIP ~ FIDELITY + yearsfirst + (1|NEST) , family = binomial, data = datiskip0m)
mSkipFidAI<- glmmTMB( SKIP ~ FIDELITY + yearsfirst + (1|YEAR) + (1|ID) , family = binomial, data = datiskip0m)
mSkipFid<- glmmTMB( SKIP ~ FIDELITY + yearsfirst    , family = binomial, data = datiskip0m)
```

We check for overdispersion, zero inflation and outliers for all models to identify the models to exclude from model selection.

```{r}
plot(simulateResiduals(mSkipFidAN))
plot(simulateResiduals(mSkipFidA))
plot(simulateResiduals(mSkipFidI))
plot(simulateResiduals(mSkipFidNI))
plot(simulateResiduals(mSkipFidANI))
plot(simulateResiduals(mSkipFidN))
plot(simulateResiduals(mSkipFidAI))
plot(simulateResiduals(mSkipFid))
```

Six models meet the correct criteria in order to be considered for the model selection.

```{r}
model.sel(mSkipFidAN,mSkipFidA,mSkipFidI,mSkipFidNI,mSkipFidN,mSkipFidAI,rank = AICc)
```

The best model structure is with `YEAR` as random intercept.

Now, we test whether the effect of mate change on the probability of skipping is significant.

```{r}
mSkipFidAnull<- glmmTMB(SKIP ~ yearsfirst + (1|YEAR), family = binomial, data = datiskip0m)
anova(mSkipFidA,mSkipFidAnull)
mSkipFidAnoyearsfirst<- glmmTMB(SKIP ~ FIDELITY + (1|YEAR), family = binomial, data = datiskip0m)
anova(mSkipFidA,mSkipFidAnoyearsfirst)# we don't find a significant effect of mate change on skipping behaviour.
```

**We do not find a significant effect of mate change on skipping behaviour** (`r round(summary(mSkipFidA)$coefficients$cond[2,1],2)`, *P = `r round(anova(mSkipFidA,mSkipFidAnull)$Pr[2],3)`*).

We plot the partial effects of the predictors on the probability of skipping and visualize the model's summary.

```{r}
plot(allEffects(mSkipFidA))
summary(mSkipFidA)
```

Now, we build the dataset to predict the probabilities of skipping breeding of females of average reproductive age (`yearsfirst=0`):

```{r}
newdatadatiskip0m<- data.frame(FIDELITY=c(0,1),yearsfirst=rep(0,2))
predsm<- predict(mSkipFid,newdatadatiskip0m, type="link", se.fit = T)
critval <- 1.96 ## approx 95% CI

mupr <- predsm$fit + (critval * predsm$se.fit)
mlwr <- predsm$fit - (critval * predsm$se.fit)
mfit <- predsm$fit

mfit2 <- plogis(mfit)
mupr2 <- plogis(mupr)
mlwr2 <- plogis(mlwr)
```

The 95%CI estimate of the **probability of skipping for males who have NOT changed mate is `r round(mfit2[1],2)` (`r round(mlwr2[1],2)` - `r round(mupr2[1],2)`)**. The 95%CI estimate of the **probability of skipping for males who have changed mate is `r round(mfit2[2],2)` (`r round(mlwr2[2],2)` - `r round(mupr2[2],2)`)**.

# Hypotheses 3 and 4 - *Breeding success and mate change*

Loading and preparing the dataset for the analysis.

```{r, warning=FALSE}
datifid<- read.table("breed.txt",header=T,sep="\t")
datifid$ID<- as.factor(datifid$ID)
datifid$sex<- as.factor(datifid$sex)
datifid$nest<- as.factor(datifid$nest)
datifid$fidelitybis<- as.factor(datifid$fidelitybis)# (0, with 1st mate; 1, with 2nd+ mate) 
datifid$bs<-as.factor(datifid$bs)# breeding success in the current year (t)
datifid$bs.before<-as.factor(datifid$bs.before)# breeding success in the previous year (t-1)
datifid$bs_change01<-as.numeric(datifid$bs_change01)# breeding success as the difference between previous and current (0 = from 1 to 0 or from 0 to 0; 1= from 0 to 1 or from 1 to 1)

# Removing the NAs:
datifid0<- datifid[,c("bs","fidelitybis","sex","ID","year","nest","yearsfirst"),]
datifid0<- datifid0[complete.cases(datifid0),]
datifid0$yearsfirst<- as.numeric(scale(datifid0$yearsfirst))# standardizing the variable with mean value = 0 and SD = 1.
datifid1<- datifid[,c("fidelitybis","sex","ID","year","nest","yearsfirst", "bs_change01"),]
datifid1<- datifid1[complete.cases(datifid1),]
datifid2<- datifid[,c("fidelitybis","sex","ID","year","nest","yearsfirst", "bs.before"),]
datifid2<- datifid2[complete.cases(datifid2),]

```

## Hypothesis 3 - First prediction 

We test a prediction of the *hypothesis 3*, that current breeding success depends on mate change.

### Females

Preparing the dataset.

```{r}
datifid0f <- datifid0[datifid0$sex=="2",]
table(datifid0f$bs,datifid0f$fidelitybis) 
```

Model selection. We compare different crossed random intercepts structures for the model with `bs` (0, breeding failure; 1, breeding success) as response variable, `fidelitybis` (0, mate retention; 1, mate change) as focal predictor and `yearsfirst` as control predictor.

```{r, warning=FALSE}
fbsa1AN<- glmmTMB(bs ~ fidelitybis + yearsfirst + (1|year) + (1|nest) , family = binomial, data = datifid0f)
fbsa1A<- glmmTMB(bs ~ fidelitybis + yearsfirst + (1|year), family = binomial, data = datifid0f)
fbsa1I<- glmmTMB(bs ~ fidelitybis + yearsfirst + (1|ID) , family = binomial, data = datifid0f)
fbsa1NI <- glmmTMB(bs ~ fidelitybis + yearsfirst + (1|ID) + (1|nest) , family = binomial, data = datifid0f)
fbsa1ANI<- glmmTMB(bs ~ fidelitybis + yearsfirst + (1|year) + (1|nest) + (1|ID) , family = binomial, data = datifid0f)
fbsa1N<- glmmTMB(bs ~ fidelitybis + yearsfirst + (1|nest) , family = binomial, data = datifid0f)
fbsa1AI<- glmmTMB(bs ~ fidelitybis + yearsfirst + (1|year) + (1|ID) , family = binomial, data = datifid0f)
fbsa1<- glmmTMB(bs ~ fidelitybis + yearsfirst , family = binomial, data = datifid0f)

```

We check for overdispersion, zero inflation and outliers for all models to identify the models to exclude from model selection.

```{r, warning=FALSE}
plot(simulateResiduals(fbsa1AN))
plot(simulateResiduals(fbsa1A))
plot(simulateResiduals(fbsa1I))
plot(simulateResiduals(fbsa1NI))
plot(simulateResiduals(fbsa1ANI))
plot(simulateResiduals(fbsa1N))
plot(simulateResiduals(fbsa1AI))
plot(simulateResiduals(fbsa1))
```

All models meet the correct criteria in order to be considered for the model selection.

```{r}
model.sel(fbsa1AN,fbsa1A,fbsa1I,fbsa1NI,fbsa1ANI,fbsa1N,fbsa1I,fbsa1AI,fbsa1,rank = AICc)
check_collinearity(fbsa1A)
```

The best model structure is with year as random intercept and there are no collinearity issues.

Now, we test whether the effect of mate change on the probability of breeding success is significant.

```{r}
fbsa1Anull<- glmmTMB(bs ~ yearsfirst + (1|year) , family = binomial, data = datifid0f)
anova(fbsa1A,fbsa1Anull)
fbsa1Anoyearsfirst<- glmmTMB(bs ~ fidelitybis + (1|year), family = binomial, data = datifid0f)
anova(fbsa1Anoyearsfirst,fbsa1A,test="chisq")
```

**We do not find a significant effect of mate change on breeding success** (`r round(summary(fbsa1A)$coefficients$cond[2,1],2)`, *P = `r round(anova(fbsa1A,fbsa1Anull)$Pr[2],3)`*), controlling for breeding experience (as proxied by `yearsfirst`). Also the effect of breeding experience (`yearsfirst`) on the breeding success is negligible.

We plot the partial effects of the predictors on the probability of skipping and visualize the model's summary.

```{r}
plot(allEffects(fbsa1A))
summary(fbsa1A)
```

### Males

Preparing the dataset.

```{r}
datifid0m <- datifid0[datifid0$sex=="1",]
table(datifid0m$bs,datifid0m$fidelitybis) 
```

Model selection. We compare different crossed random intercepts structures for the model with `bs` (0, breeding failure; 1, breeding success) as response variable, `fidelitybis` (0, mate retention; 1, mate change) as focal predictor and `yearsfirst` as control predictor.

```{r, warning=FALSE}
mbsaAN<- glmmTMB(bs ~ fidelitybis + yearsfirst + (1|year) + (1|nest) , family = binomial, data = datifid0m)
mbsaA<- glmmTMB(bs ~ fidelitybis + yearsfirst + (1|year), family = binomial, data = datifid0m)
mbsaI<- glmmTMB(bs ~ fidelitybis + yearsfirst + (1|ID) , family = binomial, data = datifid0m)
mbsaNI <- glmmTMB(bs ~ fidelitybis + yearsfirst + (1|ID) + (1|nest) , family = binomial, data = datifid0m)
mbsaANI<- glmmTMB(bs ~ fidelitybis + yearsfirst + (1|year) + (1|nest) + (1|ID) , family = binomial, data = datifid0m)
mbsaN<- glmmTMB(bs ~ fidelitybis + yearsfirst + (1|nest) , family = binomial, data = datifid0m)
mbsaAI<- glmmTMB(bs ~ fidelitybis + yearsfirst + (1|year) + (1|ID) , family = binomial, data = datifid0m)
mbsa<- glmmTMB(bs ~ fidelitybis + yearsfirst , family = binomial, data = datifid0m)

```

We check for overdispersion, zero inflation and outliers for all models to identify the models to exclude from model selection.

```{r, warning=FALSE}
plot(simulateResiduals(mbsaAN))
plot(simulateResiduals(mbsaA))
plot(simulateResiduals(mbsaI))
plot(simulateResiduals(mbsaNI))
plot(simulateResiduals(mbsaANI))
plot(simulateResiduals(mbsaN))
plot(simulateResiduals(mbsaAI))
plot(simulateResiduals(mbsa))
```

All models meet the correct criteria in order to be considered for the model selection.

```{r}
model.sel(mbsaAN,mbsaA,mbsaI,mbsaNI,mbsaANI,mbsaN,mbsaAI,mbsa,rank = AICc)
check_collinearity(mbsaA)
```

The best model structure is with year as random intercept and there are no collinearity issues.

Now, we test whether the effect of mate change on the probability of breeding success is significant.

```{r}
mbsaAnull<- glmmTMB(bs ~ yearsfirst + (1|year), family = binomial, data = datifid0m)
anova(mbsaA,mbsaAnull) 
mbsaAnoyearsfirst<- glmmTMB(bs ~ fidelitybis + (1|year), family = binomial, data = datifid0m)
anova(mbsaAnoyearsfirst,mbsaA,test="chisq")

```

**We do not find a significant effect of mate change on breeding success**(`r round(summary(mbsaA)$coefficients$cond[2,1],2)`, *P = `r round(anova(mbsaA,mbsaAnull)$Pr[2],3)`*), controlling for breeding experience (as proxied by `yearsfirst`). Also the effect of breeding experience (`yearsfirst`) on the breeding success is negligible.

We plot the partial effects of the predictors on the probability of skipping and visualize the model's summary.

```{r}
plot(allEffects(mbsaA))
summary(mbsaA)
```

## Hypothesis 3 - Second prediction

We test a prediction of the *hypothesis 3*, that the change between current and previous breeding success depends on mate change.

### Females

Preparing the dataset.

```{r}
datifid1f <- datifid1[datifid1$sex=="2",]
table(datifid1f$bs_change01,datifid1f$fidelitybis) 
```

Model selection. We compare different crossed random intercepts structures for the model with `bs` (0, breeding failure; 1, breeding success) as response variable, `fidelitybis` (0, mate retention; 1, mate change) as focal predictor and `yearsfirst` as control predictor.

```{r, warning=FALSE}
fbscANI<- glmmTMB( bs_change01 ~ fidelitybis + yearsfirst + (1|year) + (1|nest) +(1|ID) , family = binomial, data = datifid1f)
fbscAN<- glmmTMB( bs_change01 ~ fidelitybis + yearsfirst+ (1|year) + (1|nest) , family = binomial, data = datifid1f)
fbscA<- glmmTMB( bs_change01 ~ fidelitybis + yearsfirst+ (1|year) , family = binomial, data = datifid1f)
fbscNI<- glmmTMB( bs_change01 ~ fidelitybis + yearsfirst+ (1|nest) + (1|ID) , family = binomial, data = datifid1f)
fbscN<- glmmTMB( bs_change01 ~ fidelitybis + yearsfirst+ (1|nest) , family = binomial, data = datifid1f)
fbscAI<- glmmTMB( bs_change01 ~ fidelitybis + yearsfirst+ (1|year) + (1|ID) , family = binomial, data = datifid1f)
fbscI<- glmmTMB( bs_change01 ~ fidelitybis + yearsfirst+ (1|ID) , family = binomial, data = datifid1f)
fbsc<- glmmTMB( bs_change01 ~ fidelitybis + yearsfirst, family = binomial(), data = datifid1f)
```

We check for overdispersion, zero inflation and outliers for all models to identify the models to exclude from model selection.

```{r, warning=FALSE}
plot(simulateResiduals(fbscANI))
plot(simulateResiduals(fbscAN))
plot(simulateResiduals(fbscA))
plot(simulateResiduals(fbscNI))
plot(simulateResiduals(fbscN))
plot(simulateResiduals(fbscAI))
plot(simulateResiduals(fbscI))
plot(simulateResiduals(fbsc))

```

Only four models out of eight meet the correct criteria in order to be considered for the model selection.

```{r}
model.sel(fbscA,fbscNI,fbscN,fbscI,rank = AICc)
check_collinearity(fbscI)
```

The best model structure is with year and ID as random intercept and there are no collinearity issues.

Now, we test whether the effect of mate change on the probability of improving (or maintaining) the breeding success is significant.

```{r}
fbscInull<- glmmTMB(bs_change01 ~ yearsfirst +(1|ID), family = binomial(), data = datifid1f)
anova(fbscI,fbscInull) 
fbscInoyearsfirst<- glmmTMB( bs_change01 ~ fidelitybis + (1|ID), family = binomial, data = datifid1f)
anova(fbscInoyearsfirst,fbscI,test="chisq")
```

**We do not find a significant effect of mate change on breeding success change** (`r round(summary(fbscI)$coefficients$cond[2,1],2)`, *P = `r round(anova(fbscI,fbscInull)$Pr[2],3)`*), controlling for breeding experience (as proxied by `yearsfirst`). Also the effect of breeding experience (`yearsfirst`) on the breeding success is negligible.

We plot the partial effects of the predictors on the probability of skipping and visualize the model's summary.

```{r}
plot(allEffects(fbscI))
summary(fbscI)
```

### Males

Preparing the dataset.

```{r}
datifid1m <- datifid1[datifid1$sex=="1",]
table(datifid1m$bs_change01,datifid1m$fidelitybis) 
```

Model selection. We compare different crossed random intercepts structures for the model with `bs` (0, breeding failure; 1, breeding success) as response variable, `fidelitybis` (0, mate retention; 1, mate change) as focal predictor and `yearsfirst` as control predictor.

```{r, warning=FALSE}
mbscANI<- glmmTMB(bs_change01 ~ fidelitybis + yearsfirst + (1|year) + (1|nest) +(1|ID) , family = binomial, data = datifid1m)
mbscAN<- glmmTMB(bs_change01 ~ fidelitybis + yearsfirst+ (1|year) + (1|nest) , family = binomial, data = datifid1m)
mbscA<- glmmTMB(bs_change01 ~ fidelitybis + yearsfirst+ (1|year) , family = binomial, data = datifid1m)
mbscNI<- glmmTMB(bs_change01 ~ fidelitybis + yearsfirst+ (1|nest) + (1|ID) , family = binomial, data = datifid1m)
mbscN<- glmmTMB(bs_change01 ~ fidelitybis + yearsfirst+ (1|nest) , family = binomial, data = datifid1m)
mbscAI<- glmmTMB(bs_change01 ~ fidelitybis + yearsfirst+ (1|year) + (1|ID) , family = binomial, data = datifid1m)
mbscI<- glmmTMB(bs_change01 ~ fidelitybis + yearsfirst+ (1|ID), family = binomial, data = datifid1m)
mbsc<- glmmTMB(bs_change01 ~ fidelitybis + yearsfirst, family = binomial(), data = datifid1m)
```

We check for overdispersion, zero inflation and outliers for all models to identify the models to exclude from model selection.

```{r, warning=FALSE}
plot(simulateResiduals(mbscANI))
plot(simulateResiduals(mbscAN))
plot(simulateResiduals(mbscA))
plot(simulateResiduals(mbscANI))
plot(simulateResiduals(mbscN))
plot(simulateResiduals(mbscAI))
plot(simulateResiduals(mbscI))
plot(simulateResiduals(mbsc))
```

All models meet the correct criteria in order to be considered for the model selection.

```{r}
model.sel(mbscANI,mbscAN,mbscA,mbscNI,mbscN,mbscAI,mbscI,mbsc,rank = AICc)
check_collinearity(mbsc)
```

The best model structure is without random intercept and there are no collinearity issues.

Now, we test whether the effect of mate change on the probability of breeding success is significant.

```{r}
mbscnull<- glmmTMB(bs_change01 ~ yearsfirst, family = binomial(), data = datifid1m)
anova(mbsc,mbscnull)
mbscnoyearsfirst<- glmmTMB(bs_change01 ~ fidelitybis, family = binomial, data = datifid1m)
anova(mbscnoyearsfirst,mbsc,test="chisq")
```

**We do not find a significant effect of mate change on breeding success**(`r round(summary(mbsc)$coefficients$cond[2,1],2)`, *P = `r round(anova(mbsc,mbscnull)$Pr[2],3)`*), controlling for breeding experience (as proxied by `yearsfirst`). Also the effect of breeding experience (`yearsfirst`) on the breeding success is marginally significant.

We plot the partial effects of the predictors on the probability of skipping and visualize the model's summary.

```{r}
plot(allEffects(mbsc))
summary(mbsc)

```

## Hypothesis 4

We test a prediction of the hypothesis 4 (poor breeding success triggers divorce), that mate change depends increases after breeding failure in the previous season.

### Females

Preparing the dataset.

```{r}
datifid2f <- datifid2[datifid2$sex=="2",]
table(datifid2f$fidelitybis,datifid2f$bs.before) 
```

Unfortunately **this hypothesis cannot be tested with females' data**. There is only 1 case of failure for keepers, no failure for changers in the year before.

### Males

Preparing the dataset.

```{r}
datifid2m <- datifid2[datifid2$sex=="1",]
table(datifid2m$fidelitybis,datifid2m$bs.before) 
```

Model selection. We compare different crossed random intercepts structures for the model with `bs` (0, breeding failure; 1, breeding success) as response variable, `fidelitybis` (0, mate retention; 1, mate change) as focal predictor and `yearsfirst` as control predictor.

```{r, warning=FALSE}
mbsbAIN<- glmmTMB(fidelitybis ~ bs.before + yearsfirst + (1|year) + (1|ID) + (1|nest) , family = binomial, data = datifid2m)
mbsbIN<- glmmTMB(fidelitybis ~ bs.before + yearsfirst + (1|ID) + (1|nest) , family = binomial, data = datifid2m) #model convergence problem
mbsbAI<- glmmTMB(fidelitybis ~ bs.before + yearsfirst + (1|year) + (1|ID) , family = binomial, data = datifid2m)
mbsbAN<- glmmTMB(fidelitybis ~ bs.before + yearsfirst + (1|year) + (1|nest) , family = binomial, data = datifid2m)
mbsbI<- glmmTMB(fidelitybis ~ bs.before + yearsfirst + (1|ID), family = binomial, data = datifid2m)
mbsbA<- glmmTMB(fidelitybis ~ bs.before + yearsfirst + (1|year), family = binomial, data = datifid2m)
mbsbN<- glmmTMB(fidelitybis ~ bs.before + yearsfirst+ (1|nest) , family = binomial, data = datifid2m)
mbsb<- glmmTMB(fidelitybis ~ bs.before + yearsfirst, family = binomial, data = datifid2m)
```

We check for overdispersion, zero inflation and outliers for all models to identify the models to exclude from model selection.

```{r, warning=FALSE}
plot(simulateResiduals(mbsbAIN))
plot(simulateResiduals(mbsbIN))
plot(simulateResiduals(mbsbAI))
plot(simulateResiduals(mbsbAN))
plot(simulateResiduals(mbsbI))
plot(simulateResiduals(mbsbA))
plot(simulateResiduals(mbsbN))
plot(simulateResiduals(mbsb))
```

All models meet the correct criteria in order to be considered for the model selection except mbsbIN for model convergence problems.

```{r}
model.sel(mbsbAIN, mbsbIN, mbsbAN, mbsbAI, mbsbI, mbsbA, mbsbN, mbsb,rank = AICc)
```

The best model structure is with nest as random intercept, and there are no collinearity issues.

Now, we test whether the effect of mate change on the probability of breeding success is significant.

```{r}
mbsbNnull<- glmmTMB(fidelitybis ~ yearsfirst + (1|nest), family = binomial, data = datifid2m)
anova(mbsbN,mbsbNnull)
mbsbNnoyearsfirst<- glmmTMB(fidelitybis ~ bs.before + (1|nest), family = binomial, data = datifid2m)
anova(mbsbNnoyearsfirst,mbsbN,test="chisq")
```

**We do not find a significant effect of breeding success on mate change**(`r round(summary(mbsbN)$coefficients$cond[2,1],2)`, *P = `r round(anova(mbsbN,mbsbNnull)$Pr[2],3)`*).

We plot the partial effects of the predictors on the probability of skipping and visualize the model's summary.

```{r}
plot(allEffects(mbsbN))
summary(mbsbN)

```

# Divorce Rates

We calculate 95%CI and mean values of [@culina2013multievent] formula with bootstrap method for females and males. The values (low and high, `l`, `h` from the 95%CI values obtained from our lowest AICc capture-recapture models) of mortality (`Mort`calculated as 1-`Phi`) and mate change (`Psi`) are, for females (`F`) and males (`M`):

```{r}
MortFl <- 0.100507802; MortFh <- 0.205466633
MortMl <- 0.079348396; MortMh <- 0.143488436
PsiFl <- 0.053928372; PsiFh <- 0.161948775
PsiMl <- 0.085837216; PsiMh <- 0.222829484
```

### Females

```{r, warning=FALSE}
combFunEx<-function(pars){pars[[1]]-pars[[2]]}
bostF <- bootComb(distributions=c("beta","beta"),
                  qLowVect=c(PsiFl,MortMl), 
                  qUppVect=c(PsiFh,MortMh), 
                  combFun=combFunEx,
                  doPlot=TRUE,
                  method="hdi",
                  N=1e6, 
                  seed=1000);

normF<-getNormFromCI(qLow=bostF$conf.int[1],qUpp=bostF$conf.int[2])
xf<-seq(-0.1,0.15,length=1e3)
yf<-normF$d(xf)
plot(xf,yf,type="l",xlab="",ylab="density") # density plot for the fitted 
```

We calculate 95%CI and mean values of [@culina2013multievent] formula using the *delta method*.

```{r, warning=FALSE, results='hide'}
mortF<- 0.145289666#P mortality females
SEmortF<- 0.026582364#SE P mortality females
psiF<- 0.094985187#P of mate change females
SEpsiF<- 0.026772704# SE P mate change females
mortM<- 0.107271125#P mortality males
SEmortM<- 0.016236402#SE P mortality males
psiM<- 0.140952112#P of mate change males
SEpsiM<- 0.034481409# SE P mate change females

div.matrixF <- matrix( c(SEpsiF^2, 0, 0, SEmortM^2), 2, 2)
derivsF <- matrix( c(1, -1), 2, 1)
VarDivF <- t(derivsF) %*% div.matrixF %*% derivsF
divorceF <- psiF - mortM # Females' Divorce rate
SEdivF <- sqrt(VarDivF)
lcidivF <- divorceF - 1.96*SEdivF # 95% CI low
ucidivF <- divorceF + 1.96*SEdivF # 95% CI high

```

**Females' divorce rate**, truncated for negative values, is **`r round(divorceF,0)` (`r round(lcidivF,0)`-`r round(ucidivF,2)`)**.

### Males

```{r, warning=FALSE}
combFunEx<-function(pars){pars[[1]]-pars[[2]]}
bostM <- bootComb(distributions=c("beta","beta"),
                  qLowVect=c(PsiMl,MortFl), 
                  qUppVect=c(PsiMh,MortFh), 
                  combFun=combFunEx,
                  doPlot=TRUE,
                  method="hdi",
                  N=1e6, 
                  seed=352);

normM<-getNormFromCI(qLow=bostM$conf.int[1],qUpp=bostM$conf.int[2])
xm<-seq(-0.1,0.15,length=1e3)
ym<-normM$d(xm)
plot(xm,ym,type="l",xlab="",ylab="density") # density plot for the fitted normal distribution

```

We calculate 95%CI and mean values of [@culina2013multievent] formula using the *delta method*.

```{r, warning=FALSE, results='hide'}
mortF<- 0.145289666#P mortality females
SEmortF<- 0.026582364#SE P mortality females
psiF<- 0.094985187#P of mate change females
SEpsiF<- 0.026772704# SE P mate change females
mortM<- 0.107271125#P mortality males
SEmortM<- 0.016236402#SE P mortality males
psiM<- 0.140952112#P of mate change males
SEpsiM<- 0.034481409# SE P mate change females

div.matrixM <- matrix( c(SEpsiM^2, 0, 0, SEmortF^2), 2, 2)
derivsM <- matrix( c(1, -1), 2, 1)
VarDivM <- t(derivsM) %*% div.matrixM %*% derivsM
divorceM <- psiM - mortF # Males' Divorce rate
SEdivM <- sqrt(VarDivM)
lcidivM <- divorceM - 1.96*SEdivM # 95% CI low
ucidivM <- divorceM + 1.96*SEdivM # 95% CI high

```

**Males' divorce rate**, truncated for negative values, is **`r round(divorceM,0)` (`r round(lcidivM,0)`-`r round(ucidivM,2)`)**.

# Reproductive Life Expectancy (RLE)

### Females

These are the 95%CI survival estimates from the lowest AICc model without survival difference between *keepers* and *changers*, see Table 2, model F2 in main text.

```{r}
fem.phi<- 0.881653968
fem.phi.low<- 0.838120805
fem.phi.up<- 0.914671915
```

We calculate the 95%CI estimates of the *Reproductive Life Expectancy (RLE)* using the formula in [@seber1982estimation].

```{r}

RLSfems<- round(1/-log(fem.phi),2)
RLSfems.low<- round(1/-log(fem.phi.low),2)
RLSfems.up<- round(1/-log(fem.phi.up),2)

```

**Females' RLE is `r RLSfems` (`r RLSfems.low`-`r RLSfems.up`)**.

### Males

These are the 95%CI survival estimates from the lowest AICc model, see Table 1, model M1 in main text.

```{r}
mal.phi<- 0.89272888
mal.phi.low<- 0.85651156
mal.phi.up<- 0.9206516
```

We calculate the 95%CI estimates of the *Reproductive Life Expectancy (RLE)* using the formula in [@seber1982estimation].

```{r}

RLSmals<- round(1/-log(mal.phi),2)
RLSmals.low<- round(1/-log(mal.phi.low),2)
RLSmals.up<- round(1/-log(mal.phi.up),2)

```

**Males' RLE is `r RLSmals` (`r RLSmals.low`-`r RLSmals.up`)**.

# Simulation study

For each sex, we simulate a process where a cohort of breeders undergo the survival and mate change (from *keeper* to *changer*) process over 25 seasons after the reproductive recruitment. The process is stochastic and is driven by the parameters estimates (survival and mate change probabilities) from the lowest AICc capture-recapture models (model M1, Table 1 and model F1, Table 2). Then, we use the output of these simulations to calculate the expected number of individuals in each state (*keeper* or *changer*) at each season after the recruitment, which we use to calculate the **proportions of the cohort's individuals in state *keeper* and *changer***. Finally, we calculate **95%CI of the percentage of individuals (relative to the initial size of the cohort) which has changed mate over their lifetime**. The outputs of these simulations have been used to prepare Figure 2 and Figure S1.

## Constants

### Females

Survival and mate change probabilities from model F1 Table 2

```{r}
# female keeper survival P on the logit scale
logit.femk.phi<- 1.77203318
logit.femk.phi.se<- 0.21406217

# female changer survival P on the logit scale
logit.femc.phi<- 2.83462896
logit.femc.phi.se<- 0.65442989

# female mate change P on the logit scale
logit.fem.psi<- -2.25423036
logit.fem.psi.se<- 0.3114445

```

### Males

Survival and mate change probabilities from model M1 Table 1

```{r}
# male survival P on the logit scale
logit.mal.phi<- 2.11892341
logit.mal.phi.se<- 0.16954594

# male mate change P on the logit scale
logit.mal.psi<- -1.80740447
logit.mal.psi.se<- 0.28477119

```

### Simulation size values

```{r}
N<- 200 # cohort initial size, number of individuals at the beginning of the simulated process
seasons<- 50 # number of simulated seasons
n.sim<- 1000 # number of iterations 

```

### Containers to store outputs

These are the biological states in the simulated process:

-   dead = 0

-   alive with 1st mate = 1

-   alive with 2nd mate = 2

```{r}
# Container whithin each simulation iteration
fem.arrayNseasons<- matrix(NA,N,seasons)
fem.arrayNseasons[,1]<- 1 # all individuals are alive with 1st mate at the beginning of the simulated process

mal.arrayNseasons<- matrix(NA,N,seasons)
mal.arrayNseasons[,1]<- 1 # all individuals are alive with 1st mate at the beginning of the simulated process

# Containers of the simulation outputs
fem.arrayChangers<- vector()
fem.percChangers<- vector()

mal.arrayChangers<- vector()
mal.percChangers<- vector()

```

### To obtain reproducible results

```{r}
set.seed(1234)
setseeds<- sample(100000,n.sim)
```

## Females' simulation

```{r}
# SIMULATION PROCESS FEMALES
for (z in 1:n.sim){
  set.seed(setseeds[z])
for (i in 1:N){
  for (t in 2:seasons){
    # If an individual was dead (0) in t-1 assign 0 in t
    if (fem.arrayNseasons[i,t-1]==0){
      fem.arrayNseasons[i,t]<- 0
    } # if it was a 1 (alive with 1st mate), flip the coin to see if it survives to t! 
    else if (fem.arrayNseasons[i,t-1]==1) {
    logit.phi<- rnorm(1,logit.femk.phi,logit.femk.phi.se)
    phi<- plogis(logit.phi)
    aliveDead<- rbinom(1,1,phi)
    if (aliveDead==0) {# if it does not survive to t, assign 0
      fem.arrayNseasons[i,t]<- 0
    }
    else {# if it survives to t, flip the coin to see if it change mate  
      logit.psi<- rnorm(1,logit.fem.psi,logit.fem.psi.se)
      psi<- plogis(logit.psi)
      changed<- rbinom(1,1,psi)
      fem.arrayNseasons[i,t]<- changed+1 # assign 1 if it does not change mate, 2 if it changes mate
    }
  } else {# if it was a 2 (alive with 2nd mate), just flip the coin to see if it survives to t!
      logit.phi<- rnorm(1,logit.femc.phi,logit.femc.phi.se)
      phi<- plogis(logit.phi)
      aliveDead<- rbinom(1,1,phi)
      if (aliveDead==0) {# if it does not survive to t, assign 0
        fem.arrayNseasons[i,t]<- 0
      } else {# if it survives to t, assign 2
        fem.arrayNseasons[i,t]<- fem.arrayNseasons[i,t-1]
      }
    }
  }
}

  # fulfill the containers at the end of each iteration (z)
  x<- fem.arrayNseasons
  fem.arrayNseasonsChangers<- x[unique(which(x == 2 , arr.ind = TRUE)[,1]),]
  fem.arrayNseasonsChangers.1<- fem.arrayNseasonsChangers
  fem.arrayNseasonsChangers.1[fem.arrayNseasonsChangers.1 == 2]<- 1
  
  # This is for the data needed to calculate the 95% credible intervals of the percentage, relative to the cohort size, of individuals which have been changers.
  fem.percChangers[z]<- nrow(fem.arrayNseasonsChangers)/nrow(fem.arrayNseasons)
}

# Percentage and 95% credible intervals of females who change mate over their life
fem.percChange<- as.numeric(round(100*quantile(fem.percChangers,c(.5,.025,.975)),2))

```

According to the simulation, **`r fem.percChange[1]`% (`r fem.percChange[2]`% - `r fem.percChange[3]`%) of females have changed mate over their lives**.

## Males' simulation

```{r}
for (z in 1:n.sim){
  set.seed(setseeds[z])
  for (i in 1:N){
    for (t in 2:seasons){
      # If an individual was dead (0) in t-1 assign 0 in t
      if (mal.arrayNseasons[i,t-1]==0){
        mal.arrayNseasons[i,t]<- 0
      } # if it was a 1 (alive with 1st mate), flip the coin to see if it survives to t! 
      else if (mal.arrayNseasons[i,t-1]==1) {
        logit.phi<- rnorm(1,logit.mal.phi,logit.mal.phi.se)
        phi<- plogis(logit.phi)
        aliveDead<- rbinom(1,1,phi)
        if (aliveDead==0) {# if it does not survive to t, assign 0
          mal.arrayNseasons[i,t]<- 0
        }
        else {# if it survives to t, flip the coin to see if it change mate  
          logit.psi<- rnorm(1,logit.mal.psi,logit.mal.psi.se)
          psi<- plogis(logit.psi)
          changed<- rbinom(1,1,psi)
          mal.arrayNseasons[i,t]<- changed+1 # assign 1 if it does not change mate, 2 if it changes mate
        }
      } else {# if it was a 2 (alive with 2nd mate), just flip the coin to see if it survives to t!
        logit.phi<- rnorm(1,logit.mal.phi,logit.mal.phi.se)
        phi<- plogis(logit.phi)
        aliveDead<- rbinom(1,1,phi)
        if (aliveDead==0) {# if it does not survive to t, assign 0
          mal.arrayNseasons[i,t]<- 0
        } else {# if it survives to t, assign 2
          mal.arrayNseasons[i,t]<- mal.arrayNseasons[i,t-1]
        }
      }
    }
  }

  # fulfill the containers at the end of each iteration (z)
  x<- mal.arrayNseasons
  mal.arrayNseasonsChangers<- x[unique(which(x == 2 , arr.ind = TRUE)[,1]),]
  mal.arrayNseasonsChangers.1<- mal.arrayNseasonsChangers
  mal.arrayNseasonsChangers.1[mal.arrayNseasonsChangers.1 == 2]<- 1
  

  # This is for the data needed to calculate the 95% credible intervals of the percentage, relative to the cohort size, of individuals which have been changers.
  mal.percChangers[z]<- nrow(mal.arrayNseasonsChangers)/nrow(mal.arrayNseasons)
}

# Percentage and 95% credible intervals of malales who change mate over their life
mal.percChange<- as.numeric(round(100*quantile(mal.percChangers,c(.5,.025,.975)),2))

```

According to the simulation, **`r mal.percChange[1]`% (`r mal.percChange[2]`% - `r mal.percChange[3]`%) of males have changed mate over their lives**.

# Figure S1

```{r}
#| echo=FALSE, fig.cap="**Figure S1 -** Simulation results, for each sex, the proportions of the cohort's alive individuals in state *keeper* and *changer* which change over time.",
#| out.width = '100%'
knitr::include_graphics("Figure S1.png")
```



Note that the proportion of male *changers* increase faster than the proportion of female *changers* whose curve, however, decreases at a slower rate relecting their higher survival probability.



# Session Info

```{r}
sessionInfo()
```

# Bibliography