R codes for: Understanding chimpanzee nesting site selection through the lens of vegetation science1 Article published in American Journal of Primatology
12 February 2022
Abstract
Sleep is an important aspect of great ape life, and they build night nests every day to sleep. For chimpanzees, each group select a nesting site where each individual builds a nest, mostly in a tree. Some previous hypotheses include predation avoidance and thermoregulation, with measures on nest heights and tree heights. However, the precise influence of vegetation composition and structure on chimpanzee nesting site selection remains unknown, though they are known to prefer some particular tree species. Using botanical inventories around nesting sites, we show that the tree community structure drives chimpanzee nesting site selection. It was previously thought that habitat types or their abundance were the determinant factors of nesting site selection in chimpanzees. Results from this study indicate that the role played by habitats in nesting site selection is dependent on the botanical characteristics of their different clusters, including the variation in the average tree heights, the abundance of all trees, the abundance of nesting trees and the preferred nesting tree species at the habitat level. Furthermore, the selection of a nesting site by chimpanzees may be linked to the presence of their preferred nesting tree species, but not necessarily to particular vegetation composition. Also, the importance of tree characteristics (DBH and height) in chimpanzee nesting site selection is not tree specific but may characterize the general pattern of the site. Our results demonstrate how considering vegetation parameters can help understand the importance of sleep in great ape lives. We anticipated our study to be a starting point for studying spatial cognition in chimpanzees through their nesting ecology. For example, testing how chimpanzees use spatial memory to find nesting sites may be of great interest.
Contents
- 1 Introduction
- 2 Tree preference for nest building
- 3 Modelling nesting site selection (Logistic regression)
- 3.1 Data sets
- 3.2 Correlation between variables
- 3.3 Constructing models
- 3.4 Determining the confidence intervals of the coefficients, and Odd ratios
- 3.5 Overall effects of habitat and the interactions
- 3.6 Calculate the predicted values
- 3.7 Predicted values and confidence limits into probabilities
- 3.8 Plotting the GLM results
- 4 Tree DBH preference for chimpanzee nesting
- Session info
“It always seems impossible until it’s done.”
—Nelson Mandela
1 Introduction
The file showing the codes used in the article published in American Journal of Primatology
1.1 Loading packages
The codes bellow will load the required packages used for the analyses. The
function ipak is set to install the packages if they are not yet installed
and then load them.Note that some of those packages were used for the present
project.
To know which package was used for a particular function, simply press F2 on
that function and the details will open in a new windows.
package_list=c('dendextend','ggdendro','ggplot2','ggpubr','ggthemes',
'indicspecies','MuMIn','NbClust','plyr','png','questionr',
'adehabitatHS', 'knitr','aod','car','clubSandwich','cluster',
'corrr','cowplot','DescTools','factoextra','fpc',
'ggfortify','iNEXT','interactions','jtools',
'openxlsx','optpart','phyloseq','pvclust','reshape2',
'rsq','vegan','viridis','xlsxjars','rms','ggrepel',
'ExcelFunctionsR','rJava', 'tidyverse', 'sf','here',
'tmaptools','ggmap')
ipak <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if (length(new.pkg))
install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
ipak(package_list)## dendextend ggdendro ggplot2 ggpubr ggthemes
## TRUE TRUE TRUE TRUE TRUE
## indicspecies MuMIn NbClust plyr png
## TRUE TRUE TRUE TRUE TRUE
## questionr adehabitatHS knitr aod car
## TRUE TRUE TRUE TRUE TRUE
## clubSandwich cluster corrr cowplot DescTools
## TRUE TRUE TRUE TRUE TRUE
## factoextra fpc ggfortify iNEXT interactions
## TRUE TRUE TRUE TRUE TRUE
## jtools openxlsx optpart phyloseq pvclust
## TRUE TRUE TRUE TRUE TRUE
## reshape2 rsq vegan viridis xlsxjars
## TRUE TRUE TRUE TRUE TRUE
## rms ggrepel ExcelFunctionsR rJava tidyverse
## TRUE TRUE TRUE TRUE TRUE
## sf here tmaptools ggmap
## TRUE TRUE TRUE TRUE2 Tree preference for nest building
The codes bellow calculate the Manly indexes for each tree species used by
chimpanzees to build their nests. To create the data table for this analysis,
for each species and for each habitat type, we calculated the number of time
the tree species was used to build nest and we created a column for the
availability of the species, that is the abundance of that species in each
habitat type. The analysis was done considering all the tree species at once,
so that the tree species that were not used for nest building in a particular
habitat receive a score of 0 for that habitat. For this analysis, Av
indicates the availability and Chimp indicates the use for chimpanzee.
options(scipen = 1, digits = 2) #set to two decimal
## Tree_Preference
Tree_Preference=here("C:/Users/tedon/OneDrive - UGent/PHD Luc/Articles en cours/",
"Ecological factors of nesting site choice/Results/",
"Tree_Preference.xlsx")
Tree_Preference_df <- openxlsx::readWorkbook(xlsxFile = Tree_Preference,
sheet = "Tree_Preference")
df_ka=Tree_Preference_df[Tree_Preference_df$Av_KA>0,]#Remove species with null availability
df_kj=Tree_Preference_df[Tree_Preference_df$Av_KJ>0,]#Remove species with null availability
df_rap=Tree_Preference_df[Tree_Preference_df$Av_RAP>0,]#Remove species with null availability
df_rip=Tree_Preference_df[Tree_Preference_df$Av_RIP>0,]#Remove species with null availability
Av_KA=df_ka$Av_KA;names(Av_KA)=df_ka$Species
Av_KJ=df_kj$Av_KJ;names(Av_KJ)=df_kj$Species
Av_RAP=df_rap$Av_RAP;names(Av_RAP)=df_rap$Species
Av_RIP=df_rip$Av_RIP;names(Av_RIP)=df_rip$Species
Av_All=Tree_Preference_df$Av_All;names(Av_All)=Tree_Preference_df$Species
Chimp_KA=df_ka$Chimp_KA;names(Chimp_KA)=df_ka$Species
Chimp_KJ=df_kj$Chimp_KJ;names(Chimp_KJ)=df_kj$Species
Chimp_RAP=df_rap$Chimp_RAP;names(Chimp_RAP)=df_rap$Species
Chimp_RIP=df_rip$Chimp_RIP;names(Chimp_RIP)=df_rip$Species
Chimp_All=Tree_Preference_df$Chimp_All;names(Chimp_All)=Tree_Preference_df$Species
Gor_KJ=df_kj$Gor_KJ;names(Gor_KJ)=df_kj$Species
Gor_RAP=df_rap$Gor_RAP;names(Gor_RAP)=df_rap$Species
Gor_All=Tree_Preference_df$Gor_All;names(Gor_All)=Tree_Preference_df$Species
## Computation of wi (preference)
wi_Chimp_All <- widesI(Chimp_All,Av_All,alpha = 0.05);
wi_Chimp_All_df=cbind(data.frame(wi_Chimp_All$wi),
data.frame(wi_Chimp_All$se.wi),
data.frame(wi_Chimp_All$chisquwi));
colnames(wi_Chimp_All_df)=c("wi", "se.wi", "testwi","p");
wi_Chimp_All_df$Species=rownames(wi_Chimp_All_df);
wi_Chimp_All_df$Animal="Chimpanzee"
knitr::kable(wi_Chimp_All_df[wi_Chimp_All_df$wi>0,-c(5,6)])| wi | se.wi | testwi | p | |
|---|---|---|---|---|
| Antidesma laciniatum | 8.85 | 8.82 | 0.79 | 0.37 |
| Desbordesia glaucescens | 0.18 | 0.17 | 22.31 | 0.00 |
| Anthonotha macrophylla | 0.32 | 0.32 | 4.63 | 0.03 |
| Heisteria parvifolia | 27.41 | 1.50 | 310.14 | 0.00 |
| Irvingia gabonensis | 1.77 | 1.76 | 0.19 | 0.66 |
| Klainedoxa gabonensis | 0.75 | 0.75 | 0.11 | 0.73 |
| Lepidobotrys staudtii | 2.16 | 0.95 | 1.48 | 0.22 |
| Maesobotrya sp | 1.56 | 1.56 | 0.13 | 0.72 |
| Olax latifolia | 0.87 | 0.87 | 0.02 | 0.88 |
| Omphalocarpum elatum | 9.23 | 4.56 | 3.26 | 0.07 |
| Pausinystalia lane-poolei | 1.50 | 1.05 | 0.22 | 0.64 |
| Plagiostyles africana | 0.16 | 0.16 | 26.02 | 0.00 |
| Polyalthia suaveolens | 0.47 | 0.27 | 3.94 | 0.05 |
| Santiria trimera | 0.43 | 0.30 | 3.61 | 0.06 |
| Sorindeia grandifolia | 3.18 | 1.03 | 4.48 | 0.03 |
| Strombosiopsis tetranda | 0.62 | 0.62 | 0.36 | 0.55 |
| Trichilia gilgiana | 2.53 | 2.52 | 0.37 | 0.54 |
| Uapaca acuminata | 1.38 | 0.61 | 0.39 | 0.53 |
| Uapaca guineensis | 1.84 | 0.68 | 1.52 | 0.22 |
| Uapaca vanhouttei | 2.41 | 0.97 | 2.14 | 0.14 |
| Uvariastrum pierreanum | 4.08 | 4.07 | 0.57 | 0.45 |
3 Modelling nesting site selection (Logistic regression)
Following the example from this website: https://stats.idre.ucla.edu/r/dae/logit-regression/
3.1 Data sets
Here the package openxlsx is used to manage Excel files in R
options(scipen = 1, digits = 2) #set to two decimal
cluster_tree="C:/Users/tedon/OneDrive - UGent/PHD Luc/Articles en cours/Ecological factors of nesting site choice/Data/Botany/Analyses_Botany-Parcelles sans nids June 2021 (24-07-2021).xlsx"
## getSheetNames(cluster_tree)
Plot_Charact_All_For_GLM <- readWorkbook(xlsxFile = cluster_tree,
sheet = "Plot_Charact_All_For_GLM")
Plot_Charact_All_For_GLM$Liana_all=Plot_Charact_All_For_GLM$Liana_05_10+
Plot_Charact_All_For_GLM$Liana_10_more
## Re-coding variables to correspond to the codes in the manuscript
Plot_Charact_All_For_GLM$NP=Plot_Charact_All_For_GLM$Pres_Abs_Chimp
Plot_Charact_All_For_GLM$HAB=Plot_Charact_All_For_GLM$Hab_Eng
Plot_Charact_All_For_GLM$CTA=Plot_Charact_All_For_GLM$Count_Trees_All
Plot_Charact_All_For_GLM$AHA=Plot_Charact_All_For_GLM$Avrg_Height_All
Plot_Charact_All_For_GLM$CHC=Plot_Charact_All_For_GLM$Chimp_Hab_Com
Plot_Charact_All_For_GLM$CHS=Plot_Charact_All_For_GLM$Chimp_Hab_Select
Plot_Charact_All_For_GLM$CAC=Plot_Charact_All_For_GLM$Chimp_All_Com
Plot_Charact_All_For_GLM$CAS=Plot_Charact_All_For_GLM$Chimp_All_Select
Plot_Charact_All_For_GLM$ADT=Plot_Charact_All_For_GLM$Avrg_DBH_Trees
Plot_Charact_All_For_GLM$LA=Plot_Charact_All_For_GLM$Liana_all
Plot_Charact_All_For_GLM$C2x2=Plot_Charact_All_For_GLM$Count2x2SN_All
Plot_Charact_All_For_GLM$C4x4=Plot_Charact_All_For_GLM$Count4x4SN_All
## Selecting the columns used in the analysis
Chimp_GLM_Data=Plot_Charact_All_For_GLM[Plot_Charact_All_For_GLM$Chimpanzee=="Yes",
c("Type","Chimpanzee",
"NP","HAB", "CTA","AHA","CHC","CHS","CAC",
"CAS","ADT","LA","C2x2","C4x4" )]
knitr::kable(Chimp_GLM_Data[c(1:6),], format = "markdown")| Type | Chimpanzee | NP | HAB | CTA | AHA | CHC | CHS | CAC | CAS | ADT | LA | C2x2 | C4x4 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Without nest | Yes | 0 | OSF | 24 | 18 | 0 | 2 | 1 | 2 | 23 | 8 | 17 | 4 |
| Without nest | Yes | 0 | OSF | 19 | 27 | 2 | 3 | 3 | 5 | 31 | 4 | 15 | 8 |
| Without nest | Yes | 0 | OSF | 16 | 22 | 1 | 1 | 2 | 2 | 26 | 9 | 13 | 12 |
| Without nest | Yes | 0 | OSF | 27 | 22 | 1 | 2 | 6 | 2 | 26 | 2 | 9 | 13 |
| Without nest | Yes | 0 | OSF | 22 | 26 | 0 | 1 | 6 | 1 | 24 | 2 | 14 | 7 |
| Without nest | Yes | 0 | OSF | 14 | 28 | 0 | 2 | 4 | 2 | 26 | 2 | 8 | 14 |
3.2 Correlation between variables
The pearson correlation was considered
options(scipen = 1, digits = 2) #set to two decimal
Chimp_GLM_Data_Corr_Spear <- correlate(Chimp_GLM_Data[,-c(1:4)],method="spearman", quiet = TRUE)
Chimp_GLM_Data_Corr_Pears <- correlate(Chimp_GLM_Data[,-c(1:4)],method="pearson", quiet = TRUE)
Chimp_GLM_Data_Corr_Spear_df1=as.data.frame(Chimp_GLM_Data_Corr_Spear)
Chimp_GLM_Data_Corr_Pears_df1=as.data.frame(Chimp_GLM_Data_Corr_Pears)
Chimp_GLM_Data_Corr_Spear_df2=round(Chimp_GLM_Data_Corr_Spear_df1[,-1],digits=3)
Chimp_GLM_Data_Corr_Pears_df2=round(Chimp_GLM_Data_Corr_Pears_df1[,-1],digits=3)
Chimp_GLM_Data_Corr_Spear_df=cbind(Chimp_GLM_Data_Corr_Spear_df1$term,Chimp_GLM_Data_Corr_Spear_df2)
Chimp_GLM_Data_Corr_Pears_df=cbind(Chimp_GLM_Data_Corr_Pears_df1$term,Chimp_GLM_Data_Corr_Pears_df2)
colnames(Chimp_GLM_Data_Corr_Spear_df)=c("term","CTA","AHA","CHC","CHS","CAC","CAS","ADT",
"LA","C2x2","C4x4")
colnames(Chimp_GLM_Data_Corr_Pears_df)=c("term","CTA","AHA","CHC","CHS","CAC","CAS","ADT",
"LA","C2x2","C4x4")
## Saving correlation results in an excel file, pearson and spearman in different sheets
List_Chimp_GLM_Data_Corr <- list("Chimp_GLM_Data_Corr_Spear_df2"= Chimp_GLM_Data_Corr_Spear_df,
"Chimp_GLM_Data_Corr_Pears_df2"=Chimp_GLM_Data_Corr_Pears_df)
write.xlsx(List_Chimp_GLM_Data_Corr, file = "Chimp_GLM_Data_Corr_df24.xlsx")
## Printing correlation results
library(knitr)
knitr::kable(Chimp_GLM_Data_Corr_Pears_df, format = "markdown")| term | CTA | AHA | CHC | CHS | CAC | CAS | ADT | LA | C2x2 | C4x4 |
|---|---|---|---|---|---|---|---|---|---|---|
| CTA | NA | -0.05 | 0.36 | 0.21 | 0.45 | 0.41 | -0.18 | 0.22 | 0.22 | 0.31 |
| AHA | -0.05 | NA | -0.06 | 0.12 | 0.08 | 0.11 | 0.72 | 0.10 | 0.08 | 0.10 |
| CHC | 0.36 | -0.06 | NA | 0.22 | 0.50 | 0.38 | -0.06 | 0.03 | -0.02 | 0.03 |
| CHS | 0.21 | 0.12 | 0.22 | NA | 0.25 | 0.62 | 0.09 | -0.02 | -0.03 | -0.02 |
| CAC | 0.45 | 0.08 | 0.50 | 0.25 | NA | 0.24 | 0.00 | 0.09 | 0.05 | 0.13 |
| CAS | 0.41 | 0.11 | 0.38 | 0.62 | 0.24 | NA | 0.06 | 0.06 | 0.14 | 0.12 |
| ADT | -0.18 | 0.72 | -0.06 | 0.09 | 0.00 | 0.06 | NA | 0.14 | -0.01 | 0.00 |
| LA | 0.22 | 0.10 | 0.03 | -0.02 | 0.09 | 0.06 | 0.14 | NA | 0.13 | 0.17 |
| C2x2 | 0.22 | 0.08 | -0.02 | -0.03 | 0.05 | 0.14 | -0.01 | 0.13 | NA | 0.34 |
| C4x4 | 0.31 | 0.10 | 0.03 | -0.02 | 0.13 | 0.12 | 0.00 | 0.17 | 0.34 | NA |
Here we can see that none of the correlation was higher than 0.7.
3.3 Constructing models
options(scipen = 1, digits = 2) #set to two decimal
## Convert habitat variable into factors
Chimp_GLM_Data$HAB=factor(Chimp_GLM_Data$HAB)
CGData=Chimp_GLM_Data
## Constructing models
NP12=Chimp_logit_wHab <- glm(NP ~ CTA+AHA+CHC+CHS+CAC+CAS+LA+C2x2+C4x4+ADT,
data=CGData,family="binomial")
NP11=Chimp_logit_Sig2 <- glm(NP ~ CTA+AHA+CHC+CHS +C2x2+C4x4,
data=CGData,family="binomial")
NP09=Chimp_logit_Sig <- glm(NP ~ CTA+AHA+CHC+CHS,data=CGData,family="binomial")
NP10=Chimp_logit_Sig3 <- glm(NP ~ CTA+AHA+CHC+CHS+LA,data=CGData,family="binomial")
NP07=Chimp_logit_wHab2 <- glm(NP ~ CTA+AHA+CHC+CHS+CAC+CAS+LA+C2x2+C4x4+ADT+HAB,
data=CGData,family="binomial")
NP08=Chimp_logit_all_AHA <- glm(NP ~ CTA+AHA+CHC+CHS+CAC+CAS+ADT+HAB+HAB:AHA,
data=CGData,family="binomial")
NP06=Chimp_logit_all_CHS <- glm(NP ~ CTA+AHA+CHC+CHS+CAC+CAS+ADT+HAB+HAB:CHS,
data=CGData,family="binomial")
NP05=Chimp_logit_all_CTA <- glm(NP ~ CTA+AHA+CHC+CHS+CAC+CAS+ADT+HAB+HAB:CTA,
data=CGData,family="binomial")
NP03=Chimp_logit_all_CHC <- glm(NP ~ CTA+AHA+CHC+CHS+CAC+CAS+ADT+HAB+HAB:CHC,
data=CGData,family="binomial")
NP04=Chimp_logit_all_CHC_AHA <- glm(NP ~ CTA+AHA+CHC+CHS+CAC+CAS+ADT+HAB+HAB:CHC+HAB:AHA,
data=CGData,family="binomial")## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredNP02=Chimp_logit_all_CHC_CHS <- glm(NP ~ CTA+AHA+CHC+CHS+CAC+CAS+ADT+HAB+HAB:CHC+HAB:CHS,
data=CGData,family="binomial")
NP01=Chimp_logit_all_CHC_CTA <- glm(NP ~ CTA+AHA+CHC+CHS+CAC+CAS+ADT+HAB+HAB:CHC+HAB:CTA,
data=CGData,family="binomial")
NP07=Chimp_logit_all <- glm(NP ~ CTA+CHC+CHS+CAC+LA+CAS+C2x2+C4x4+AHA+ADT+HAB,
data=CGData,family="binomial")
## Creating the summary of all models
NP12_sum=c(NP12$call,NP12$family$family,NP12$family$link,
NP12$aic,NP12$method,NP12$converged,(rsq=rsq(NP12,adj=FALSE,type=c('lr'))),
(rsq=rsq(NP12,adj=FALSE,type=c('n'))),NP12$df.residual,NP12$df.null,
NP12$null.deviance,NP12$deviance,"NP12")
NP11_sum=c(NP11$call,NP11$family$family,NP11$family$link,
NP11$aic,NP11$method,NP11$converged,(rsq=rsq(NP11,adj=FALSE,type=c('lr'))),
(rsq=rsq(NP11,adj=FALSE,type=c('n'))),NP11$df.residual,NP11$df.null,
NP11$null.deviance,NP11$deviance,"NP11")
NP10_sum=c(NP10$call,NP10$family$family,NP10$family$link,
NP10$aic,NP10$method,NP10$converged,(rsq=rsq(NP10,adj=FALSE,type=c('lr'))),
(rsq=rsq(NP10,adj=FALSE,type=c('n'))),NP10$df.residual,NP10$df.null,
NP10$null.deviance,NP10$deviance,"NP10")
NP09_sum=c(NP09$call,NP09$family$family,NP09$family$link,
NP09$aic,NP09$method,NP09$converged,(rsq=rsq(NP09,adj=FALSE,type=c('lr'))),
(rsq=rsq(NP09,adj=FALSE,type=c('n'))),NP09$df.residual,NP09$df.null,
NP09$null.deviance,NP09$deviance,"NP09")
NP08_sum=c(NP08$call,NP08$family$family,NP08$family$link,
NP08$aic,NP08$method,NP08$converged,(rsq=rsq(NP08,adj=FALSE,type=c('lr'))),
(rsq=rsq(NP08,adj=FALSE,type=c('n'))),NP08$df.residual,NP08$df.null,
NP08$null.deviance,NP08$deviance,"NP08")
NP07_sum=c(NP07$call,NP07$family$family,NP07$family$link,
NP07$aic,NP07$method,NP07$converged,(rsq=rsq(NP07,adj=FALSE,type=c('lr'))),
(rsq=rsq(NP07,adj=FALSE,type=c('n'))),NP07$df.residual,NP07$df.null,
NP07$null.deviance,NP07$deviance,"NP07")
NP06_sum=c(NP06$call,NP06$family$family,NP06$family$link,
NP06$aic,NP06$method,NP06$converged,(rsq=rsq(NP06,adj=FALSE,type=c('lr'))),
(rsq=rsq(NP06,adj=FALSE,type=c('n'))),NP06$df.residual,NP06$df.null,
NP06$null.deviance,NP06$deviance,"NP06")
NP05_sum=c(NP05$call,NP05$family$family,NP05$family$link,
NP05$aic,NP05$method,NP05$converged,(rsq=rsq(NP05,adj=FALSE,type=c('lr'))),
(rsq=rsq(NP05,adj=FALSE,type=c('n'))),NP05$df.residual,NP05$df.null,
NP05$null.deviance,NP05$deviance,"NP05")
NP04_sum=c(NP04$call,NP04$family$family,NP04$family$link,
NP04$aic,NP04$method,NP04$converged,(rsq=rsq(NP04,adj=FALSE,type=c('lr'))),
(rsq=rsq(NP04,adj=FALSE,type=c('n'))),NP04$df.residual,NP04$df.null,
NP04$null.deviance,NP04$deviance,"NP04")
NP03_sum=c(NP03$call,NP03$family$family,NP03$family$link,
NP03$aic,NP03$method,NP03$converged,(rsq=rsq(NP03,adj=FALSE,type=c('lr'))),
(rsq=rsq(NP03,adj=FALSE,type=c('n'))),NP03$df.residual,NP03$df.null,
NP03$null.deviance,NP03$deviance,"NP03")
NP02_sum=c(NP02$call,NP02$family$family,NP02$family$link,
NP02$aic,NP02$method,NP02$converged,(rsq=rsq(NP02,adj=FALSE,type=c('lr'))),
(rsq=rsq(NP02,adj=FALSE,type=c('n'))),NP02$df.residual,NP02$df.null,
NP02$null.deviance,NP02$deviance,"NP02")
NP01_sum=c(NP01$call,NP01$family$family,NP01$family$link,
NP01$aic,NP01$method,NP01$converged,(rsq=rsq(NP01,adj=FALSE,type=c('lr'))),
(rsq=rsq(NP01,adj=FALSE,type=c('n'))),NP01$df.residual,NP01$df.null,
NP01$null.deviance,NP01$deviance,"NP01")
## Creating a dataframe for the summary table of all models
GLM_Summary_Dataframe=data.frame(rbind(NP12_sum,NP11_sum,NP10_sum,NP09_sum,NP08_sum,NP07_sum,
NP06_sum, NP05_sum, NP04_sum, NP03_sum,
NP02_sum,NP01_sum))
colnames(GLM_Summary_Dataframe)=c("Call","family","link","AIC","method","converged",
"RsqrLR","RsqrNag_Cragg.Uhler","df.residual","df.null",
"null.deviance","deviance","Model")3.3.1 Printing the GLM results: Model expressions
options(scipen = 1, digits = 2) #set to two decimal
knitr::kable(GLM_Summary_Dataframe[,c(1,13)], format = "markdown")| Call | Model | |
|---|---|---|
| NP12_sum | glm(formula = NP ~ CTA + AHA + CHC + CHS + CAC + CAS + LA + C2x2 + , C4x4 + ADT, family = “binomial”, data = CGData) | NP12 |
| NP11_sum | glm(formula = NP ~ CTA + AHA + CHC + CHS + C2x2 + C4x4, family = “binomial”, , data = CGData) | NP11 |
| NP10_sum | glm(formula = NP ~ CTA + AHA + CHC + CHS + LA, family = “binomial”, , data = CGData) | NP10 |
| NP09_sum | glm(formula = NP ~ CTA + AHA + CHC + CHS, family = “binomial”, , data = CGData) | NP09 |
| NP08_sum | glm(formula = NP ~ CTA + AHA + CHC + CHS + CAC + CAS + ADT + , HAB + HAB:AHA, family = “binomial”, data = CGData) | NP08 |
| NP07_sum | glm(formula = NP ~ CTA + CHC + CHS + CAC + LA + CAS + C2x2 + , C4x4 + AHA + ADT + HAB, family = “binomial”, data = CGData) | NP07 |
| NP06_sum | glm(formula = NP ~ CTA + AHA + CHC + CHS + CAC + CAS + ADT + , HAB + HAB:CHS, family = “binomial”, data = CGData) | NP06 |
| NP05_sum | glm(formula = NP ~ CTA + AHA + CHC + CHS + CAC + CAS + ADT + , HAB + HAB:CTA, family = “binomial”, data = CGData) | NP05 |
| NP04_sum | glm(formula = NP ~ CTA + AHA + CHC + CHS + CAC + CAS + ADT + , HAB + HAB:CHC + HAB:AHA, family = “binomial”, data = CGData) | NP04 |
| NP03_sum | glm(formula = NP ~ CTA + AHA + CHC + CHS + CAC + CAS + ADT + , HAB + HAB:CHC, family = “binomial”, data = CGData) | NP03 |
| NP02_sum | glm(formula = NP ~ CTA + AHA + CHC + CHS + CAC + CAS + ADT + , HAB + HAB:CHC + HAB:CHS, family = “binomial”, data = CGData) | NP02 |
| NP01_sum | glm(formula = NP ~ CTA + AHA + CHC + CHS + CAC + CAS + ADT + , HAB + HAB:CHC + HAB:CTA, family = “binomial”, data = CGData) | NP01 |
3.3.2 Printing the GLM results: Model summaries
options(scipen = 1, digits = 2) #set to two decimal
knitr::kable(GLM_Summary_Dataframe[,-c(1,13)], format = "markdown")| family | link | AIC | method | converged | RsqrLR | RsqrNag_Cragg.Uhler | df.residual | df.null | null.deviance | deviance | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| NP12_sum | binomial | logit | 278 | glm.fit | TRUE | 0.33 | 0.58 | 564 | 574 | 485 | 256 |
| NP11_sum | binomial | logit | 277 | glm.fit | TRUE | 0.32 | 0.56 | 568 | 574 | 485 | 263 |
| NP10_sum | binomial | logit | 274 | glm.fit | TRUE | 0.32 | 0.56 | 569 | 574 | 485 | 262 |
| NP09_sum | binomial | logit | 274 | glm.fit | TRUE | 0.32 | 0.56 | 570 | 574 | 485 | 264 |
| NP08_sum | binomial | logit | 250 | glm.fit | TRUE | 0.37 | 0.64 | 561 | 574 | 485 | 222 |
| NP07_sum | binomial | logit | 248 | glm.fit | TRUE | 0.37 | 0.65 | 561 | 574 | 485 | 220 |
| NP06_sum | binomial | logit | 231 | glm.fit | TRUE | 0.39 | 0.68 | 561 | 574 | 485 | 203 |
| NP05_sum | binomial | logit | 226 | glm.fit | TRUE | 0.39 | 0.69 | 561 | 574 | 485 | 198 |
| NP04_sum | binomial | logit | 179 | glm.fit | TRUE | 0.45 | 0.78 | 558 | 574 | 485 | 145 |
| NP03_sum | binomial | logit | 177 | glm.fit | TRUE | 0.44 | 0.78 | 561 | 574 | 485 | 149 |
| NP02_sum | binomial | logit | 177 | glm.fit | TRUE | 0.45 | 0.79 | 558 | 574 | 485 | 143 |
| NP01_sum | binomial | logit | 164 | glm.fit | TRUE | 0.46 | 0.81 | 558 | 574 | 485 | 130 |
3.3.3 Choosing the best model
We considered the best model as the model with the smallest AIC
options(scipen = 1, digits = 2) #set to two decimal
best_model = GLM_Summary_Dataframe[GLM_Summary_Dataframe$AIC==min(unlist(GLM_Summary_Dataframe$AIC)),c(1,13)]
knitr::kable(best_model, format = "markdown")| Call | Model | |
|---|---|---|
| NP01_sum | glm(formula = NP ~ CTA + AHA + CHC + CHS + CAC + CAS + ADT + , HAB + HAB:CHC + HAB:CTA, family = “binomial”, data = CGData) | NP01 |
The best model is NP01.
3.4 Determining the confidence intervals of the coefficients, and Odd ratios
options(scipen = 1, digits = 2) #set to two decimal
## Creating the different data frames for the summary table for the log and odd ratio
NP01_confit_log=round(as.data.frame((cbind(LogOdds = coef(NP01), confint(NP01)))),
digits=3)
NP01_confit_Odds=round(as.data.frame(exp(cbind(OddRatio = coef(NP01),
confint(NP01)))),digits=3)
## Combining the data frames into one
NP01_confit_All1=cbind(NP01_confit_Odds,NP01_confit_log)
NP01_confit_All1$Names=rownames(NP01_confit_All1)
colnames(NP01_confit_All1)=c("OddRatio","OR2.5","OR97.5","LogOdds","LogOR2.5","LogOR97.5","Names" )
rownames(NP01_confit_All1)=NULL
## Changing the order of columns
NP01_confit_All=NP01_confit_All1[,c("Names","OddRatio","OR2.5","OR97.5","LogOdds","LogOR2.5","LogOR97.5")]
## Importing data into an excel file.
List_GLM_Results_confit <- list("NP01_confit_All"= NP01_confit_All)
write.xlsx(List_GLM_Results_confit, file = "NP01_confit_All.xlsx")
## Printing the coefficient table
knitr::kable(NP01_confit_All, format = "markdown")| Names | OddRatio | OR2.5 | OR97.5 | LogOdds | LogOR2.5 | LogOR97.5 |
|---|---|---|---|---|---|---|
| (Intercept) | 0.22 | 0.00 | 44.12 | -1.53 | -6.96 | 3.79 |
| CTA | 0.82 | 0.68 | 0.97 | -0.20 | -0.39 | -0.03 |
| AHA | 0.76 | 0.64 | 0.90 | -0.27 | -0.45 | -0.10 |
| CHC | 11.86 | 5.44 | 33.58 | 2.47 | 1.69 | 3.51 |
| CHS | 7.24 | 3.77 | 15.75 | 1.98 | 1.33 | 2.76 |
| CAC | 0.82 | 0.60 | 1.08 | -0.19 | -0.51 | 0.08 |
| CAS | 0.96 | 0.52 | 1.70 | -0.04 | -0.65 | 0.53 |
| ADT | 1.16 | 0.99 | 1.37 | 0.15 | -0.01 | 0.32 |
| HABRIP | 0.44 | 0.00 | 86.84 | -0.81 | -5.92 | 4.46 |
| HABSW | 0.04 | 0.00 | 3.51 | -3.17 | -7.86 | 1.26 |
| HABYSF | 0.37 | 0.00 | 46.37 | -0.99 | -5.85 | 3.84 |
| CHC:HABRIP | 0.01 | 0.00 | 0.09 | -4.31 | -7.26 | -2.38 |
| CHC:HABSW | 0.66 | 0.06 | 7.84 | -0.41 | -2.77 | 2.06 |
| CHC:HABYSF | 0.10 | 0.03 | 0.24 | -2.31 | -3.46 | -1.42 |
| CTA:HABRIP | 1.52 | 1.00 | 2.52 | 0.42 | 0.00 | 0.92 |
| CTA:HABSW | 1.65 | 1.29 | 2.17 | 0.50 | 0.26 | 0.78 |
| CTA:HABYSF | 1.20 | 0.93 | 1.55 | 0.18 | -0.07 | 0.44 |
3.5 Overall effects of habitat and the interactions
options(scipen = 1, digits = 2) #set to two decimal
## Note: 40 represents the magnitude of the data in CTA
Chimp_GLM_Data_CTA <- with(Chimp_GLM_Data, data.frame(CTA = rep(seq(from = 0, to = 40, length.out = 100),4),
AHA = mean(AHA), CHC = mean(CHC),
CHS = mean(CHS), CAC = mean(CAC),
CAS = mean(CAS), ADT = mean(ADT),
HAB = factor(rep(c("OSF","RIP","SW","YSF"), each = 100))))
Chimp_GLM_Data_AHA <- with(Chimp_GLM_Data, data.frame(AHA = rep(seq(from = 0, to = 40, length.out = 100),4),
CTA = mean(CTA), CHC = mean(CHC),
CHS = mean(CHS), CAC = mean(CAC),
CAS = mean(CAS), ADT = mean(ADT),
HAB = factor(rep(c("OSF","RIP","SW","YSF"), each = 100))))
Chimp_GLM_Data_CHC <- with(Chimp_GLM_Data, data.frame(CHC = rep(seq(from = 0, to = 5, length.out = 100),4),
AHA = mean(AHA), CTA = mean(CTA),
CHS = mean(CHS), CAC = mean(CAC),
CAS = mean(CAS), ADT = mean(ADT),
HAB = factor(rep(c("OSF","RIP","SW","YSF"), each = 100))))
Chimp_GLM_Data_CHS <- with(Chimp_GLM_Data, data.frame(CHS = rep(seq(from = 0, to = 8, length.out = 100),4),
AHA = mean(AHA), CHC = mean(CHC),
CTA = mean(CTA), CAC = mean(CAC),
CAS = mean(CAS), ADT = mean(ADT),
HAB = factor(rep(c("OSF","RIP","SW","YSF"), each = 100))))
knitr::kable(Chimp_GLM_Data_CTA[c(1:10),], format = "markdown")| CTA | AHA | CHC | CHS | CAC | CAS | ADT | HAB |
|---|---|---|---|---|---|---|---|
| 0.00 | 21 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF |
| 0.40 | 21 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF |
| 0.81 | 21 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF |
| 1.21 | 21 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF |
| 1.62 | 21 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF |
| 2.02 | 21 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF |
| 2.42 | 21 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF |
| 2.83 | 21 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF |
| 3.23 | 21 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF |
| 3.64 | 21 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF |
3.6 Calculate the predicted values
Chimp_GLM_Data_CTA_pred_1 <- cbind(Chimp_GLM_Data_CTA,
predict(NP01, newdata = Chimp_GLM_Data_CTA,
type = "link",se = TRUE))
Chimp_GLM_Data_AHA_pred_1 <- cbind(Chimp_GLM_Data_AHA,
predict(NP01, newdata = Chimp_GLM_Data_AHA,
type = "link",se = TRUE))
Chimp_GLM_Data_CHC_pred_1 <- cbind(Chimp_GLM_Data_CHC,
predict(NP01, newdata = Chimp_GLM_Data_CHC,
type = "link",se = TRUE))
Chimp_GLM_Data_CHS_pred_1 <- cbind(Chimp_GLM_Data_CHS,
predict(NP01, newdata = Chimp_GLM_Data_CHS,
type = "link",se = TRUE))3.7 Predicted values and confidence limits into probabilities
Chimp_GLM_Data_CTA_pred <- within(Chimp_GLM_Data_CTA_pred_1, {
PredictedProb <- plogis(fit)
LL <- plogis(fit - (1.96 * se.fit))
UL <- plogis(fit + (1.96 * se.fit))
})
Chimp_GLM_Data_AHA_pred <- within(Chimp_GLM_Data_AHA_pred_1, {
PredictedProb <- plogis(fit)
LL <- plogis(fit - (1.96 * se.fit))
UL <- plogis(fit + (1.96 * se.fit))
})
Chimp_GLM_Data_CHC_pred <- within(Chimp_GLM_Data_CHC_pred_1, {
PredictedProb <- plogis(fit)
LL <- plogis(fit - (1.96 * se.fit))
UL <- plogis(fit + (1.96 * se.fit))
})
Chimp_GLM_Data_CHS_pred <- within(Chimp_GLM_Data_CHS_pred_1, {
PredictedProb <- plogis(fit)
LL <- plogis(fit - (1.96 * se.fit))
UL <- plogis(fit + (1.96 * se.fit))
})
knitr::kable(Chimp_GLM_Data_AHA_pred[c(1:10),], format = "markdown")| AHA | CTA | CHC | CHS | CAC | CAS | ADT | HAB | fit | se.fit | residual.scale | UL | LL | PredictedProb |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.00 | 18 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF | 3.3 | 1.9 | 1 | 1 | 0.39 | 0.96 |
| 0.40 | 18 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF | 3.2 | 1.9 | 1 | 1 | 0.38 | 0.96 |
| 0.81 | 18 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF | 3.1 | 1.8 | 1 | 1 | 0.37 | 0.96 |
| 1.21 | 18 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF | 3.0 | 1.8 | 1 | 1 | 0.36 | 0.95 |
| 1.62 | 18 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF | 2.9 | 1.8 | 1 | 1 | 0.35 | 0.95 |
| 2.02 | 18 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF | 2.8 | 1.7 | 1 | 1 | 0.34 | 0.94 |
| 2.42 | 18 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF | 2.6 | 1.7 | 1 | 1 | 0.33 | 0.93 |
| 2.83 | 18 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF | 2.5 | 1.7 | 1 | 1 | 0.33 | 0.93 |
| 3.23 | 18 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF | 2.4 | 1.6 | 1 | 1 | 0.32 | 0.92 |
| 3.64 | 18 | 1.6 | 0.56 | 3.6 | 0.95 | 27 | OSF | 2.3 | 1.6 | 1 | 1 | 0.31 | 0.91 |
3.8 Plotting the GLM results
The codes below will create the plot of the predicted probabilities of chimpanzee
nest presence in response to the variables (CTA, AHA, CHC and CHS) and in
relation to the different habitat types (OSF, RIP, SW and YSF). The 4
plots produced will be arranged into one, and the final plot will be saved into
different file formats (tif, pdf and eps).
Chimp_GLM_Data_CTA_pred$probability="Probability"
## Plot the results without interaction
Chimp_GLM_CTA_gg_wint=ggplot(Chimp_GLM_Data_CTA_pred, aes(x = CTA, y = PredictedProb)) +
#geom_ribbon(aes(ymin = LL,ymax = UL), alpha = 0.2) +
geom_line(aes(colour=HAB, linetype=HAB),size = 1)+theme_bw()+
labs(y = "Predicted probability", x="Count of all tree stems")+
labs(fill = "Habitat",color="Habitat")+theme(legend.position = "none")+ # Expand y range
scale_colour_manual(values = c("#1B9E77", "#D95F02", "#E6AB02", "#666666", "#7570B3","#F8766D",
"#A6761D","#a6d854",
"#00BF7D","#E7298A"))
Chimp_GLM_AHA_gg_wint=ggplot(Chimp_GLM_Data_AHA_pred, aes(x = AHA, y = PredictedProb)) +
#geom_ribbon(aes(ymin = LL,ymax = UL), alpha = 0.2) +
geom_line(aes(colour = HAB, linetype=HAB),size = 1)+theme_bw()+
labs(y = "Predicted probability", x="Average height of all trees")+
labs(fill = "Habitat",color="Habitat")+theme(legend.position = "none")+ # Expand y range
scale_colour_manual(values = c("#1B9E77", "#D95F02", "#E6AB02", "#666666", "#7570B3","#F8766D",
"#A6761D","#a6d854",
"#00BF7D","#E7298A"))
Chimp_GLM_CHC_gg_wint=ggplot(Chimp_GLM_Data_CHC_pred, aes(x = CHC, y = PredictedProb)) +
#geom_ribbon(aes(ymin = LL,ymax = UL), alpha = 0.2) +
geom_line(aes(colour = HAB, linetype=HAB),size = 1)+theme_bw()+
labs(y = "Predicted probability", x="Number of neutral trees per habitat")+
labs(fill = "Habitat",color="Habitat")+theme(legend.position = "none")+ # Expand y range
scale_colour_manual(values = c("#1B9E77", "#D95F02", "#E6AB02", "#666666", "#7570B3","#F8766D",
"#A6761D","#a6d854",
"#00BF7D","#E7298A"))
Chimp_GLM_CHS_gg_wint=ggplot(Chimp_GLM_Data_CHS_pred, aes(x = CHS, y = PredictedProb)) +
#geom_ribbon(aes(ymin = LL,ymax = UL), alpha = 0.2) +
labs(y = "Predicted probability", x="Number of selected trees per habitat")+
geom_line(aes(colour = HAB, linetype=HAB),size = 1)+theme_bw()+
labs(fill = "Habitat",color="Habitat",linetype="Habitat")+theme(legend.position = c(0.7, 0.38))+
theme(#legend.key.height= unit(2, 'cm'),
legend.key.width= unit(1.5, 'cm'))+ # Expand y range
scale_colour_manual(values = c("#1B9E77", "#D95F02", "#E6AB02", "#666666", "#7570B3","#F8766D",
"#A6761D","#a6d854",
"#00BF7D","#E7298A"))
###################################################
## Arranging the 4 plots into 1
Chimp_GLM_gg_arrange_wint <- ggarrange(Chimp_GLM_CHS_gg_wint,
Chimp_GLM_CHC_gg_wint,
Chimp_GLM_AHA_gg_wint,
Chimp_GLM_CTA_gg_wint,
heights = c(1,1,1,1),labels = c("(a)", "(b)", "(c)", "(d)"),
ncol = 2, nrow = 2,
font.label = list(size = 9, color = "black",
face = "bold",family = NULL))
## Saving the plot into different file formats
setwd("C:/Users/tedon/OneDrive - UGent/PHD Luc/Articles en cours/Ecological factors of nesting site choice/Results")
## Plots for clusters
tiff(file = "Chimp_GLM_gg_arrange_wint3.tif", width = 17.5, height = 13, units = "cm", pointsize = 10,
bg = "white", res = 300)
Chimp_GLM_gg_arrange_wint
dev.off()## png
## 2HAb_dir1=setwd("C:/Users/tedon/OneDrive - UGent/PHD Luc/Articles en cours/Ecological factors of nesting site choice/Results")
ggsave("Chimp_GLM_gg_arrange_wint3.pdf",Chimp_GLM_gg_arrange_wint, path = HAb_dir1,
scale = 1, width = 17.5, height = 13, units = "cm",
dpi = 600)
ggsave("Chimp_GLM_gg_arrange_wint3.eps",Chimp_GLM_gg_arrange_wint, path = HAb_dir1,
scale = 1, width = 17.5, height = 13, units = "cm",
dpi = 600)
Figure 1: Variation of the probability of chimpanzee nesting site selection in relation to the scores of the variables with significant effects in the logistic regression
OSF old secondary forest, YSF young secondary forest, SW swamp, RIP riparian forest.
4 Tree DBH preference for chimpanzee nesting
4.1 Preparing the data
options(scipen = 1, digits = 2) #set to two decimal
library(adehabitatHS)
Chimp_avec_Nid="C:/Users/tedon/OneDrive - UGent/PHD Luc/Articles en cours/Ecological factors of nesting site choice/Data/Botany/Analyses_Botany-Parcelles avec nids June 2021 (24-07-2021).xlsx"
S1_AT_NP_AT_SP=readWorkbook(xlsxFile = Chimp_avec_Nid, sheet = "DBH_Used_SelectionType")
S2_NTS_NP_NTS_SP=readWorkbook(xlsxFile = Chimp_avec_Nid, sheet = "DBH_Used_SelectionType_NestSpec")
S3_TU_NB_AT_NP_Global<- readWorkbook(xlsxFile = Chimp_avec_Nid, sheet = "DBH_Used_SelectionGlobal")
S3_TU_NB_AT_NP_Tree=readWorkbook(xlsxFile = Chimp_avec_Nid, sheet = "DBH_Used_SelectionTree")
S4_UT_TA_NTS_NP=readWorkbook(xlsxFile = Chimp_avec_Nid, sheet = "DBH_Used_SelectionType_NestTree")
## Arranging the data
## S1_AT_NP_AT_SP: All trees in nesting plots vs all trees in systematic plots
S1_AT_NP_AT_SP_Global=S1_AT_NP_AT_SP[S1_AT_NP_AT_SP$Habitat=="Global",];rownames(S1_AT_NP_AT_SP_Global)=S1_AT_NP_AT_SP_Global$Classe_DBH70;S1_AT_NP_AT_SP_Global_Use=S1_AT_NP_AT_SP_Global$PAN;S1_AT_NP_AT_SP_Global_avail=S1_AT_NP_AT_SP_Global$PSN;names(S1_AT_NP_AT_SP_Global_Use)=S1_AT_NP_AT_SP_Global$Classe_DBH70;names(S1_AT_NP_AT_SP_Global_avail)=S1_AT_NP_AT_SP_Global$Classe_DBH70
## S3 = S3_TU_NB_AT_NP: Trees used for nest building vs all trees available in nesting plots
S3_TU_NB_AT_NP_Glob=S3_TU_NB_AT_NP_Global$Used;S3_TU_NB_AT_NP_Glob_avail=S3_TU_NB_AT_NP_Global$Avail;names(S3_TU_NB_AT_NP_Glob)=S3_TU_NB_AT_NP_Global$DBH;names(S3_TU_NB_AT_NP_Glob_avail)=S3_TU_NB_AT_NP_Global$DBH
## S2 = S2_NTS_NP_NTS_SP: Nesting tree species in nesting plots vs nesting tree species in systematic plots
S2_NTS_NP_NTS_SP_Global=S2_NTS_NP_NTS_SP[S2_NTS_NP_NTS_SP$Habitat=="Global",];rownames(S2_NTS_NP_NTS_SP_Global)=S2_NTS_NP_NTS_SP_Global$Classe_DBH70;S2_NTS_NP_NTS_SP_Global_Use=S2_NTS_NP_NTS_SP_Global$PAN;S2_NTS_NP_NTS_SP_Global_avail=S2_NTS_NP_NTS_SP_Global$PSN;names(S2_NTS_NP_NTS_SP_Global_Use)=S2_NTS_NP_NTS_SP_Global$Classe_DBH70;names(S2_NTS_NP_NTS_SP_Global_avail)=S2_NTS_NP_NTS_SP_Global$Classe_DBH70
## S4 = S4_UT_TA_NTS_NP: Used trees vs trees available of nesting tree species in nesting plots
S4_UT_TA_NTS_NP_Global=S4_UT_TA_NTS_NP[S4_UT_TA_NTS_NP$Habitat=="Global",];rownames(S4_UT_TA_NTS_NP_Global)=S4_UT_TA_NTS_NP_Global$Classe_DBH70;S4_UT_TA_NTS_NP_Global_Use=S4_UT_TA_NTS_NP_Global$Used;S4_UT_TA_NTS_NP_Global_avail=S4_UT_TA_NTS_NP_Global$Avail;names(S4_UT_TA_NTS_NP_Global_Use)=S4_UT_TA_NTS_NP_Global$Classe_DBH70;names(S4_UT_TA_NTS_NP_Global_avail)=S4_UT_TA_NTS_NP_Global$Classe_DBH70
## Showing an example
knitr::kable(S4_UT_TA_NTS_NP_Global, format = "markdown")| Habitat | Classe_DBH70 | No | Avail | Used | |
|---|---|---|---|---|---|
| [10-20[ | Global | [10-20[ | 1205 | 297 | 27 |
| [20-30[ | Global | [20-30[ | 725 | 154 | 26 |
| [30-40[ | Global | [30-40[ | 437 | 83 | 20 |
| [40-50[ | Global | [40-50[ | 176 | 47 | 19 |
| [50-60[ | Global | [50-60[ | 53 | 12 | 4 |
| [60-70[ | Global | [60-70[ | 69 | 15 | 4 |
| [70-more[ | Global | [70-more[ | 10 | 7 | 1 |
4.2 Calculating the DBH preference with the function widesI2 package adehabitatHS
options(scipen = 1, digits = 2) #set to two decimal
## S1_AT_NP_AT_SP: All trees in nesting plots vs all trees in systematic plots
wi_S1_AT_NP_AT_SP_Glob <- widesI(S1_AT_NP_AT_SP_Global_Use,S1_AT_NP_AT_SP_Global_avail,alpha = 0.05); wi_S1_AT_NP_AT_SP_Glob_df=cbind(data.frame(wi_S1_AT_NP_AT_SP_Glob$wi), data.frame(wi_S1_AT_NP_AT_SP_Glob$se.wi*2), data.frame(wi_S1_AT_NP_AT_SP_Glob$chisquwi)); colnames(wi_S1_AT_NP_AT_SP_Glob_df)=c("wi", "se.wix2", "testwi","p"); wi_S1_AT_NP_AT_SP_Glob_df$DBH=rownames(wi_S1_AT_NP_AT_SP_Glob_df); wi_S1_AT_NP_AT_SP_Glob_df$Animal="Chimpanzee"
wi_S1_AT_NP_AT_SP_Glob_df$upper=wi_S1_AT_NP_AT_SP_Glob_df$wi+wi_S1_AT_NP_AT_SP_Glob_df$se.wix2;wi_S1_AT_NP_AT_SP_Glob_df$lower=wi_S1_AT_NP_AT_SP_Glob_df$wi-wi_S1_AT_NP_AT_SP_Glob_df$se.wix2
wi_S1_AT_NP_AT_SP_Glob_df$Scenario="S1_AT_NP_AT_SP";wi_S1_AT_NP_AT_SP_Glob_df$Habitat="Global";wi_S1_AT_NP_AT_SP_Glob_df$Selection=ifelse(wi_S1_AT_NP_AT_SP_Glob_df$upper>1 & wi_S1_AT_NP_AT_SP_Glob_df$lower>1,"Selected",ifelse(!wi_S1_AT_NP_AT_SP_Glob_df$upper>1 & !wi_S1_AT_NP_AT_SP_Glob_df$lower>1,"Avoided","Common"))
## S3 = S3_TU_NB_AT_NP: Trees used for nest building vs all trees available in nesting plots
wi_S3_TU_NB_AT_NP_Glob <- widesI(S3_TU_NB_AT_NP_Glob,S3_TU_NB_AT_NP_Glob_avail,alpha = 0.05); wi_S3_TU_NB_AT_NP_Glob_df=cbind(data.frame(wi_S3_TU_NB_AT_NP_Glob$wi), data.frame(wi_S3_TU_NB_AT_NP_Glob$se.wi*2), data.frame(wi_S3_TU_NB_AT_NP_Glob$chisquwi)); colnames(wi_S3_TU_NB_AT_NP_Glob_df)=c("wi", "se.wix2", "testwi","p"); wi_S3_TU_NB_AT_NP_Glob_df$DBH=rownames(wi_S3_TU_NB_AT_NP_Glob_df); wi_S3_TU_NB_AT_NP_Glob_df$Animal="Chimpanzee"
wi_S3_TU_NB_AT_NP_Glob_df$upper=wi_S3_TU_NB_AT_NP_Glob_df$wi+wi_S3_TU_NB_AT_NP_Glob_df$se.wix2;wi_S3_TU_NB_AT_NP_Glob_df$lower=wi_S3_TU_NB_AT_NP_Glob_df$wi-wi_S3_TU_NB_AT_NP_Glob_df$se.wix2
wi_S3_TU_NB_AT_NP_Glob_df$Scenario="S3_TU_NB_AT_NP";wi_S3_TU_NB_AT_NP_Glob_df$Habitat="Global";wi_S3_TU_NB_AT_NP_Glob_df$Selection=ifelse(wi_S3_TU_NB_AT_NP_Glob_df$upper>1 & wi_S3_TU_NB_AT_NP_Glob_df$lower>1,"Selected",ifelse(!wi_S3_TU_NB_AT_NP_Glob_df$upper>1 & !wi_S3_TU_NB_AT_NP_Glob_df$lower>1,"Avoided","Common"))
## S2 = S2_NTS_NP_NTS_SP: Nesting tree species in nesting plots vs nesting tree species in systematic plots
wi_S2_NTS_NP_NTS_SPGlob <- widesI(S2_NTS_NP_NTS_SP_Global_Use,S2_NTS_NP_NTS_SP_Global_avail,alpha = 0.05); wi_S2_NTS_NP_NTS_SPGlob_df=cbind(data.frame(wi_S2_NTS_NP_NTS_SPGlob$wi), data.frame(wi_S2_NTS_NP_NTS_SPGlob$se.wi*2), data.frame(wi_S2_NTS_NP_NTS_SPGlob$chisquwi)); colnames(wi_S2_NTS_NP_NTS_SPGlob_df)=c("wi", "se.wix2", "testwi","p"); wi_S2_NTS_NP_NTS_SPGlob_df$DBH=rownames(wi_S2_NTS_NP_NTS_SPGlob_df); wi_S2_NTS_NP_NTS_SPGlob_df$Animal="Chimpanzee"
wi_S2_NTS_NP_NTS_SPGlob_df$upper=wi_S2_NTS_NP_NTS_SPGlob_df$wi+wi_S2_NTS_NP_NTS_SPGlob_df$se.wix2;wi_S2_NTS_NP_NTS_SPGlob_df$lower=wi_S2_NTS_NP_NTS_SPGlob_df$wi-wi_S2_NTS_NP_NTS_SPGlob_df$se.wix2
wi_S2_NTS_NP_NTS_SPGlob_df$Scenario="S2_NTS_NP_NTS_SP";wi_S2_NTS_NP_NTS_SPGlob_df$Habitat="Global";wi_S2_NTS_NP_NTS_SPGlob_df$Selection=ifelse(wi_S2_NTS_NP_NTS_SPGlob_df$upper>1 & wi_S2_NTS_NP_NTS_SPGlob_df$lower>1,"Selected",ifelse(!wi_S2_NTS_NP_NTS_SPGlob_df$upper>1 & !wi_S2_NTS_NP_NTS_SPGlob_df$lower>1,"Avoided","Common"))
## S4 = S4_UT_TA_NTS_NP: Used trees vs trees available of nesting tree species in nesting plots
wi_S4_UT_TA_NTS_NPGlob <- widesI(S4_UT_TA_NTS_NP_Global_Use,S4_UT_TA_NTS_NP_Global_avail,alpha = 0.05); wi_S4_UT_TA_NTS_NPGlob_df=cbind(data.frame(wi_S4_UT_TA_NTS_NPGlob$wi), data.frame(wi_S4_UT_TA_NTS_NPGlob$se.wi*2), data.frame(wi_S4_UT_TA_NTS_NPGlob$chisquwi)); colnames(wi_S4_UT_TA_NTS_NPGlob_df)=c("wi", "se.wix2", "testwi","p"); wi_S4_UT_TA_NTS_NPGlob_df$DBH=rownames(wi_S4_UT_TA_NTS_NPGlob_df); wi_S4_UT_TA_NTS_NPGlob_df$Animal="Chimpanzee"
wi_S4_UT_TA_NTS_NPGlob_df$upper=wi_S4_UT_TA_NTS_NPGlob_df$wi+wi_S4_UT_TA_NTS_NPGlob_df$se.wix2;wi_S4_UT_TA_NTS_NPGlob_df$lower=wi_S4_UT_TA_NTS_NPGlob_df$wi-wi_S4_UT_TA_NTS_NPGlob_df$se.wix2
wi_S4_UT_TA_NTS_NPGlob_df$Scenario="S4_UT_TA_NTS_NP";wi_S4_UT_TA_NTS_NPGlob_df$Habitat="Global";wi_S4_UT_TA_NTS_NPGlob_df$Selection=ifelse(wi_S4_UT_TA_NTS_NPGlob_df$upper>1 & wi_S4_UT_TA_NTS_NPGlob_df$lower>1,"Selected",ifelse(!wi_S4_UT_TA_NTS_NPGlob_df$upper>1 & !wi_S4_UT_TA_NTS_NPGlob_df$lower>1,"Avoided","Common"))
## Showing an example
knitr::kable(wi_S4_UT_TA_NTS_NPGlob_df, format = "markdown")| wi | se.wix2 | testwi | p | DBH | Animal | upper | lower | Scenario | Habitat | Selection | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| [10-20[ | 0.55 | 0.18 | 23.97 | 0.00 | [10-20[ | Chimpanzee | 0.74 | 0.37 | S4_UT_TA_NTS_NP | Global | Avoided |
| [20-30[ | 1.03 | 0.35 | 0.03 | 0.87 | [20-30[ | Chimpanzee | 1.38 | 0.68 | S4_UT_TA_NTS_NP | Global | Common |
| [30-40[ | 1.47 | 0.59 | 2.53 | 0.11 | [30-40[ | Chimpanzee | 2.05 | 0.88 | S4_UT_TA_NTS_NP | Global | Common |
| [40-50[ | 2.46 | 1.02 | 8.25 | 0.00 | [40-50[ | Chimpanzee | 3.48 | 1.44 | S4_UT_TA_NTS_NP | Global | Selected |
| [50-60[ | 2.03 | 1.99 | 1.07 | 0.30 | [50-60[ | Chimpanzee | 4.02 | 0.04 | S4_UT_TA_NTS_NP | Global | Common |
| [60-70[ | 1.62 | 1.59 | 0.61 | 0.43 | [60-70[ | Chimpanzee | 3.22 | 0.03 | S4_UT_TA_NTS_NP | Global | Common |
| [70-more[ | 0.87 | 1.73 | 0.02 | 0.88 | [70-more[ | Chimpanzee | 2.60 | -0.86 | S4_UT_TA_NTS_NP | Global | Common |
4.3 Ploting the result
## Preparing the data
wi_DBH_ALL_Global=rbind(wi_S1_AT_NP_AT_SP_Glob_df, wi_S4_UT_TA_NTS_NPGlob_df, wi_S2_NTS_NP_NTS_SPGlob_df,
wi_S3_TU_NB_AT_NP_Glob_df)
wi_DBH_ALL_Global$Scenario2=wi_DBH_ALL_Global$Scenario
wi_DBH_ALL_Global$Selection2=wi_DBH_ALL_Global$Selection
t2.rect1_DBH <- data.frame (xmin=3.5, xmax=4.5, ymin=1, ymax=3.5)
t2.rect2_DBH <- data.frame (xmin=0.5, xmax=1.5, ymin=0, ymax=0.85)
wi_DBH_ALL_Global$Scenario<-ifelse(wi_DBH_ALL_Global$Scenario2=="S1_AT_NP_AT_SP", "S1",
ifelse(wi_DBH_ALL_Global$Scenario2=="S4_UT_TA_NTS_NP", "S4",
ifelse(wi_DBH_ALL_Global$Scenario2=="S2_NTS_NP_NTS_SP", "S2",
ifelse(wi_DBH_ALL_Global$Scenario2=="S3_TU_NB_AT_NP", "S3", NA))))
wi_DBH_ALL_Global$Selection<-ifelse(wi_DBH_ALL_Global$Selection2=="Common", "Neutral",
ifelse(wi_DBH_ALL_Global$Selection2=="Selected", "Selected",
ifelse(wi_DBH_ALL_Global$Selection2=="Avoided", "Avoided", NA)))
## Plots for DBH selection
pd <- position_dodge(0.5)
Plot_wi_DBH_ALL_Global <- ggplot(wi_DBH_ALL_Global, aes(x=DBH, y=wi, colour=Scenario, group=Scenario)) +
geom_errorbar(aes(ymin=lower, ymax=upper),size=.6, width=.6, position=pd) +
geom_line(position=pd, size=0.6) +
geom_rect(data=t2.rect1_DBH, aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax), fill="gray20", alpha=0.1, inherit.aes = FALSE)+
geom_rect(data=t2.rect2_DBH, aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax), fill="gray20", alpha=0.1, inherit.aes = FALSE)+
geom_point(position=pd, size=1.5) + # 21 is filled circle
xlab("Diameter at breast height (cm)")+ geom_point(aes(shape=Selection), size=2) +
scale_shape_manual(values=c(6, 8, 17))+
ylab("Selection ratio (+/- CI)") +
scale_colour_hue(name="Scenario", # Legend label, use darker colors
l=50) +
expand_limits(y=-1) + # Expand y range
scale_colour_manual(values = c("#1B9E77", "#D95F02", "#E6AB02", "#666666", "#7570B3","#F8766D",
"#A6761D","#a6d854",
"#00BF7D","#E7298A"))+
geom_hline(aes(yintercept=1),size=0.6, colour="#7570B3", linetype="solid") +
theme(strip.text.x = element_text(size=8, face="plain"),
strip.text.y = element_text(size=8, face="plain"))+
theme(plot.title = element_text(face = "plain", size=9))+
theme_light()+
theme(legend.text = element_text(size = 8),
legend.title = element_text(size = 9),
#legend.position = "right", legend.direction = "vertical",
legend.position = c(.01, .99),
legend.justification = c("left", "top"),
legend.box.just = "left",
legend.margin = margin(6, 6, 6, 6))+
theme(plot.subtitle = element_text(vjust = 1),
plot.caption = element_text(vjust = 1),
axis.title = element_text(size = 9,
face = "bold"), axis.text = element_text(size = 8,
face = "bold"), axis.text.x = element_text(size = 8),
legend.text = element_text(size = 8),
legend.title = element_text(size = 9,
face = "bold"), legend.background = element_rect(size = 1.2))+
theme(strip.text.x = element_text(size=8, face="bold"),
strip.text.y = element_text(size=8, face="bold"), plot.subtitle = element_text(vjust = 1),
plot.caption = element_text(vjust = 1),
axis.title = element_text(size = 9,face = "bold"),
axis.text = element_text(size = 8, face = "plain", colour = "black"),
axis.text.x = element_text(size = 10),
legend.text = element_text(size = 8),
legend.title = element_text(size = 8,
face = "bold"), legend.background = element_rect(size = 1))+
theme(#legend.key.height= unit(2, 'cm'),
legend.key.width= unit(1.5, 'cm'))
## Saving the plot in different file formats (tif, pdf and eps)
setwd("C:/Users/tedon/OneDrive - UGent/PHD Luc/Articles en cours/Ecological factors of nesting site choice/Results")
tiff(file = "Plot_wi_DBH_ALL_Global4.tif", width = 17.5, height = 13, units = "cm", pointsize = 10,
bg = "white", res = 300)
Plot_wi_DBH_ALL_Global
dev.off()## png
## 2HAb_dir=setwd("C:/Users/tedon/OneDrive - UGent/PHD Luc/Articles en cours/Ecological factors of nesting site choice/Results")
ggsave("Plot_wi_DBH_ALL_Global4.pdf",Plot_wi_DBH_ALL_Global, path = HAb_dir,
scale = 1, width = 17.5, height = 13, units = "cm",
dpi = 600)
## Displaying the graph
"Plot_wi_DBH_ALL_Global"## [1] "Plot_wi_DBH_ALL_Global"
Figure 2: DBH representation in different scenarios of chimpanzee nesting and nesting site selection
S1: All trees in nesting plots vs all trees in systematic plots, S2: Nesting tree species in nesting plots vs nesting tree species in systematic plots, S3: Trees used for nest building vs all trees available in nesting plots, S4: trees used for nest building vs trees available of nesting tree species in nesting plots
Session info
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United Kingdom.1252
## [2] LC_CTYPE=English_United Kingdom.1252
## [3] LC_MONETARY=English_United Kingdom.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United Kingdom.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggmap_3.0.0 tmaptools_3.1-1 here_1.0.1
## [4] sf_1.0-5 forcats_0.5.1 stringr_1.4.0
## [7] dplyr_1.0.7 purrr_0.3.4 readr_2.1.2
## [10] tidyr_1.2.0 tibble_3.1.6 tidyverse_1.3.1
## [13] ExcelFunctionsR_0.1.4 ggrepel_0.9.1 rms_6.2-0
## [16] SparseM_1.81 Hmisc_4.6-0 Formula_1.2-4
## [19] survival_3.2-13 xlsxjars_0.6.1 rJava_1.0-6
## [22] viridis_0.6.2 viridisLite_0.4.0 vegan_2.5-7
## [25] lattice_0.20-45 rsq_2.2 reshape2_1.4.4
## [28] pvclust_2.2-0 phyloseq_1.38.0 optpart_3.0-3
## [31] plotrix_3.8-2 labdsv_2.0-1 mgcv_1.8-38
## [34] nlme_3.1-155 openxlsx_4.2.5 jtools_2.1.4
## [37] interactions_1.1.5 iNEXT_2.0.20 ggfortify_0.4.14
## [40] fpc_2.2-9 factoextra_1.0.7 DescTools_0.99.44
## [43] cowplot_1.1.1 corrr_0.4.3 cluster_2.1.2
## [46] clubSandwich_0.5.5 car_3.0-12 carData_3.0-5
## [49] aod_1.3.1 knitr_1.37 adehabitatHS_0.3.15
## [52] adehabitatHR_0.4.19 adehabitatLT_0.3.25 CircStats_0.2-6
## [55] boot_1.3-28 MASS_7.3-55 deldir_1.0-6
## [58] adehabitatMA_0.3.14 ade4_1.7-18 sp_1.4-6
## [61] questionr_0.7.7 png_0.1-7 plyr_1.8.6
## [64] NbClust_3.0 MuMIn_1.43.17 indicspecies_1.7.9
## [67] permute_0.9-7 ggthemes_4.2.4 ggpubr_0.4.0
## [70] ggplot2_3.3.5 ggdendro_0.1.22 dendextend_1.15.2
## [73] BiocStyle_2.22.0
##
## loaded via a namespace (and not attached):
## [1] prabclus_2.3-2 multcomp_1.4-18 data.table_1.14.2
## [4] rpart_4.1.16 RCurl_1.98-1.5 generics_0.1.2
## [7] BiocGenerics_0.40.0 TH.data_1.1-0 polspline_1.1.19
## [10] proxy_0.4-26 tzdb_0.2.0 xml2_1.3.3
## [13] lubridate_1.8.0 httpuv_1.6.5 assertthat_0.2.1
## [16] xfun_0.29 hms_1.1.1 jquerylib_0.1.4
## [19] evaluate_0.14 promises_1.2.0.1 DEoptimR_1.0-10
## [22] fansi_1.0.2 readxl_1.3.1 dbplyr_2.1.1
## [25] igraph_1.2.11 DBI_1.1.2 htmlwidgets_1.5.4
## [28] stats4_4.1.2 ellipsis_0.3.2 backports_1.4.1
## [31] bookdown_0.24 vctrs_0.3.8 Biobase_2.54.0
## [34] quantreg_5.87 abind_1.4-5 withr_2.4.3
## [37] robustbase_0.93-9 checkmate_2.0.0 mclust_5.4.9
## [40] ape_5.6-1 crayon_1.4.2 labeling_0.4.2
## [43] units_0.7-2 pkgconfig_2.0.3 GenomeInfoDb_1.30.1
## [46] nnet_7.3-17 rlang_1.0.1 diptest_0.76-0
## [49] lifecycle_1.0.1 miniUI_0.1.1.1 MatrixModels_0.5-0
## [52] sandwich_3.0-1 modelr_0.1.8 dichromat_2.0-0
## [55] rprojroot_2.0.2 cellranger_1.1.0 Matrix_1.4-0
## [58] Rhdf5lib_1.16.0 zoo_1.8-9 reprex_2.0.1
## [61] base64enc_0.1-3 rjson_0.2.21 rootSolve_1.8.2.3
## [64] bitops_1.0-7 KernSmooth_2.23-20 rhdf5filters_1.6.0
## [67] pander_0.6.4 Biostrings_2.62.0 classInt_0.4-3
## [70] jpeg_0.1-9 rstatix_0.7.0 S4Vectors_0.32.3
## [73] ggsignif_0.6.3 scales_1.1.1 magrittr_2.0.2
## [76] zlibbioc_1.40.0 compiler_4.1.2 RColorBrewer_1.1-2
## [79] lme4_1.1-27.1 cli_3.1.1 XVector_0.34.0
## [82] htmlTable_2.4.0 tidyselect_1.1.1 stringi_1.7.6
## [85] highr_0.9 yaml_2.2.2 latticeExtra_0.6-29
## [88] grid_4.1.2 sass_0.4.0 tools_4.1.2
## [91] lmom_2.8 RgoogleMaps_1.4.5.3 parallel_4.1.2
## [94] rstudioapi_0.13 foreach_1.5.2 foreign_0.8-82
## [97] roperators_1.1.0 gridExtra_2.3 gld_2.6.4
## [100] farver_2.1.0 Rtsne_0.15 stars_0.5-5
## [103] digest_0.6.29 BiocManager_1.30.16 shiny_1.7.1
## [106] Rcpp_1.0.8 broom_0.7.12 lwgeom_0.2-8
## [109] later_1.3.0 httr_1.4.2 kernlab_0.9-29
## [112] Deriv_4.1.3 colorspace_2.0-2 XML_3.99-0.8
## [115] rvest_1.0.2 fs_1.5.2 IRanges_2.28.0
## [118] splines_4.1.2 expm_0.999-6 multtest_2.50.0
## [121] Exact_3.1 flexmix_2.3-17 xtable_1.8-4
## [124] jsonlite_1.7.3 nloptr_2.0.0 modeltools_0.2-23
## [127] R6_2.5.1 pillar_1.7.0 htmltools_0.5.2
## [130] mime_0.12 glue_1.6.1 fastmap_1.1.0
## [133] minqa_1.2.4 class_7.3-20 codetools_0.2-18
## [136] mvtnorm_1.1-3 utf8_1.2.2 bslib_0.3.1
## [139] zip_2.2.0 rmarkdown_2.11 biomformat_1.22.0
## [142] munsell_0.5.0 e1071_1.7-9 rhdf5_2.38.0
## [145] GenomeInfoDbData_1.2.7 iterators_1.0.14 labelled_2.9.0
## [148] haven_2.4.3 gtable_0.3.0