---
title: "Phylogenetic logistic regression modeling the probability of performing the broken-wing display as a function of ecological and life-history traits"
author: "L. de Framond"
date: "3/19/2021"
output: html_document
---


# Setup


```{r setup, include=TRUE}
knitr::opts_chunk$set(echo = TRUE, warning = F)

# necessary packages 
library(nlme)
library(ape)
library(geiger)
library(phylolm) # this is package that has the models we use

library(effects)
library(lattice)
library(ggeffects)
library(dplyr)

library(performance)
library(ggplot2)


# functions

#### function that converts names into Jetz names
convert_names = function(data, names){
  n_sp = length(data)
  newname = character()
  sapply(1:n_sp, function(x){
    line_name = which(names[,1]==data[x])
    if(length(line_name)>0){
      newname = names[line_name, 2]
    }else{
      newname = data[x]
    }
    return(newname)
  })
}

#### prepare table for data
modelselect = data.frame(Sample_size = NA, Model = NA, Alpha = NA, AIC = NA)
nmod = 1

modelselection = function(model){
  Sample_size = model$n
  Model = as.character(model$formula)[3]
  Alpha = model$alpha
  AIC = model$aic
  return(c(Sample_size, Model, Alpha, AIC))
}

# number of bootstrap replicate
rep = 100
```

### Loading data & tree.

```{r}


# load data and tree

tree_out<- read.tree("Stage2_Hackett_MCC_no_neg.tre")

dat <-read.csv("BWD_data_final_forDryad.csv", header=T)


```

## Models

### Models M0, lattitude and habitat

```{r}

###restrict dataset to species with mass and match with tree, then resort data. This needs to happen every time the data is subsetted to consider a new variable.

### model with mass, whole data
lat_dat <- subset(dat, dat$Lat != "NA")

lat_nc <- treedata(tree_out, lat_dat, sort=T)

str("Sample size of M0 :") ; length(lat_dat$species)

## make sure values match up
head(cbind(lat_nc$phy$tip.label, row.names(lat_dat), row.names(lat_nc$data)))

lat_dat <- lat_dat[lat_nc$phy$tip.label,]


# global model based on this sample
lat_phy <- phyloglm(performs ~ Lat + hab_scores + mobbing_recoded + colonial_recoded,
                     data=lat_dat,phy=lat_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                    boot = rep)
summary(lat_phy)

modelselect[nmod,] = modelselection(lat_phy)
nmod = nmod+1


# colonial_recoded is uniformative, remove this predictor from the model
l1_phy <- phyloglm(performs ~ Lat + hab_scores + mobbing_recoded,
                   data=lat_dat,phy=lat_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                   boot = rep)
summary(l1_phy)

modelselect[nmod,] = modelselection(l1_phy)
nmod = nmod+1


```

### Models M1, with body mass

```{r}


### model with mass
mass_dat <- subset(lat_dat, lat_dat$g_mean_dunning != "NA")
mass_nc <- treedata(tree_out, mass_dat, sort=T)

str("Sample size of M1 :") ; length(mass_dat$species)

## make sure values match up
head(cbind(mass_nc$phy$tip.label, row.names(mass_dat), row.names(mass_nc$data)))

mass_dat <- mass_dat[mass_nc$phy$tip.label,]

# make mass a numeric variable
mass_dat$g_mean_dunning <- as.numeric(mass_dat$g_mean_dunning)

# global model based on this sample
mass_phy <- phyloglm(performs ~ log(g_mean_dunning) + Lat + hab_scores + mobbing_recoded,
                     data=mass_dat,phy=mass_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL,
                     boot = rep)
summary(mass_phy)

modelselect[nmod,] = modelselection(mass_phy)
nmod = nmod+1


# log of mass is uniformative, remove

m1_phy <- phyloglm(performs ~ Lat + hab_scores + mobbing_recoded ,
                   data=mass_dat,phy=mass_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                   boot = rep)
summary(m1_phy)

modelselect[nmod,] = modelselection(m1_phy)
nmod = nmod+1

# mobbing behaviour is uniformative, remove

m2_phy <- phyloglm(performs ~ Lat + hab_scores ,
                   data=mass_dat,phy=mass_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                   boot = rep)
summary(m2_phy)

modelselect[nmod,] = modelselection(m2_phy)
nmod = nmod+1

```

### Models M2, with nest cover and nest location

```{r} 

# subset the dataset and keep only the species for which data for all predictors is available
nest_dat <- subset(lat_dat, lat_dat$nest_cover != "NA")
which(is.na(nest_dat$nestL))
which(is.na(nest_dat$nest_type))
nest_dat <- subset(nest_dat, nest_dat$nest_type != "NA")
length(nest_dat$nest_cover)

nest_nc <- treedata(tree_out, nest_dat, sort=T)

str("Sample size of M2 :") ; length(nest_dat$species)


## make sure values match up
head(cbind(nest_nc$phy$tip.label, row.names(nest_dat), row.names(nest_nc$data)))

nest_dat <- nest_dat[nest_nc$phy$tip.label,]


# global model based on this sample

nest_phy <- phyloglm(performs ~ Lat + hab_scores + nestL + nest_cover + nest_protec + nest_conspic,
                     data=nest_dat, phy=nest_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                     boot = rep)
summary(nest_phy)
check_collinearity(nest_phy)
modelselect[nmod,] = modelselection(nest_phy)
nmod = nmod+1

# remove nest protection

n1_phy <- phyloglm(performs ~ Lat + hab_scores + nestL + nest_cover + nest_conspic,
                     data=nest_dat, phy=nest_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                     boot = rep)
summary(n1_phy)
check_collinearity(n1_phy)
modelselect[nmod,] = modelselection(n1_phy)
nmod = nmod+1

# remove nest conspicuousness

n2_phy <- phyloglm(performs ~ Lat + hab_scores + nestL + nest_cover ,
                     data=nest_dat, phy=nest_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                     boot = rep)
summary(n2_phy)
check_collinearity(n2_phy)
modelselect[nmod,] = modelselection(n2_phy)
nmod = nmod+1



# nestL (location) is uniformative, remove

n3_phy <- phyloglm(performs ~ Lat + hab_scores + nest_cover,
                   data=nest_dat, phy=nest_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                   boot = rep)
summary(n3_phy)
modelselect[nmod,] = modelselection(n3_phy)
nmod = nmod+1

```


In this models, we included the predictors that were significant in M1 (lattitude, habitat score) and predictores related to the nest location (nest location and nest cover)
In the final model, only lattitude, habitat scores and nest cover had an influence on broken wing display.

### Models M3, with clutch size

```{r}


# subset the dataset
head(nest_dat$neggs_analys)

negg_dat <- subset(nest_dat, nest_dat$neggs_analys != "NA")
length(negg_dat$neggs_analys)
negg_nc <- treedata(tree_out, negg_dat, sort=T)

## make sure values match up
head(cbind(negg_nc$phy$tip.label, row.names(negg_dat), row.names(negg_nc$data)))

negg_dat <- negg_dat[negg_nc$phy$tip.label,]

str("Sample size of M3 :") ; length(negg_dat$species)


# global model based on this sample
negg_phy <- phyloglm(performs ~ Lat + hab_scores + nest_cover + neggs_analys,
                     data=negg_dat, phy=negg_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL,
                     boot = rep)
summary(negg_phy)
modelselect[nmod,] = modelselection(negg_phy)
nmod = nmod+1

# number of egg per nest is uniformative, remove

eg1_phy <- phyloglm(performs ~ Lat + hab_scores + nest_cover,
                    data=negg_dat,phy=negg_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL,
                    boot = rep)
summary(eg1_phy)
modelselect[nmod,] = modelselection(eg1_phy)
nmod = nmod+1

```

In this model, we included the number of eggs per clutch as a predictor, which we realised had no influence on broken wing display. Therefore, we do not consider this parameter in the next analysis, and all species for which we do not have data about clutch size are included back in the dataset. 

### Model M4 with precociality

```{r}

# subset data
head(nest_dat$chick_type_analys)

chick_dat <- subset(nest_dat, nest_dat$chick_type_analys != "NA")
length(chick_dat$chick_type_analys)
chick_nc <- treedata(tree_out, chick_dat, sort=T)

## make sure values match up
head(cbind(chick_nc$phy$tip.label, row.names(chick_dat), row.names(chick_nc$data)))

chick_dat <- chick_dat[chick_nc$phy$tip.label,]
str("Sample size of M4 :") ; length(chick_dat$species)

# global model based on this sample
chick_phy <- phyloglm(performs ~ Lat + hab_scores + nest_cover + chick_type_analys,
                      data=chick_dat, phy=chick_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL,
                      boot = rep)
summary(chick_phy)
check_collinearity(chick_phy)
modelselect[nmod,] = modelselection(chick_phy)
nmod = nmod+1

# lattitude is uniformative, remove

c1_phy <- phyloglm(performs ~  hab_scores + nest_cover + chick_type_analys,
                      data=chick_dat, phy=chick_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL,
                   boot = rep)
summary(c1_phy)
modelselect[nmod,] = modelselection(c1_phy)
nmod = nmod+1

```

In this model, we include the chick precociality in the predictors. With this model, chick precociality is informative, and lattitude becomes uninformative. Therefore, we drop it from the model from now on.

### Model M5 with incubation duty

```{r}

# subset data

head(nest_dat$inc_duty)

duty_dat <- subset(chick_dat, chick_dat$inc_duty != "NA")
length(duty_dat$inc_duty)
duty_nc <- treedata(tree_out, duty_dat, sort=T)

## make sure values match up
head(cbind(duty_nc$phy$tip.label, row.names(duty_dat), row.names(duty_nc$data)))

duty_dat <- duty_dat[duty_nc$phy$tip.label,]
str("Sample size of M5 :") ; length(duty_dat$species)


# global model based on this sample
duty_phy <- phyloglm(performs ~ hab_scores + nest_cover +  chick_type_analys + duty,
                     data=duty_dat,phy=duty_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                     boot = rep)
summary(duty_phy)
check_collinearity(duty_phy)
modelselect[nmod,] = modelselection(duty_phy)
nmod = nmod+1

# global model based on this sample
d1_phy <- phyloglm(performs ~ hab_scores + nest_cover +   duty,
                     data=duty_dat,phy=duty_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                     boot = rep)
summary(d1_phy)
check_collinearity(d1_phy)
modelselect[nmod,] = modelselection(d1_phy)
nmod = nmod+1


```

In these models, we included the incubation duty (whether the nest in taken care of by two parents (1) or just one parent (0) ). This parameter is significant, and chick type is not significant anymore. Therefore, all the species for which we don't have data about chick type are included back into the anlysis

### Model M6 with incubation duration

```{r}

# subset dataset

head(nest_dat$inc_mean)

inc_dat <- subset(nest_dat, nest_dat$inc_mean != "NA")

which(is.na(inc_dat$duty))
if(length(which(is.na(inc_dat$duty)))>0){
  inc_dat <- subset(inc_dat, inc_dat$duty != "NA")
}

length(inc_dat$inc_mean)
inc_nc <- treedata(tree_out, inc_dat, sort=T)

str("Sample size of M6 :") ; length(inc_dat$species)

## make sure values match up
head(cbind(inc_nc$phy$tip.label, row.names(inc_dat), row.names(inc_nc$data)))

inc_dat <- inc_dat[inc_nc$phy$tip.label,]


# global model based on this sample
inc_phy <- phyloglm(performs ~ hab_scores + nest_cover + duty + inc_mean,
                    data=inc_dat, phy=inc_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                    boot = rep)
summary(inc_phy)
modelselect[nmod,] = modelselection(inc_phy)
nmod = nmod+1

# remove duty
i2_phy <- phyloglm(performs ~ hab_scores + nest_cover + inc_mean,
                    data=inc_dat, phy=inc_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                    boot = rep)
summary(i2_phy)
modelselect[nmod,] = modelselection(i2_phy)
nmod = nmod+1

# To be sure that the species with very long incubation duration do not have a too important leverage effect on the model result, we run the sme model without the species with incubation durations longer than 40 days
outlierfree_incdat = subset(inc_dat, inc_dat$inc_mean < 40)
inc_nc <- treedata(tree_out, inc_dat, sort=T)

outlierfree_i2_phy <- phyloglm(performs ~ hab_scores + nest_cover + inc_mean,
                    data=inc_dat, phy=inc_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                    boot = rep)
saveRDS(outlierfree_i2_phy, "outlierfree_model_incubationdur_boot100.RDS")

```


### Model M7 with mass of eggs, mass of clutch and relative  mass of clutch


```{r}

# subset the dataset
head(nest_dat$inc_mean)

megg_dat <- subset(nest_dat, nest_dat$clutch_percent != "NA")
which(is.na(megg_dat$mass_egg_analys))
which(is.na(megg_dat$tot_mass_clutch))
which(is.na(megg_dat$inc_mean))
length(megg_dat$clutch_percent)
megg_dat <- subset(megg_dat, megg_dat$inc_mean != "NA")
megg_nc <- treedata(tree_out, megg_dat, sort=T)

## make sure values match up
head(cbind(megg_nc$phy$tip.label, row.names(megg_dat), row.names(megg_nc$data)))

megg_dat <- megg_dat[megg_nc$phy$tip.label,]
str("Sample size of M7 :") ; length(megg_dat$species)


# global model based on this sample
megg_phy <- phyloglm(performs ~ hab_scores + nest_cover + inc_mean +  mass_egg_analys + tot_mass_clutch + clutch_percent,
                    data=megg_dat,phy=megg_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                    boot = rep)
summary(megg_phy)
modelselect[nmod,] = modelselection(megg_phy)
nmod = nmod+1
check_collinearity(megg_phy)

# remove mass of egg
me2_phy <- phyloglm(performs ~ hab_scores + nest_cover + inc_mean +  tot_mass_clutch + clutch_percent,
                     data=megg_dat,phy=megg_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL,
                    boot = rep)
check_collinearity(me2_phy)
summary(me2_phy)
modelselect[nmod,] = modelselection(me2_phy)
nmod = nmod+1

# remove mass of clutch
me3_phy <- phyloglm(performs ~ hab_scores + nest_cover  +inc_mean + clutch_percent,
                    data=megg_dat,phy=megg_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                    boot = rep)
check_collinearity(me3_phy)
summary(me3_phy)
modelselect[nmod,] = modelselection(me3_phy)
nmod = nmod+1

# remove mass of clutch relative to the birds body mass
me4_phy <- phyloglm(performs ~ hab_scores + nest_cover + inc_mean   ,
                    data=megg_dat,phy=megg_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                    boot = rep)
check_collinearity(me4_phy)
summary(me4_phy)
modelselect[nmod,] = modelselection(me4_phy)
nmod = nmod+1

# remove mass of clutch relative to the birds body mass
me5_phy <- phyloglm(performs ~ hab_scores + nest_cover  ,
                    data=megg_dat,phy=megg_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                    boot = rep)
check_collinearity(me5_phy)
summary(me5_phy)
modelselect[nmod,] = modelselection(me5_phy)
nmod = nmod+1

```

In this model we test the effect of clutch mass (egg mass, total clutch mass and total clutch mass relative to the body mass of the bird). In this subset of the data, none of these parameters are informative.

### Model M8 with double brooding


```{r}

# subset the dataset
head(nest_dat$multibrood)

brood_dat <- subset(nest_dat, nest_dat$multibrood != "NA")
length(brood_dat$multibrood)
brood_nc <- treedata(tree_out, brood_dat, sort=T)


str("Sample size of M8 :") ; length(brood_dat$species)

## make sure values match up
head(cbind(brood_nc$phy$tip.label, row.names(brood_dat), row.names(brood_nc$data)))

brood_dat <- brood_dat[brood_nc$phy$tip.label,]


# global model based on this sample
brood_phy <- phyloglm(performs ~ hab_scores + nest_cover + multibrood,
                     data=brood_dat, phy=brood_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                     boot = rep)
summary(brood_phy)
modelselect[nmod,] = modelselection(brood_phy)
nmod = nmod+1

```

This model test if the number of broods raised per season impacts the broken-wing display. It does, so we keep it  for the next step of the analysis.

### Models M9 with longevity

```{r}
# subset dataset
head(brood_dat$longevity_y)

life_dat <- subset(brood_dat, brood_dat$longevity_y != "NA")
length(life_dat$longevity_y)
life_nc <- treedata(tree_out, life_dat, sort=T)

str("Sample size of M9 :") ; length(life_dat$species)

## make sure values match up
head(cbind(life_nc$phy$tip.label, row.names(life_dat), row.names(life_nc$data)))

life_dat <- life_dat[life_nc$phy$tip.label,]

# global model based on this sample
life_phy <- phyloglm(performs ~ hab_scores + nest_cover + multibrood + longevity_y,
                      data=life_dat, phy=life_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                     boot = rep)
summary(life_phy)
modelselect[nmod,] = modelselection(life_phy)
nmod = nmod+1


# lifespan not informative
li2_phy <- phyloglm(performs ~ hab_scores + nest_cover + multibrood,
                     data=life_dat, phy=life_nc$phy, method = "logistic_IG10", btol=30, start.beta = NULL, 
                   boot = rep)
summary(li2_phy)
modelselect[nmod,] = modelselection(li2_phy)
nmod = nmod+1

```


This model examines the effect of the maximum longevity on the broken wing display. It is not informative. I stopped there because the number of species dropped under 200 and we decided to stop at this point to avoid that the thinning of dataset obscures some previously shown significant effects. 

### Write down the summary table

```{r}

write.csv(modelselect, "ModelSelectionTable_btol30_boot100_2021-03-19.csv")

```

### Save all results and datasets for ploting

```{r}

goodmodel = list(lat_phy, l1_phy, mass_phy, m1_phy, m2_phy, nest_phy, n1_phy, n2_phy, n3_phy, negg_phy, eg1_phy, chick_phy, c1_phy, duty_phy, d1_phy, inc_phy,i2_phy, megg_phy, me2_phy, me3_phy, me4_phy, me5_phy, brood_phy, life_phy, li2_phy, dat, lat_dat, mass_dat, nest_dat, negg_dat, chick_dat, duty_dat, inc_dat, megg_dat, brood_dat, life_dat)
names(goodmodel) = c("M0_Lat-habscores-mobbing-colonial",
                     "M0_Lat-habscores-mobbing", 
                     "M1_Mass-Lat-mobbing-habscores",
                     "M1_Lat-mobbing-habscores",
                     "M1_Lat-habscores",
                     "M2_Lat-habscores-nestcover-nestloc-nestconspic-nestprotec",
                     "M2_Lat-habscores-nestcover-nestloc-nestconspic",
                     "M2_Lat-habscores-nestcover-nestloc", 
                     "M2_Lat-habscores-nestcover", 
                     "M3_Lat-habscores-nestcover-negg", 
                     "M3_Lat-habscores-nestcover", 
                     "M4_Lat-habscores-nestcover-precocial", 
                     "M4_habscores-nestcover-precocial", 
                     "M5_habscores-nestcover-precocial-duty",
                     "M5_habscores-nestcover-duty", 
                     "M6_habscores-nestcover-duty-incub", 
                     "M6_habscores-nestcover-incub", 
                     "M7_habscores-nestcover-incub-megg-mclutch-percentclutch", 
                     "M7_habscores-nestcover-incub-mclutch-clutchpercent",
                     "M7_habscores-nestcover-incub-clutchpercent",
                     "M7_habscores-nestcover-incub",
                     "M7_habscores-nestcover",
                     "M8_habscores-nestcover-multibrood", 
                     "M9_habscores-nestcover-multibrood-lifespan",
                     "M9_habscores-nestcover-multibrood", 
                     "dat",
                     "lat_dat",
                     "mass_dat",
                     "nest_dat", 
                     "negg_dat", 
                     "chick_dat", 
                     "duty_dat", 
                     "inc_dat", 
                     "megg_dat", 
                     "brood_dat", 
                     "life_dat")

saveRDS(goodmodel, "finalmodels_boot100_2021-03-19.RDS")

# good = readRDS("finalmodels.RDS")

```


