---
title: "R Markdown: Links between personality, reproductive success and re-pairing
  patterns in a long-lived seabird"
output:
  word_document:
    toc: yes
  html_document:
    toc: yes
    toc_float: yes
    theme: yeti
    highlight: haddock
  pdf_document:
    toc: yes
---

```{r, echo=FALSE}
knitr::read_chunk("R Markdown re-pairing script.R")

```

## Kiitiwake re-pairing model script

This R Markdown document details the code used to generate the results for a study on the links between personality (boldness), breeding outcome and re-pairing in black-legged kittiwakes. It covers the selection criteria used to create multiple datasets suitable for analysis in various different models used to analyse the data. 

## Data Setup, Descriptive Statistics and Variable Refinement

The following packages are used in this script. The data are uploaded to r studio. The categorical variables are defined as factors. 

```{r setup2, message=FALSE, warning=FALSE, results = 'hide'}

```

## M1: male boldness_female boldness- Assortative Mating

This code shows we tested for assortative mating by boldness at the Grumant and Pyramiden colonies. First the model is set up by transforming the male boldness scores with a square root transformation (following reflection) to correct for negative skew (and then returning the scores to the original direction) and scaling the female boldness scores. We also check for colinearity between the continuous variables.  

```{r assort_mating_setup, message=FALSE, warning=FALSE}
```

We then build a linear mixed model testing the relationships between male boldness score and female boldness score within pairs. We use a random effect of female ID to account for repeated measures.Colony and the interaction between colony and female boldness were included as fixed effects. The R squared values for the model are also produced. 

```{r assort_mating_model, message=FALSE, warning=FALSE}
```

The model selection criteria (see paper) are applied. This involves building every possible model and ranking them according to AICc. We remove models where delta is < 2 and nested models. 

```{r mating_dredge, message=FALSE, warning=FALSE}
```

Two models (with all fixed effects) is retained. This suggests there is evidence of assortative mating at Grumant, while mating at Pyramiden is weakly disassortative.

The data are extracted from the model for visualisation.

```{r mating_graph_data, message=FALSE, warning=FALSE}
```

A graph (Figure 2) is then constructed which shows how mating patterns differ between colonies, using original data and estimates from the best fitting model. 

```{r mating_graph, message=FALSE, warning=FALSE}
```

## M2: boldness_breeding outcome- Boldness and Reproductive Outcome

In this model the influence of the boldness metrics on reproductive success in year t are evaluated.To begin, the data are filtered so that any missing data from the variables required for this model are removed and the continuous variables are scaled.  

`male_pers_sc` = the scaled male's boldness score
`female_pers_sc` = the scaled female's boldness score
`absol_pers_sc` = the scaled absolute difference between pair members' boldness scores

```{r model_2_setup, message=FALSE, warning=FALSE, results = 'hide'}
```

Exploring the data gives us information on the boldness variables within and between sexes and colonies. We also retrieve details of the sample sizes included in the model.

```{r model_2_des, message=FALSE, warning=FALSE}
```

```{r model_2_sample_sizes, message=FALSE, warning=FALSE, results = 'hide'}
```

We check for multi-colinearity again. Any VIF scores greater than 5 would be a concern. These are all below that threshold.
```{r model_2_colin, message=FALSE, warning=FALSE}
```

We then build the global model. Here, `chick_sur` represents the binary breeding outcome for year t. We also generate R squared values for the model.     

```{r model_2_gm, message=FALSE, warning=FALSE}
```

We then employ model selection. Three models are retained. The first highlights absolute difference in boldness score, the second highlights female boldness score and the third is the null model.

*HOWEVER as the second model (highlighting female boldness) has an AICc equal to the null model, it can be labelled as nested and so we do not take it into consideration.*

```{r model_2_sel, message=FALSE, warning=FALSE}
```

Simplified models (based on the model selection criteria) are built to generate the data needed for visualisation. 

```{r model_2_sim_mod, message=FALSE, warning=FALSE, results = 'hide'}
```

To create the prediction graphs (Figure 3) the relevant data are extracted using sjPlot and used in tandem with the original data to show the trends. 

```{r model_2_sj, message=FALSE, warning=FALSE}
```

The model data and the original data and then combined to make the plot.

```{r model_2_g, message=FALSE, warning=FALSE}
```

## M3: breeding outcome_re-pairing outcome- Reproductive outcome and re-pairing probability 

In this model the influence of reproductive success in year t on re-pairing likelihood in year t + 1 is evaluated.To begin, the data are filtered so that any missing data from the variables required for this model are removed.

```{r model_3_setup, message=FALSE, warning=FALSE, results = 'hide'}
```

Exploring the data gives us details of the sample sizes included in the model.

```{r model_3_des, message=FALSE, warning=FALSE, results = 'hide'}
```

We then build the global model. Here, `chick_sur` represents the binary breeding outcome for year t. We also generate R squared values for the model.     

```{r model_3_gm, message=FALSE, warning=FALSE}
```

We then employ model selection. 

```{r model_3_sel, message=FALSE, warning=FALSE}
```

The best fitting model retains only chick survival as an influence on re-pairing probability. 

We build a simplified model which reflect the results of our model selection criteria.This will be used for visualisation. 

```{r model3_g, message=FALSE, warning=FALSE, results = 'hide'}
```

We use sjPlot to extract the data used in the graph (Figure 4).

```{r model3_sj, message=FALSE, warning=FALSE}
```

Variables are renamed to make the final graph clearer.A new scale is created for the response variable so it is clear on the final graph. 

```{r model3_rename, message=FALSE, warning=FALSE}
```

The final graph is produced. 

```{r model3_graph, message=FALSE, warning=FALSE}
```

## M4: re-pairing outcome_breeding outcome- Re-pairing and future reproductive outcomes

In this model the influence of re-pairing outcome on reproductive success in year t + 1 is evaluated.To begin, the data are filtered so that any missing data from the variables required for this model are removed. the categorical variables are turned into factors.

```{r mod4_setup, message=FALSE, warning=FALSE}
```

The sample sizes may be determined using the below code.

```{r mod4_ss, message=FALSE, warning=FALSE, results = 'hide'}
```

The global M4: re-pairing outcome_ breeding outcome is built. In the this model `next_season` is the breeding success in year t + 1 and `season_outcome` is the re-pairing outcome at the start of year t + 1.  

```{r mod4_global, message=FALSE, warning=FALSE}
```

Model selection takes place as before.

```{r model4_dredge, message=FALSE, warning=FALSE}
```

The null model is retained as the best fitting model. This suggests there is little evidence to suggest re-pairing at the start of year t impacts breeding success in that year. 

## M5: Boldness_re-pairing outcome- boldness and re-pairing probability

In this model the influence of pair members' boldness and absolute difference in boldness within pairs on re-pairing probability is evaluated. To begin, the data are filtered so that any missing data from the variables required for this model are removed. M5: boldness_re-pairing outcome is independent of breeding success in year t and thus the variable `ulti_outcome` can be used (i.e. did the pair dissolve at any point in the data collection period? Y/N). The `missed_seasons` variable accounts for the extra opportunities these birds had to re-pair in between sightings (see paper). Boldness metrics are defined as before.

```{r model5_setup, message=FALSE, warning=FALSE}
```

Descriptive reveal similar patterns to those reported under M2: boldness_breeding outcome. Males are bolder than females and birds are bolder at Grumant than Pyramiden.

```{r model5_des, message=FALSE, warning=FALSE}
```

Sample sizes can be determined using the below code. 

```{r model5_ss, message=FALSE, warning=FALSE, results = 'hide'}
```                                        

We test for colinearity between continuous variables. And VIF scores > 5 are problematic. There are no problems highlighted here. 

```{r model5_colin, message=FALSE, warning=FALSE}
``` 

The global M5: boldness_re-pairing outcome is built.

```{r mod5_global, message=FALSE, warning=FALSE}
``` 

Model selection occurs as before.

```{r mod5_ms, message=FALSE, warning=FALSE}
``` 

Only `missed_seasons` is retained and therefore there is little veidnece that the boldness metrics have an impact on re-pairing outcome.

## Permutation tests

These permutation tests assess if re-pairing birds choose a new partner which is more similar to themselves compared to their old partner. First, to determine the observed or ‘true’ difference in boldness scores between new partnerships and original partnerships, the absolute difference in boldness within each new partnership was subtracted from the absolute difference in boldness within the original partnership. A positive value indicated that the new partner was more similar in boldness to the focal bird than its original partner. 

```{r true_diff, message=FALSE, warning=FALSE}
``` 

The colonies are divided to account for the different personality distributions of potential mates at the two sites. 
```{r perm_sep, message=FALSE, warning=FALSE}
```

The loop parameters are set. This includes setting the number of permutations, creating new vectors and calculating the mean true difference. 

```{r loop_per, message=FALSE, warning=FALSE}
```

The loop is built. Boldness differences are calculated within randomised pairings. Within each permutation, every focal bird in the population was paired with a new partner boldness score, which is randomly sampled from the distribution of boldness scores available within each population. The absolute difference between the focal birds’ boldness and the randomly sampled personalities is calculated and subtracted from the original within pair difference in boldness. This value (within original pair difference – within random pair difference) is then averaged across the population to give a randomised mean difference (difference between the original and randomised pair’s absolute difference). This process is repeated 10,000 times per colony. The code `set.seed(100)` is used to ensure that the randomsiations return the same result every time for the purposes of repeatability.  

```{r random_gru, message=FALSE, warning=FALSE}
```

By calculating the proportion of cases where randomised personality values are equal or closer to focal birds' personality than true partners' personality (i.e. the null hypothesis)

```{r create_p_gru, message=FALSE, warning=FALSE}
```

These results can then be visualised and edited outside r for clarity (Figure 5). 

```{r perm_g_gru, message=FALSE, warning=FALSE}
```

This result suggests that birds at Grumant  choose new partners that were only slightly more similar to themselves than their previous mates, but that this was still significantly different from random, because the birds were mated assortatively to begin with. 

This process is then repeated using data for the Pyramiden birds.

```{r repeat_pyra, message=FALSE, warning=FALSE}
```

This result suggests that birds at Pyramiden may have re-paired at random. 