Group_D_Project_Code_Main.Rmd

---
title: "R Notebook"
output: html_notebook
---
```{r setup}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = "G:/My Drive/UofT/2022-2023 FallWinter/EEB313/Project/Amineote life history/Data_Files")
getwd()
```
```{r}
setwd("G:/My Drive/UofT/2022-2023 FallWinter/EEB313/Project/Amineote life history/Data_Files")
getwd()
```
#Hypotheses draft

1. The reproductive output, in this case the average offspring produced per year, has a correlation with morphological traits. 

```{r}
library(tidyverse)
library(ggplot2)
library(ggfortify)
library(vegan)
library(dplyr)
library(lme4)
library(car)
library(MuMIn)
library(psych)
library(Hmisc)
library(gg.gap)
library(data.table)
```

#Data manipulations
```{r}
amine_avg <- read.csv("Amniote_Database_Aug_2015.csv") 
ncol(amine_avg)
count(amine_avg, order)
```

```{r}

aves_avg_bodysize <- amine_avg %>% #gather all parameters for class Aves 
  filter(class == "Aves") %>% 
  select(class, order, family, genus, species, adult_body_mass_g, longevity_y, litter_or_clutch_size_n, adult_svl_cm, egg_mass_g, litters_or_clutches_per_y,female_maturity_d,male_maturity_d, birth_or_hatching_weight_g, incubation_d, fledging_age_d) %>% 
  filter(adult_body_mass_g != -999.000, longevity_y != -999.000, litter_or_clutch_size_n != -999, adult_svl_cm != -999.00, egg_mass_g != -999.00, litters_or_clutches_per_y != -999.0,female_maturity_d != -999.000, male_maturity_d != -999.000,birth_or_hatching_weight_g != -999.00,incubation_d!= -999.00,fledging_age_d != -999.00,species != "lathami") %>% 
  mutate(avg_offspring_per_y = litter_or_clutch_size_n * litters_or_clutches_per_y)

write.csv(aves_avg_bodysize, "G:/My Drive/UofT/2022-2023 FallWinter/EEB313/Project/Aves_traits.csv", row.names=FALSE)
```

```{r}
mammal_avg_bodysize <- amine_avg %>% #gather all parameters for class Mammalia
  filter(class == "Mammalia") %>% 
  select(class, order, family, genus, species, adult_body_mass_g, longevity_y, litter_or_clutch_size_n,adult_svl_cm, birth_or_hatching_weight_g,birth_or_hatching_svl_cm,  litters_or_clutches_per_y, gestation_d, weaning_d, weaning_weight_g, female_maturity_d, male_maturity_d) %>% 
    filter(adult_body_mass_g != -999.000, longevity_y != -999.000, litter_or_clutch_size_n != -999, adult_svl_cm != -999.00, litters_or_clutches_per_y != -999.0, gestation_d != -999.000, weaning_d != -999.000, weaning_weight_g != -999.000, female_maturity_d != -999.000, birth_or_hatching_weight_g != -999.000, birth_or_hatching_svl_cm != -999.000, male_maturity_d != -999.000) %>% 
   mutate(avg_offspring_per_y = litter_or_clutch_size_n * litters_or_clutches_per_y)

write.csv(mammal_avg_bodysize, "G:/My Drive/UofT/2022-2023 FallWinter/EEB313/Project/Mammalia_traits.csv", row.names=FALSE)
```

```{r}
reptile_avg_bodysize <- amine_avg %>% #gather all parameters for class Reptilia
  filter(class == "Reptilia") %>% 
  select(class, order, family, genus, species, adult_body_mass_g, longevity_y, litter_or_clutch_size_n, adult_svl_cm, egg_mass_g,litters_or_clutches_per_y,female_maturity_d, male_maturity_d,birth_or_hatching_weight_g, incubation_d) %>% 
  filter(adult_body_mass_g != -999.000, longevity_y != -999.000, litter_or_clutch_size_n != -999.000, adult_svl_cm != -999.00,egg_mass_g != -999.00,litters_or_clutches_per_y != -999.0, female_maturity_d != -999.000, male_maturity_d != -999.000, birth_or_hatching_weight_g != -999.00, incubation_d != -999) %>% 
  mutate(avg_offspring_per_y = litter_or_clutch_size_n * litters_or_clutches_per_y)

write.csv(reptile_avg_bodysize, "G:/My Drive/UofT/2022-2023 FallWinter/EEB313/Project/Reptilia_traits.csv", row.names=FALSE)
```

################################################################################################
#Class: Aves

##Correlation test among variables
```{r}
corr.test(aves_avg_bodysize[,6:16], method = "pearson")

```

##Aves PCA without body-size correlated variables (egg mass, svl, birth weight)

###PCA
```{r}
# remove all the categorical variables so we can perform PCA
PCA_aves_no.corr <- aves_avg_bodysize %>%  
  select(adult_body_mass_g, longevity_y,  female_maturity_d, male_maturity_d, incubation_d, fledging_age_d)

#PCA analysis with Aves
aves.pca_no.corr <- prcomp(PCA_aves_no.corr, scale = T)
summary(aves.pca_no.corr)

#extract loadings from the PCA
aves.pca_no.corr$rotation 
```
####PC1 and PC2 explained the most variances. With body mass, female maturity time, male maturity time and fledging age contribute the most to PC1 (but also longevity and incubation period also contribute to PC1 at a level pretty close to 4 mains), and incubation time contribute the most to the PC2.

###plotting the PCA
```{r}
#plots of PCA of dataset without body-size correlated variables
biplot(aves.pca_no.corr)
autoplot(aves.pca_no.corr, loadings = T, loadings.label = T,loadings.colour = "red") +
  theme_classic()
plot(aves.pca_no.corr)
```
###Insert PC1 and PC2 back into dataset
```{r}
#Insert PC1 and PC2 back into the Aves dataset 
aves_avg_bodysize$PC1_no.corr <- aves.pca_no.corr$x[,1]
aves_avg_bodysize$PC2_no.corr <- aves.pca_no.corr$x[,2]
```
###Aves model with PC axes
```{r}
#model of Aves dataset with PC axes 
model_aves.pca_no.corr <- lm(data = aves_avg_bodysize, avg_offspring_per_y ~ PC1_no.corr + PC2_no.corr)
summary(model_aves.pca_no.corr)

```
####PC1 has a significant impact on the average offspring produced per year, but PC2 has a insignificant but relatively weak impact on the average offspring produced per year

###Ploting the model
```{r}
#plot of model with PC1 
plot_aves_no.corr_PC1 <- ggplot(data = aves_avg_bodysize, aes(x = PC1_no.corr, y = avg_offspring_per_y))+
  geom_point(size = 1, alpha = 0.5, colour = "red")+
  geom_smooth(se = F)+
  labs(x = "PC1", y = "average number of offsprings produced per year")+
  scale_x_continuous(breaks = seq(-5, 15, by = 1))+
  theme_classic()
plot_aves_no.corr_PC1 

#plot of model with PC2
plot_aves_no.corr_PC2 <- ggplot(data = aves_avg_bodysize, aes(x = PC2_no.corr, y = avg_offspring_per_y))+
  geom_point(size = 1, alpha = 0.5, colour = "blue")+
   geom_smooth(se = F)+
  labs(x = "PC2", y = "average number of offsprings produced per year")+
  scale_x_continuous(breaks = seq(-7, 10, by = 1))+
  theme_classic()
plot_aves_no.corr_PC2

```
#####the annual offspring produced is negatively correlated to PC1, in an rather linear mannar. The annual offspring produced will increase with increase PC2 value, reach a peak around PC2=0.5, and then decline rapidly.



###########################################################################################################

#Class: Mammalia

##Correlation test
```{r}
corr.test(mammal_avg_bodysize[,6:16], method = "pearson")
```

##Mammal PCA without highly correlated variables (birth weight/weaning weight/body mass, gestation and maturity, adult svl/birth svl)

###PCA
```{r}
# remove all the categorical variables so we can perform PCA
PCA_mammal_dataset_no.corr <- mammal_avg_bodysize %>% 
  select(adult_body_mass_g, longevity_y, adult_svl_cm, female_maturity_d, male_maturity_d, weaning_d)

#PCA analysis with Mammalia
mammal.pca_no.corr <- prcomp(PCA_mammal_dataset_no.corr, scale = T)
summary(mammal.pca_no.corr)

#extract loadings from the PCA
mammal.pca_no.corr$rotation
```
####PC1 and PC2 explained most variances. All varialbes contributed to PC1 relative equally, with body mass and weaning time contribute slightly less. Body mass contribute the most to the PC2 and it's negative

###ploting the result of PCA
```{r}
biplot(mammal.pca_no.corr)
autoplot(mammal.pca_no.corr, loadings = T, loadings.label = T,loadings.colour = "red") +
  theme_classic()
plot(mammal.pca_no.corr)
```

###Insert PC1 and PC2 back into the Mammal dataset
```{r}
mammal_avg_bodysize$PC1_no.corr <- mammal.pca_no.corr$x[,1]
mammal_avg_bodysize$PC2_no.corr <- mammal.pca_no.corr$x[,2]

```

###Mammalia model with PC axes
```{r}
model_mammal.pca_no.corr <- lm(data = mammal_avg_bodysize, avg_offspring_per_y ~ PC1_no.corr + PC2_no.corr)
summary(model_mammal.pca_no.corr)
```
####PC1 has a significant negative impact on the average offspring produced per year, PC2 has a insignificant but weak negative impact on the average offspring produced per year

###Ploting the mammal model
```{r}
plot_mammal_no.corr_PC1 <- ggplot(data = mammal_avg_bodysize, aes(x = PC1_no.corr, y = avg_offspring_per_y))+
  geom_point(size = 1, alpha = 0.5, colour = "red")+
  labs(x = "PC1", y = "average number of offsprings produced per year")+
  scale_x_continuous(breaks = seq(-5, 15, by = 1))+
  theme_classic()
plot_mammal_no.corr_PC1 

plot_mammal_no.corr_PC2 <- ggplot(data = mammal_avg_bodysize, aes(x = PC2_no.corr, y = avg_offspring_per_y))+
  geom_point(size = 1, alpha = 0.5, colour = "blue")+
  labs(x = "PC2", y = "average number of offsprings produced per year")+
  scale_x_continuous(breaks = seq(-15, 10, by = 1))+
  theme_classic()
plot_mammal_no.corr_PC2
```
#### The annual offspring produced is negatively correlated to PC1, it first experiences a rapid drop when PC1 is around 0 and then slowly decline. The annual offsprings produced concentrated in the region -1< PC2 <1, with most of species in this region and the value of offsprings produced per year scatter from 0 to the 30+, and with some other species scatter on the positive side of the PC2 axis and all with relative low annual offspring produced. 

#############################################################################################

#Class: Reptilia

##Correlation test
```{r}
corr.test(reptile_avg_bodysize[,6:15], method = "pearson")
```

##Reptilia PCA without highly correlated variables (body mass/svl/birth mass/egg mass)

###PCA
```{r}
# remove all the categorical variables so we can perform PCA
PCA_reptile_dataset_no.corr <- reptile_avg_bodysize %>% 
  select(adult_body_mass_g, longevity_y, female_maturity_d, male_maturity_d, incubation_d)

#PCA analysis with Aves
reptile.pca_no.corr <- prcomp(PCA_reptile_dataset_no.corr, scale = T)
summary(reptile.pca_no.corr)

#extract loadings from the PCA
reptile.pca_no.corr$rotation
```
####PC1 and PC2 explained most of the variance. female and male maturity time and adult body mass contribute the most to PC1, and longevity contribute the most to PC2

###Ploting the result of PCA
```{r}
biplot(reptile.pca_no.corr)
autoplot(reptile.pca_no.corr, loadings = T, loadings.label = T,loadings.colour = "red") +
  theme_classic()
plot(reptile.pca_no.corr)
```

###Insert PC1 and PC2 back into the Reptilia dataset
```{r}
reptile_avg_bodysize$PC1_no.corr <- reptile.pca_no.corr$x[,1]
reptile_avg_bodysize$PC2_no.corr <- reptile.pca_no.corr$x[,2]

```

###Reptilia model with PC axes
```{r}
model_reptile.pca_no.corr <- lm(data = reptile_avg_bodysize, avg_offspring_per_y ~ PC1_no.corr + PC2_no.corr)
summary(model_reptile.pca_no.corr)

```
####PC1 has a insignificant but weak impact on the average offspring produced per year, and PC2 has a significant impact on the average offspring produced per year.

###Ploting the reptile model
```{r}
plot_reptile_no.corr_PC1 <- ggplot(data = reptile_avg_bodysize, aes(x = PC1_no.corr, y = avg_offspring_per_y))+
  geom_point(size = 1, alpha = 0.5, colour = "red")+
  geom_smooth(se = F)+
  labs(x = "PC1", y = "average number of offsprings produced per year")+
  scale_x_continuous(breaks = seq(-5, 7, by = 1))+
  theme_classic()
plot_reptile_no.corr_PC1 

plot_reptile_no.corr_PC2 <- ggplot(data = reptile_avg_bodysize, aes(x = PC2_no.corr, y = avg_offspring_per_y))+
  geom_point(size = 1, alpha = 0.5, colour = "blue")+
  geom_smooth(se = F)+
  labs(x = "PC2", y = "average number of offsprings produced per year")+
  scale_x_continuous(breaks = seq(-3, 2, by = 0.2))+
  theme_classic()
plot_reptile_no.corr_PC2
```
####annual offspring produced is positively correlated to the PC1 and PC2.

####################################################################################