---
title: "Exploration of the Broken-wing-display dataset"
author: "L. de Framond"
date: "12/02/2021"
output:
  html_document:
    toc: true
---
\  
\  

## Introduction


This piece of code is a first insight in the dataset. It provides an overview of the sample size and distribution of each variable. Variables that have been gathered in different formats are put together. The outputs of this codes are:
1- an html file summarizing all the distribution plots for each variable
2- the final dataset that is then used in the models

\  
\  

## I - Exploration of data, recoding of variables, and calculation of egg mass


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = F, warning = F, fig.width = 5, fig.height = 3)

# Useful packages

library(ggplot2)
library(knitr)
library(dplyr)
ts=12

# Useful functions


calc_egg_mass = function(length, breadth){
  # egg length and largest diameter in mm, returns mass in g
  0.548*(length/10)*((breadth/10)^2)
}

calc_egg_mass(42,30)



complete_data = function(avrg, min, max){
  # function that completes the column with average data with the mean of the range data when available
  
  eggs = data.frame(avrg = is.na(avrg), max = is.na(max))
  
  sapply(1:length(avrg), function(x){
  
  if(eggs$avrg[x]==TRUE){
    if(eggs$max[x]==TRUE){
      dat = NA
    }else{
      dat = (max[x]+ min[x])/2
    }
  }else{
    dat = avrg[x]
  }
  
  return(dat)
})
}


# read raw data

bwd_data = read.csv("BWD_data_raw_forDryad.csv")
```

\ 

### 1 - Habitats


```{r, echo=FALSE}
ggplot(data = bwd_data, aes(x=hab_scores))+
  theme_classic()+
  theme(text=element_text(size = ts))+
  xlab("Habitat score (0 = desert,  8 = dense forest)")+
  ylab("number of species")+
  geom_histogram(bins = 25, fill  ="darkblue")
```

**Fig. 1: Distribution of the habitats scores, determined using the occurrence of words in the habitat description**

\ 


### 2 - Number of clutches, replacement clutches

\ 

```{r}
# count the number and percentage of species for which we have the data about the brood number
nas_broods = sum(is.na(bwd_data$nbrood_max))
brood_percent = round(100*(570-nas_broods)/570)
```

We have data about the maximum number of brood raised every year for `r brood_percent` % of the species in the dataset

\ 

```{r}

ggplot(bwd_data, aes(x=nbrood_max))+
  theme_classic()+
  theme(text=element_text(size = ts))+
  xlab("Maximum number of brood raised per year")+
  ylab("number of species")+
  geom_histogram( stat = "count", fill  ="darkblue")

```

**Fig. 2: Distribution of the maximum number of broods raised per year**

\ 

### 3 - Number of eggs per clutch

```{r}
# complete the "typical number of egg" column with the average of the range, and then counts the number of species for which we have the number of eggs

bwd_data$neggs_analys = complete_data(avrg = bwd_data$neggs_typ, 
                                      min = bwd_data$neggs_min, 
                                      max = bwd_data$neggs_max)

nas_neggs = sum(is.na(bwd_data$neggs_analys))
neggs_percent = round(100*(570-nas_neggs)/570)

```

We have information about the average number of eggs laid per clutch for `r neggs_percent` % of the species (including species for which we could only collect the maximum and minimum egg number and species for which I found the average number of eggs per nest)

\ 

```{r}


ggplot(bwd_data, aes(x=as.factor(round(neggs_analys))))+
  theme_classic()+
  theme(text=element_text(size = ts))+
  xlab("Typical number of eggs per brood")+
  ylab("number of species")+
  geom_histogram( stat = "count", fill  ="darkblue")

```

**Fig. 3: Distribution of Typical number of eggs per clutch **

\ 


```{r}


dataforplot = data.frame(egg_pres = bwd_data$neggs_analys, neggs_max = bwd_data$neggs_max, neggs_min = bwd_data$neggs_min,  sp = bwd_data$species)
dataforplot = arrange(dataforplot, egg_pres)
dataforplot$index = c(1:(570-nas_neggs),rep(NA, nas_neggs))

ggplot(dataforplot, aes(y= egg_pres, x = index))+
  theme_classic()+
  theme(text=element_text(size = ts), axis.text.x = element_blank())+
  xlab("Sorted species")+
  ylab("Number of eggs laid per clutch \n (dots = average, grey bars = range)")+
  geom_errorbar(aes(ymin = neggs_min, ymax = neggs_max), col = "grey40", alpha = .5)+
  geom_point(col = "darkblue", na.rm=TRUE)
  

  
```

**Fig. 4: Number of eggs per clutch (average and range) for all species in the dataset, sorted from smallest to largest clutch**

\ 

### 4 - Egg dimensions and mass

```{r}
# complete average data and counts species with data

bwd_data$length_egg_analys = complete_data(avrg = bwd_data$length_egg_avrg_mm,
                                           min = bwd_data$length_egg_min_mm,
                                           max = bwd_data$length_egg_max_mm)

bwd_data$breadth_egg_analys = complete_data(avrg = bwd_data$breadth_egg_avrg_mm,
                                           min = bwd_data$breadth_egg_min_mm,
                                           max = bwd_data$breadth_egg_max_mm)


nas_dim = sum(is.na(bwd_data$length_egg_analys))
dim_percent = round(100*(570-nas_dim)/570)

# invert breadth and length for species where I have typed it wrong (in general length > breadth)
inverted = which(bwd_data$breadth_egg_analys>bwd_data$length_egg_analys)

for(i in 1:length(inverted)){
  len = bwd_data$breadth_egg_analys[inverted[i]]
  bwd_data$breadth_egg_analys[inverted[i]] = bwd_data$length_egg_analys[inverted[i]]
  bwd_data$length_egg_analys[inverted[i]] = len
}

# calculate egg mass with the formula from Hoyt  (1979)

bwd_data$mass_egg_calc = calc_egg_mass(bwd_data$length_egg_analys, bwd_data$breadth_egg_analys)

# complete egg mass data with midpoint of the range of egg masses

bwd_data$mass_egg_compl = complete_data(avrg = bwd_data$mass_egg_avrg_g,
                                           min = bwd_data$mass_egg_min_g,
                                           max = bwd_data$mass_egg_max_g)

nas_mass = sum(is.na(bwd_data$mass_egg_compl))
mass_percent = round(100*(570-nas_mass)/570)

nspfem = 570-sum(is.na(bwd_data$mass_percent_female))
percentmassfem = 100*(nspfem/570)


# complete egg mass with calculated egg mass 

bwd_data$mass_egg_analys = sapply(1:length(bwd_data$mass_egg_avrg_g), function(x){
  eggs = is.na(bwd_data$mass_egg_compl)
  if(eggs[x]==TRUE){
    dat = bwd_data$mass_egg_calc[x]
  }else{
    dat = bwd_data$mass_egg_compl[x]
  }
  return(dat)
})

nas2 = sum(is.na(bwd_data$mass_egg_analys))

```

We have information about the egg dimensions for `r dim_percent` % of the species, and the egg mass for `r mass_percent` % of the species (the egg mass data can be completed for the species which egg dimension we have, with Hoyts formula (1979)). In addition, we have the egg mass relative to the females body mass for `r nspfem` species (`r round(percentmassfem)` % of all species)

\ 

```{r}
dataforplot = data.frame(egg_pres = bwd_data$mass_egg_analys, neggs_max = bwd_data$mass_egg_max_g, neggs_min = bwd_data$mass_egg_min_g,  sp = bwd_data$species, calc = bwd_data$mass_egg_calc, meas= bwd_data$mass_egg_compl)
dataforplot = arrange(dataforplot, egg_pres)
dataforplot$index = c(1:(570-nas2),rep(NA, nas2))

ggplot(dataforplot, aes(y= egg_pres, x = index))+
  theme_classic()+
  theme(text=element_text(size = ts), axis.text.x = element_blank())+
  xlab("Sorted species")+
  ylab("Egg mass (g)")+
  geom_text(aes(x = 10, y = 160, label = "A"))+
  geom_errorbar(aes(ymin = neggs_min, ymax = neggs_max), col = "grey40", alpha = .5)+
  geom_point(col = "darkblue")+
  geom_point(col = "red", alpha=.2, aes(y = calc, x=index))

  
```

**Fig. 5: Egg mass (average and range), sorted from the smallest to largest egg. ** The blue dots are the data that figure in the final dataset. The red dots are the values calculated with Hoyt's formula. Grey bars: range of egg mass (when available). 

\  

```{r}

summary(lm(bwd_data$mass_egg_calc~bwd_data$mass_egg_avrg_g))

ggplot(bwd_data, aes(y= mass_egg_calc, x = mass_egg_avrg_g))+
  theme_classic()+
  theme(text=element_text(size = ts))+
  xlab("Measured egg mass (g)")+
  ylab("Calculated egg mass (g)\n Mass = .548 x Length x Breadth ²")+
  geom_abline(slope=1, intercept = 0, col = "red")+
  geom_point(col = "darkblue",na.rm=TRUE)

  
```

**Fig. 6: Comparison between the calculated egg mass (with Hoyt's formula, Hoyt 1979) and the measured egg mass** 

The egg mass calculated from egg dimensions is very similar to the egg mass measured, therefore, completing the missing data points with the calculated values is reasonable.

```{r}

# calculate total mass of clutch and percet of body mass
bwd_data$tot_mass_clutch  = bwd_data$neggs_analys * bwd_data$mass_egg_analys
bwd_data$clutch_percent = 100*bwd_data$tot_mass_clutch / bwd_data$g_mean_dunning


nas_clutm = sum(is.na(bwd_data$clutch_percent))
clutm_percent = round(100*((570-nas_clutm)/570))

```

I calculated the total mass of a clutch, and the mass of a clutch relative to the body mass of the bird (we have data for `r clutm_percent` % of species in the dataset.)


```{r}


ggplot(bwd_data, aes(x=clutch_percent))+
  theme_classic()+
  theme(text=element_text(size = ts))+
  xlab("Total mass of a clutch relative to the bird`s body mass (%)")+
  ylab("Number of species")+
  geom_histogram(fill  ="darkblue", bins = 45)

```

**Fig. 7: Distribution of the average mass of a total clutch relative to the body size of the female** 

\ 

### 5 - Longevity

```{r}
# count species for which we have longevity

nas_life = sum(is.na(bwd_data$longevity_y))
life_percent = round(100*((570-nas_life)/570))

```

We have information about the longevity of `r 570-nas_life` species of the dataset (`r life_percent` % of all species in the dataset). 

```{r}

ggplot(data.frame(bwd_data), aes(x=round(longevity_y)))+
  theme_classic()+
  theme(text=element_text(size = ts))+
  xlab("Longevity (years)")+
  ylab("Number of species")+
  geom_histogram( stat = "count", fill  ="darkblue")


```

**Fig. 8: Distribution of longevity in the dataset**

\ 

### 6 - Colonial species

```{r}
nas_col = sum(is.na(bwd_data$colonial_YN))
col_percent = round(100*((570-nas_col)/570))
bwd_data$colonial_YN = as.factor(bwd_data$colonial_YN)

```

In the raw data (before recoding), I have information on coloniality for `r col_percent` % of the species.

I recoded the description of coloniality in the following way:
*- always colonial  = 2
*- semi-colonial, sometimes colonial, loose groups = 1
*- solitary, no description, sometimes in loose groups or colonies = 0.

See the following table summarizing the codes I used for each description:

\ 

```{r}

colony = data.frame(phrase = levels(bwd_data$colonial_YN), code = c(0, 0, 2, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0))
kable(colony)
 x=3

bwd_data$colonial_recoded = sapply(1:length(bwd_data$colonial_YN), function(x){
  if(is.na(bwd_data$colonial_YN[x])){
    dat = 0
  }else{
    phrase = which(colony$phrase==bwd_data$colonial_YN[x])
    dat = colony[phrase, "code"]
  }
  return(dat)
})

ggplot(bwd_data, aes(x=as.factor(colonial_recoded)))+
  theme_classic()+
  theme(text=element_text(size = ts))+
  xlab("Coloniality")+
  ylab("number of species")+
  geom_histogram( stat = "count", fill  ="darkblue")

``` 


**Fig. 9: Distribution of colonial, semi colonial and solitary species**

\ 

### 7 - Mobbing behaviour


```{r}
nas_mob = sum(is.na(bwd_data$mobbing_YN))
mob_percent = round(100*((570-nas_mob)/570))
bwd_data$mobbing_YN = as.factor(bwd_data$mobbing_YN)

```

Before recoding, I have information about mobbing behaviour for `r 570-mob_percent` % of the species. We assume that the default value therefore was 0 (no mobbing) if noting is mentioned. 

Here after, you can find the list of descriptions I copied when it was not mentioned directly "mobbing":

Of these, only "alarm calls", "mimics calls of predators", and "very vocal" were considered as "not mobbing". 

```{r}

mob = data.frame(phrase = levels(bwd_data$mobbing_YN), code = c(0, 0, 1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,0,1))
kable(mob)

bwd_data$mobbing_recoded = sapply(1:length(bwd_data$mobbing_YN), function(x){
  if(is.na(bwd_data$mobbing_YN[x])){
    dat = 0
  }else{
    phrase = which(mob$phrase==bwd_data$mobbing_YN[x])
    dat = mob[phrase, "code"]
  }
  return(dat)
})

ggplot(bwd_data, aes(x=as.factor(mobbing_recoded)))+
  theme_classic()+
  theme(text=element_text(size = ts))+
  xlab("Mobbing behaviour")+
  ylab("number of species")+
  geom_histogram( stat = "count", fill  ="darkblue")
```

**Fig. 10: Distribution of mobbing and non-mobbing species**

\ 

```{r}
### some last recoding 


### complete incubation duty mean with the midpoint between min and max if available

bwd_data$inc_length <- ifelse(is.na(bwd_data$inc_mean), mean(c(bwd_data$inc_max, bwd_data$inc_min)), bwd_data$inc_mean)
# plot(bwd_data$inc_mean~bwd_data$inc_min)

### recode single brooded vs multi brooded
bwd_data$multibrood <- ifelse(bwd_data$nbrood_max=="1", 0, 1)

### recode incubation duty in 0 and 1
bwd_data$duty = bwd_data$inc_duty-1

### make a new latitude variable reflecting maximum distance from equator, which could be minimum or maximum

bwd_data$Lat <- ifelse(abs(bwd_data$MinLat) > abs(bwd_data$MaxLat), abs(bwd_data$MinLat), abs(bwd_data$MaxLat))


```

## II - Data summary

\ 

In this section, you will find a brief summary of the sample size available for each variable:

```{r}
variables = c("Number of broods", "egg number", "egg mass", "clutch mass relative to adult mass", "longevity", "coloniality", "nest defense")

data_type = c("numeric", "numeric", "grams", "%", "years", "index (0/1/2)", "index (0/1)")

n_species = 570-c(nas_broods, nas_neggs, nas_mass, nas_clutm, nas_life, nas_col, nas_mob)

percent_species = c(brood_percent, neggs_percent, mass_percent, clutm_percent, life_percent, col_percent, mob_percent)


bwd_summary = data.frame(Variable = variables, Type = data_type,  'Sample size' = n_species, 'Sample size (percent)' = percent_species)

kable(bwd_summary)

## This last line will save the current dataset if uncommented:

# write.csv(bwd_data, "BWD_data_final.csv")
```

\ 
