X_tropicalis_SDMs
Colin
March 10, 2020
Installing/ loading the packages
# install.packages('sf')
# install.packages('RStoolbox')
library(raster)
library(rgdal)
library(sf) # This package allows you to import shapefiles
library(RStoolbox) #This package allows you to perform PCA on rasters or raster stacks
library(ggplot2)
library(reshape2)DOWNLOAD the bioclimatic variables
This automatically downloads all 19 of the bioclim variables as a raster stack. You can download it in 3 resolutions: 2.5, 5, and 10. 10 returns grids that are about 18km x 18km. Obviously, we want more detail, but to save on computation time, we will start all models with the coarsest resolution
bioclim10 <- getData("worldclim",var="bio",res=10)Importing presence data
These presence data have already been cleaned
x_trop_pres <- read.csv('data/x_trop_cleaned_pres_df.csv')Importing the native range extent
Because we are using the native range to generate background data, we need to use the native range to ‘clip’ the bioclim raster. This is easy to do by creating an account in the IUCN redlist, and then just downloading a shapefile of the geographic range. Then I saved the shapefile in the working directory, and this command
x_trop_native_range_shapefile <- shapefile('data/x_trop_sf.shp')
projection(bioclim10)## [1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"projection(x_trop_native_range_shapefile)## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"Cropping the Bioclim raster
Here I am cropping the the bioclim raster stack based on the native range shapefile I imported above. Then I am plotting all the bioclim variables to demonstrate they are clipped in the extent of the imported shapefile
bio_10_native_range_crop_xtrop <- crop(bioclim10, extent(x_trop_native_range_shapefile))
plot(bio_10_native_range_crop_xtrop)
Running the PCA
This is a specific command used to perform PCA on a raster. The ‘spca = TRUE’ just sets it to use the correlation matrix rather than the covariance matrix. This is important because variables with higher values will otherwise be weighted more (e.g. if you have a variable that ranges from 0.1-0.9, and another variable that ranges from 100-1000, this latter variable will end up being much more heavily weighted in the PCA)
x_trop_ntv_rng_pca <- rasterPCA(bio_10_native_range_crop_xtrop, spca = TRUE)Looking at the model
This just pulls up how much variance in the bioclim variables is explained by each of the principal components
summary(x_trop_ntv_rng_pca$model)## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Standard deviation 3.1907181 1.9794760 1.440647 1.12437163 0.68022522
## Proportion of Variance 0.5358254 0.2062277 0.109235 0.06653745 0.02435297
## Cumulative Proportion 0.5358254 0.7420530 0.851288 0.91782544 0.94217841
## Comp.6 Comp.7 Comp.8 Comp.9
## Standard deviation 0.60938472 0.49051927 0.411684483 0.319054179
## Proportion of Variance 0.01954472 0.01266364 0.008920216 0.005357662
## Cumulative Proportion 0.96172313 0.97438677 0.983306985 0.988664647
## Comp.10 Comp.11 Comp.12 Comp.13
## Standard deviation 0.285988071 0.263409072 0.162886308 0.1111557576
## Proportion of Variance 0.004304694 0.003651807 0.001396418 0.0006502949
## Cumulative Proportion 0.992969340 0.996621147 0.998017566 0.9986678607
## Comp.14 Comp.15 Comp.16 Comp.17
## Standard deviation 0.1043206038 0.0828659588 0.0680617939 0.0438694309
## Proportion of Variance 0.0005727783 0.0003614088 0.0002438109 0.0001012909
## Cumulative Proportion 0.9992406390 0.9996020478 0.9998458588 0.9999471497
## Comp.18 Comp.19
## Standard deviation 3.168843e-02 1.233123e-08
## Proportion of Variance 5.285033e-05 8.003124e-18
## Cumulative Proportion 1.000000e+00 1.000000e+00As you can see, about 91% of the variance is explained by the first 4 PCs
Looking at the loadings
The loadings tell you the correlations of each PC to each of your variables. This means that when you see under 0.17222 under Bio1 of P.comp 1, this is telling you that the correlation coefficient between PC1 and Bio1 is 0.17222
x_trop_ntv_rng_pca$model$loadings##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9 Comp.10
## bio1 0.172 0.412 0.212
## bio2 0.263 -0.165 0.104 0.522 0.113 -0.362
## bio3 -0.280 -0.179 0.104 0.237 -0.732 0.295 0.213
## bio4 0.290 0.200 -0.287 0.239 0.287
## bio5 0.291 0.117 0.125 0.227 0.164
## bio6 -0.211 0.344 -0.107 -0.132 -0.226 0.167 0.148
## bio7 0.286 -0.135 0.132 0.129 0.259 -0.100
## bio8 0.161 0.353 0.309 -0.177 -0.207 0.509 -0.146 -0.188 -0.565
## bio9 0.468 -0.151 -0.192 -0.179 -0.755 0.286
## bio10 0.257 0.264 0.120 0.126 0.123 0.211 0.210
## bio11 0.466 0.391 -0.182 -0.181 0.194 0.391
## bio12 -0.268 0.331 0.101 0.107
## bio13 -0.169 0.564 -0.107 -0.254 -0.102
## bio14 -0.219 0.589 0.139 -0.386
## bio15 0.271 0.230 0.100 -0.375 -0.223 0.367
## bio16 -0.176 0.564 -0.187
## bio17 -0.241 0.533 -0.201
## bio18 -0.265 0.259 0.321 0.602 0.386 -0.128 0.439
## bio19 -0.226 0.262 -0.182 0.815 0.201 -0.239
## Comp.11 Comp.12 Comp.13 Comp.14 Comp.15 Comp.16 Comp.17 Comp.18 Comp.19
## bio1 0.662 0.462 0.226
## bio2 -0.135 0.104 -0.653
## bio3 0.127 0.117 0.180 -0.111 0.268
## bio4 0.400 0.190 0.107 -0.290 -0.112 -0.145 0.552
## bio5 0.130 -0.106 -0.472 -0.353 -0.142 0.450 -0.436
## bio6 0.126 -0.124 -0.521 -0.355 -0.265 0.453
## bio7 0.407 0.778
## bio8 -0.144
## bio9 0.127
## bio10 0.254 0.160 0.158 -0.116 -0.105 -0.760
## bio11 -0.297 0.198 -0.328 -0.236 0.245
## bio12 0.359 0.610 0.107 -0.233 0.259 -0.386 -0.110
## bio13 -0.634 0.239 -0.302
## bio14 -0.146 0.118 -0.378 0.481 -0.112
## bio15 -0.625 0.279 -0.203
## bio16 0.115 0.187 -0.166 -0.419 0.590
## bio17 -0.192 0.393 -0.561 -0.113 0.244
## bio18 -0.134
## bio19 -0.244
##
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9
## SS loadings 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## Proportion Var 0.053 0.053 0.053 0.053 0.053 0.053 0.053 0.053 0.053
## Cumulative Var 0.053 0.105 0.158 0.211 0.263 0.316 0.368 0.421 0.474
## Comp.10 Comp.11 Comp.12 Comp.13 Comp.14 Comp.15 Comp.16 Comp.17
## SS loadings 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## Proportion Var 0.053 0.053 0.053 0.053 0.053 0.053 0.053 0.053
## Cumulative Var 0.526 0.579 0.632 0.684 0.737 0.789 0.842 0.895
## Comp.18 Comp.19
## SS loadings 1.000 1.000
## Proportion Var 0.053 0.053
## Cumulative Var 0.947 1.000Creating a dataframe of the first 4 PCS
The purpose of this is to create a dataframe of the loadings on all bioclim variables, for the first four PCs. This way, we can export them, and link them to the bioclim variable descriptions. This will help us more intuitively interpret and describe what each PC seems to be capturing
first_pcs<-data.frame("Bioclim_variable" = c('BIO1', 'BIO2','BIO3','BIO4',
'BIO5','BIO6','BIO7','BIO8',
'BIO9','BIO10','BIO11','BIO12',
'BIO13','BIO14','BIO15','BIO16',
'BIO17','BIO18','BIO19'),
"PC1" = x_trop_ntv_rng_pca$model$loadings[1:19],
"PC2" = x_trop_ntv_rng_pca$model$loadings[20:38],
"PC3" = x_trop_ntv_rng_pca$model$loadings[39:57],
"PC4" = x_trop_ntv_rng_pca$model$loadings[58:76])
write.csv(first_pcs, 'data/first_pcs.csv')