---
title: "Waal Intrinsic Variables and Parental Care Coordination Model Script"
output: 
  html_document:
   toc : true
   toc_float : true
   theme: yeti
   highlight: haddock
---

```{r, echo=FALSE}
knitr::read_chunk("NEW anon WAAL model R markdown script.R")

```

# Wandering Albatross Model Script

This R Markdown document details the code used to generate the results for a study on the impact of intrinsic variables on parental care coordination in wandering albatrosses. It covers the refinement of the data to create a data set suitable for analysis and the mixed 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'}

```

The full data set is divided into the incubation and brooding breeding stages. Each of these subsets will be analysed in a separate model. 

```{r div, message=FALSE, warning=FALSE, echo=TRUE}

```

The following code relates to generating the descriptive statistics for the incubation and brooding subsets, respectively. This includes examining sample sizes within these subsets and examining measures of centrality. 

### Descriptive Statistics: Incubation 
Raw sample sizes for the incubation subset:
```{r des-incu, message=FALSE, warning=FALSE, echo=TRUE, results = 'hide'}

```

Sample sizes per year and age group for the incubation subset:
```{r sam-sp-in, message=FALSE, warning=FALSE, echo=TRUE}

```

Measures of centrality and group comparison for the incubation subset:
```{r cent-incu, message=FALSE, warning=FALSE, echo=TRUE}

```

### Descriptive Statistics: Brooding
Raw sample sizes for the brooding subset:

```{r des-brooding, message=FALSE, warning=FALSE, echo=TRUE, results = 'hide'}

```

Sample sizes per year and age group for the brooding subset:
```{r sam-sp-br, message=FALSE, warning=FALSE, echo=TRUE}

```

Measures of centrality and group comparison for the brooding subset:
```{r cent-brood, message=FALSE, warning=FALSE, echo=TRUE}

```

### Quadratic Variables, Scaling and Partner's Previous Trip Adjustment

Certain variables need to be adjusted for analysis. The instructions below below include making a new data frame, calculating a quadratic age variable and adjusting the partner's previous trip variable to account for each individual's mean previous trip duration. 

### Variable Refinement: Incubation
Step 1: First, a new data frame is created to separate alterations from original data set. A check for NAs is carried out simultaneously (there should be none).The following actions are labeled as steps 1-4 to make comparison with the brooding code easier.

```{r incu-data, message=FALSE, warning=FALSE, echo=TRUE}

```

Step 2: A quadratic variable representing age must be created. To prevent low sample sizes at the extremes of the age spectrum causing inaccuracies in the analyses, the highest and lowest age categories are collapsed so there is a minimum of 5 individuals at each end of the spectrum (see methods for further details and justification). This is done for both the focal bird age and also the partner bird age variables.

```{r incu-age, message=FALSE, warning=FALSE, echo=TRUE}

```

Comparing the new samples size per age category with those above shows there are now at least five individuals in each of the most extreme age categories. 

```{r check-age, message=FALSE, warning=FALSE, echo=TRUE}

```

step 3: In this important section, the partner's previous bout variable is refined to account for each individual's mean previous foraging trip duration. This helps us to ensure that the focal bird is responding to its partner's most recent trip, rather than their general foraging pattern. This is accomplished by creating a new variable which represents the differences between each partner's previous trip and their mean previous trip length. Step-by-step instructions for each line of code are included in the comments. 

```{r incu-ppb, message=FALSE, warning=FALSE, echo=TRUE}

```
### Variable Refinement: Brooding

Steps 1-3 above are then repeated, this time using the brooding data subset. 

Step 1: a new data frame is created:
```{r brood-data, message=FALSE, warning=FALSE, echo=TRUE}

```

Step 2: the quadratic age and partner age variables are created.
```{r brood-age, message=FALSE, warning=FALSE, echo=TRUE}

```

and checked:
```{r brood-check, message=FALSE, warning=FALSE, echo=TRUE}

```

Step 3: the new variable representing partner's previous trip duration is created.
```{r brood-ppb, message=FALSE, warning=FALSE, echo=TRUE}

```


The two breeding stage subsets and recombined so that the variables can be scaled across the entire dataset. This will allow coordination strength to be directly compared between models and thus, between breeding stages. The data are then redivided into the subsets for analysis. 
```{r data-scale, message=FALSE, warning=FALSE, echo=TRUE}

```



## Assumption Testing 

Both data subsets are tested individually to ensure they meet the assumptions required by linear mixed modeling. This involves examining:

The normality in the response variable
Homogeneity of variance 
Multicolinearity 

### Assumption Testing: Incubation

As in the variable refinement section, the assumption tests for incubation will be labeled as steps 1-3 to make comparison with the brooding code easier to follow. 

Step 1: the response variable representing the length of the focal bird's foraging trip (`no_hours`) is not normally distributed and requires transformation.
This can be seen using an initial histogram and is corrected by creating a square root variable (`no_hours_sqrt`).

```{r norm-incu, message=FALSE, warning=FALSE, echo=TRUE}

```

The variable `no_hours_sqrt` will be used as the response variable in the models.

Step 2: the following code is used to check each variable for homogeneity of variance. Output plots are omitted here for brevity.  

```{r incu-hov, message=FALSE, warning=FALSE, echo=TRUE, results = FALSE, fig.show='hide'}

```

Step 3: this code searches for evidence of multicolinearity between potentially correlated continuous fixed variables. VIF scores > 3 are considered problematic. Code for the production of graphs has been provided but outputs are omitted for brevity.  
```{r co-in, message=FALSE, warning=FALSE, echo=TRUE, fig.show = 'hide'}

```

### Assumption Testing: Brooding

Once again, steps 1-3 above are repeated this time for the brooding subset. 

Step 1: transformation of the response variable into `no_hour_sqrt`. 

```{r brood-norm, message=FALSE, warning=FALSE, echo=TRUE}

```

The variable `no_hours_sqrt` will be used as the response variable in the models.

Step 2: check each variable for homogeneity of variance. Output plots are omitted here for brevity.  

```{r brood-hov, message=FALSE, warning=FALSE, echo=TRUE, results = FALSE, fig.show='hide'}

```

Step 3: examine  evidence of multicolinearity between potentially correlated continuous fixed variables. VIF scores > 3 are considered problematic. Code for the production of graphs has been provided but outputs are omitted for brevity.  

```{r brood-co, message=FALSE, warning=FALSE, echo=TRUE, fig.show = 'hide'}

```

## Create a total dataset and save files
An updated, combined data set of both the incubation and brooding data is created. 

```{r combine, message=FALSE, warning=FALSE, echo=TRUE}

```

Sample size and centrality measures can be calculated for the record.

```{r total_sample, message=FALSE, warning=FALSE, echo=TRUE, results = FALSE, fig.show='hide'}

```

All files (incubation, brooding and total observations) are saved for later use in graph creation.
```{r save, message=FALSE, warning=FALSE, echo=TRUE}

```

## Linear Mixed Models and Model Selection 

The following code outlines how the linear mixed models (LMMs) were created and how information theory (specifically AIC and model averaging) was applied for the purpose of model selection. Two models, one for the incubation subset and one for the brooding subset, were created. Steps 1-6 will be used to aid comparison between models. 

### Incubation Models

Step 1: As suggested by Bolker et al (2009), a basic linear model was created ahead of the LMM.

```{r inlinmod, message=FALSE, warning=FALSE, echo=TRUE}

```

Step 2: The global model is then created. It features fixed effects, relevant interactions and random intercepts and slopes (see methods for details). Optimix optimizer is included to aid convergence. The residuals are also checked at this point. 

```{r incomod, message=FALSE, warning=FALSE, echo=TRUE}

```

Step 3: A model set is created from the global model. This creates all possible versions of the model to establish the best fit. As the focal bird age and its quadratic equivalent (`age_sc` and `age_sq_sc`) must be retained or discarded together, additional code is added to ensure this is reflected in the model set. Similarly, the partner bird age variable and its quadratic equivalent(`p_age_sc` and `p_age_sq_sc`) receive the same treatment. This portion of the code may take several hours to run.  

```{r in-dredge, message=FALSE, warning=FALSE, echo=TRUE, results = FALSE}

```

Step 4: The models are ranked according to AIC. Next, those models with delta < 2 are subsetted. Another piece of code indicates how many models have been retained following this action. 

```{r mod-sel, message=FALSE, warning=FALSE, echo=TRUE}

```
Step 5: The nested models are removed (see methods for justification).

```{r in-nest, message=FALSE, warning=FALSE, echo=TRUE}

```

Step 6: The remaining then undergo model averaging to create a final coefficient table.

```{r in-mod-ave, message=FALSE, warning=FALSE, echo=TRUE}

```


The variables listed immediately above are highlighted as having the most significant impact on the response variable following model averaging. 

### Brooding Models

The process is then repeated in full with the brooding subset. 

Step 1: A linear model is created for initial analysis.

```{r bolk-brood, message=FALSE, warning=FALSE, echo=TRUE}

```

Step 2: The global model is created. Sex:pers interaction is included in the brooding model. 

```{r coml-brood, message=FALSE, warning=FALSE, echo=TRUE}

```

Step 3: The model set is created and the age and partner age variables are linked to their quadratic counterparts. As the brooding data set is larger than the incubation data set, this action may take much longer than in the previous section. 

```{r br-dredge, message=FALSE, warning=FALSE, echo=TRUE, results = 'hide'}

```

Step 4: Model selection of those delta < 2

```{r br-sel2, message=FALSE, warning=FALSE, echo=TRUE}

```

Step 5: Nested models are removed.
```{r br-nest, message=FALSE, warning=FALSE, echo=TRUE}

```

Step 6: Final coefficient table created via model averaging. 
```{r br-ave, message=FALSE, warning=FALSE, echo=TRUE}

```

The variables listed above are highlighted as the most influential on the response variable following model averaging.