Ejercicio guiado 10.3  (opcional). GBIF  y especies invasoras

Propósito del ejercicio

 

En este ejercicio final queremos recopilar y recordar algunos de los contenidos que hemos trabajado, para aplicarlo a un nuevo problema real, pero ligeramente simpificado. Vamos a intentar predecir la capacidad de invasión potencial de una especie de planta invasora catalogada en España.

Para ello, nuestros datos van a proceder de GBIF para demostrar el potencial que tiene el sistema para facilitar

El objetivo es que escojáis una especie diferente a la del ejemplo y hagáis una estimación parecida a la que se realiza en este ejercicio, introduciendo, si os sentís cómodos, las variantes que os parezcan en cuanto a análisis, visualización, etc. En el ejercicio trabajamos con una especie nativa de Australia. Si trabajáis por ejemplo con una de Sudáfrica cuando restrinjamos los datos, deberéis hacer cambios como "South Africa" en lugar de "Australia". Para simplificar el ejercicio recomendamos actuar como si una especie solo fuera nativa de un solo país aunque en realidad no sea así (a menos que queráis realizar más operaciones, en cuyo caso estamos disponibles para ayudaros).

El dataset invasoras.csv os da una lista de las plantas invasoras de la península ibérica. Para poder tener una lista más sencilla podemos hacer lo siguiente:

invasoras<-read.csv("invasoras.csv")
levels(as.factor(invasoras$aespecie))

Tened en cuenta que una especie sacada de esta lista tiene que corregirse y separar género y especie: En vez de ‘Oenothera_biennis’ por ejemplo, el nombre de la especie a utilizar debe ser ‘Oenothera biennis’.

Los datos del origen de cualquier especie lo podemos sacar bien de internet o bien de un manual técnico como es el Atlas de Plantas Invasoras. Así, podemos ver que especies como Oenothera biennis sería nativa de Canadá o Estados Unidos; Opuntia ficus indica, la chumbera, es nativa de México. Ailanthus altissima, de China. La estética Carpobrotus edulis, de Sudáfrica.

Vamos a realizar una serie de pasos para descargar en remoto desde un repositorio de datos de distribución de especies ciertamente conocido, la Global Biodiversity Information Facility (GBIF). Esta red ha recopilado la mayor cantidad de bases de datos de distribución de organismos, desde colecciones de museos y herbarios, colecciones particulares, artículos científicos, bases de datos anexas y observaciones de ciencia ciudadana en la naturaleza. Todo junto conforma un repositorio abierto donde podemos acceder a los datos de observaciones de cualquier especie del mundo. Muchas de estas observaciones contienen las coordenadas de la observación, por lo que nos permite hacernos un perfil de la distribución mundial de cualquier especie de planta, animal, hongo…

Por lo tanto, esta va a ser nuestra herramienta que vamos a poder acceder desde R y nuestro punto de partida.

También podéis acceder a los datos de GBIF desde está web: https://www.gbif.org/

 

Buscando una especie con R en GBIF

El primer paso es preparar los paquetes a cargar y que tendremos que instalar si no disponemos de ellos.

install.packages("raster","rgbif","tidyverse")

Carga de paquetes

packages<-c("raster","factoextra","FactoMineR","tidyverse","rgbif","pander")
sapply(packages,require,character.only=TRUE)
##     raster factoextra FactoMineR  tidyverse      rgbif     pander 
##       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE

Supongamos que la especie que nos interesa es Acacia dealbata de Australia porque es una peligrosa invasora de la península y queremos conocer su distribución dentro y fuera de su área de origen. Vamos a valorar si la especie crece en áreas similares a nivel climático en España y Australia, o bien al invadir este país ha cambiado su adaptación a las condiciones climáticas. Después, vamos a elaborar un mapa que nos permita predecir todos los sitios susceptibles de ser colonizados por esta invasora.

 

Descargando los datos de GBIF

Dentro de este paquete vamos a usar preferentemente la función occ_search() y varios argumentos para simplificar los resultados de lo que descarguemos:

O bien podemos usar la función occ_data(), de la cual tendremos que extraer la parte de los datos contenidos en la parte ‘data’ del objeto que nos genera. Esta función es recomendable si tenéis problemas con occ_search().

 

Descarga de los datos

#RGBIF_occ2<-occ_data(scientificName = "Acacia dealbata", hasCoordinate=TRUE,limit=1000)
RGBIF_occ2<-occ_search(scientificName = "Acacia dealbata", hasCoordinate=TRUE,limit=1000)
RGBIF_occ2<-RGBIF_occ2$data
RGBIF_occ2<-RGBIF_occ2[c("decimalLongitude", "decimalLatitude","country")]

Ahora vamos a utilizar la suite tidyverse para limpiar y seleccionar solo los datos que nos interesan de esta especie, aquellos de Australia (área nativa) y España (área invadida).

La función filter() nos selecciona aquellas filas que cumplen una condición, en este caso que el país sea o bien Australia o bien España. Es una función idéntica a subset(). Podéis ver dos sintaxis parecidas, la primera con las funciones del paquete dplyr de tidyverse.

Quitaremos los missing data para que no nos dé problemas con na.omit(). Os muestro aquí la operación tanto de la forma habitual como con el operador %>%.

GBIF_occ<-RGBIF_occ2%>%
   filter(country=="Spain" | country== "Australia")%>%na.omit()

La segunda con las funciones convencionales:

 GBIF_occ<-subset(RGBIF_occ2,country== "Spain" | country== "Australia")%>%na.omit()
 

Utilizando mapas climáticos mundiales

En esta sección vamos a utilizar el paquete raster para mapas de pixeles para utilizar información bioclimática de la base de datos mundial Worldclim (https://www.worldclim.org/). Tenéis de estos tres mapas que cubren todo el globo con extensión .tif:

Los tres mapas, de temperatura media anual (bio1, expresado en décimas de grado, por ejemplo, una temperatura media de 14.5ºC serán 145), rango de temperaturas (bio2) y precipitación anual total (bio12), se pueden manejar y gestionar con las funciones del paquete raster. La función genérica raster() es capaz de cargar este tipo de mapas a partir de archivos del disco duro, entre otras capacidades, y crear un objeto de tipo raster tal y como vimos en el ejercicio guiado de visualización de gráficos.

En el siguiente código os mostramos cómo también podemos ser capaces de crear objetos en R a partir de datos que están colgados en repositorios virtuales. Este código lo podéis encontrar, sin simplificar, en el respositorio de lenguaje y código abierto GitHub.

Podríais utilizarlo para descargar otros mapas climáticos que puedan ser relevantes para la distribución de las especies en lugar de las 3 que tenéis disponibles para el ejercicio. En total hay mapas de temperatura (media, máxima y mínima) así como precipitación, para cada mes. Como ese tipo de variables cartografiadas no son particularmente útiles para explicar la distribución de los seres vivos, se han calculado las variables bioclimáticas", como la temperatura media anual o la precipitación total que sí son más relevantes. En total hay 19 de estas variables, de las cuales hemos dejado disponibles en la carpeta datasets 3 que utilizaremos.

Si queréis explorarlas todas, están disponibles en la base de datos worldclim.

download.file(url = "http://biogeo.ucdavis.edu/data/climate/worldclim/1_4/grid/cur/bio_10m_esri.zip", 
              destfile = "Datasets/current_bioclim_10min.zip", 
              method = "auto")
  unzip( zipfile = "Datasets/current_bioclim_10min.zip", 
        exdir = "Datsets", 
        overwrite = T)

En realidad como ya comentábamos ya disponemos de los tres archivos en formato .tif en la carpeta de datasets, así que seguimos con ello.

Para ello vamos a cargarlos y a crear un objeto de rasterstack, que no es sino una simple compilación de los mapas que tengamos cargados, que coinciden en su extensión y en la resolución de los píxeles.

bio1<-raster("bio_1.tif")
bio2<-raster("bio_2.tif")
bio12<-raster("bio_12.tif")
 
bioclim_world <- stack(bio1,bio2,bio12)
class(bioclim_world)
## [1] "RasterStack"
## attr(,"package")
## [1] "raster"
 

Visualización de los mapas raster

La función plot(), una vez que los paquetes apropiados están instalados, va a ser capaz de cargar los gráficos requeridos, en este caso, mapas de píxeles, como ya visteis los que miraseis el material adicional de R como GIS del tema 5. Podemos visualizarlos todos, o solo uno de ellos. Podemos acceder a cada capa o mapa con el símbolo de $ como si accedieramos a una columna de un dataframe.

plot(bioclim_world)

plot(bioclim_world$bio_1)

Podemos visualizar la ocurrencia de los datos seleccionados de España y Australia sobre el mapa. Los coloreamos de manera diferente según el país, dato que tenemos que convertir a factor.

plot.new()
plot(bio1,main="Average temperature") 
points(x=GBIF_occ$decimalLongitude,y=GBIF_occ$decimalLatitude,col=as.factor(GBIF_occ$country),pch=16)

 

 

Extrayendo información de los mapas

A continuación vamos a utilizar la funcion extract() del paquete raster para cruzar los datos de las coordenadas descargadas de GBIF con los mapas, y extraer los valores de las tres variables y almacenarlo en un marco de datos. La función extract() aparece en más de un paquete, por lo que hay que especificarle que queremos la del paquete raster, de la siguiente manera: raster::extract(). El argumento df=TRUE especifica que queremos que nos devuelva los datos como si fuera un data.frame.

Veremos que nos devuelve los valores para cada coordenada, así que a continuación le pegaremos las columnas de los datos que ya teníamos y necesitamos (coordenadas y país).

bio_data<-GBIF_occ[1:2]%>%raster::extract(bioclim_world,.,df=TRUE)
 head(bio_data)
##   ID bio_1 bio_2 bio_12
## 1  1   152    70    716
## 2  2   132   119    920
## 3  3   154   126    499
## 4  4   154   126    499
## 5  5   155   129    456
## 6  6   134   102   1085
 bio_data<-cbind(bio_data,GBIF_occ)

Ahora que disponemos de un dataframe que nos permite tener la información climática de cada coordenada, podemos empezar a sacar algunos análisis, aunque aprovechamos para quitar algunas columnas no informativas:

bio_data$ID<-NULL 
summary(bio_data)
##      bio_1           bio_2            bio_12       decimalLongitude 
##  Min.   : 98.0   Min.   : 61.00   Min.   : 337.0   Min.   :-16.448  
##  1st Qu.:134.0   1st Qu.: 72.00   1st Qu.: 444.0   1st Qu.: -3.758  
##  Median :139.0   Median :102.00   Median : 613.0   Median : -2.134  
##  Mean   :142.2   Mean   : 94.06   Mean   : 676.2   Mean   : 23.188  
##  3rd Qu.:155.0   3rd Qu.:106.00   3rd Qu.: 810.0   3rd Qu.:  2.631  
##  Max.   :181.0   Max.   :148.00   Max.   :1704.0   Max.   :151.681  
##  decimalLatitude    country         
##  Min.   :-42.87   Length:225        
##  1st Qu.: 37.93   Class :character  
##  Median : 40.45   Mode  :character  
##  Mean   : 27.81                     
##  3rd Qu.: 41.56                     
##  Max.   : 43.54

Proponemos utilizar dos funciones del tidyverse para explorar las medias de las tres variables climáticas, mediante la concatenación de group_by(country) y summarise(). Vendrían a realizar la misma función que aggregate, que tenéis debajo. Pero os enseñamos esta forma tanto para practicar con el operador %>% como para mostraros group_by(). Es una función que genera grupos dentro de un dataframe de tipo tibble de tidyverse. Y es muy importante porque marca mucho el comportamiento de los datos. Todo lo que hagamos lo hará por grupos hasta que los desagrupemos.

str(bio_data)
## 'data.frame':    225 obs. of  6 variables:
##  $ bio_1           : num  152 132 154 154 155 134 110 110 133 151 ...
##  $ bio_2           : num  70 119 126 126 129 102 114 114 78 69 ...
##  $ bio_12          : num  716 920 499 499 456 ...
##  $ decimalLongitude: num  2.86 146.67 145.26 145.26 144.84 ...
##  $ decimalLatitude : num  41.7 -36.3 -36.2 -36.2 -36.1 ...
##  $ country         : chr  "Spain" "Australia" "Australia" "Australia" ...
bio_data%>%
  group_by(country)%>%
  summarise(bio1=mean(bio_1,na.rm=TRUE),bio2=mean(bio_2,na.rm=TRUE),bio12=mean(bio_12,na.rm=TRUE))
## # A tibble: 2 x 4
##   country    bio1  bio2 bio12
##   <chr>     <dbl> <dbl> <dbl>
## 1 Australia  130. 115.   955.
## 2 Spain      145.  89.8  620.
aggregate(bio_data,by=list(bio_data$country),FUN=mean,na.rm=TRUE)

Tanto en los datos numéricos como en los gráficos que visualizamos a continuación podemos anticipar que hay diferencias en las variables de temperatura (bio1 y bio2), pero quizás no en bio12 (precipitación anual). Vamos a generar tres gráficos con ggplot para ver si es verdad:

ggplot(bio_data,aes(x=country,y=bio_1))+stat_boxplot(geom="errorbar")+geom_boxplot()+ggtitle("Bio 1")

ggplot(bio_data,aes(x=bio_2,fill=country))+geom_density()+ggtitle("Bio 2")

ggplot(bio_data,aes(x=bio_12,fill=country,alpha=0.8))+geom_density()+ggtitle("Total rainfall")

 

Análisis de las diferencias en el nicho o preferencias climáticas de las especies

Los datos que disponemos nos permitirían hacer un análisis de las diferencias en el nicho climático, o preferencias climáticas, de la especie de estudio, en un país o en otro.

 

T de student o de Wald

Vamos a asumir(erróneamente) la normalidad de los datos y comparar las tres variables con un test de la t. Cargamos el paquete pander para obtener mejores tablas. Con `

t.test(bio_1~country,data=bio_data)%>%pander()

t.test(bio_2~country,data=bio_data)%>%pander()

t.test(bio_12~country,data=bio_data)%>%pander()
 

Análisis con un modelo lineal

En el modelo lineal asumimos también distribución normal en el test.

lm(bio_1~country,data=bio_data)%>%anova()%>%pander()
lm(bio_2~country,data=bio_data)%>%anova()%>%pander()
lm(bio_12~country,data=bio_data)%>%anova()%>%pander()

Los p-valores nos indican que hay un efecto estadísticamente significativo de las tres variables en las diferencias climáticas entre ambos territorios.

 

Análisis de componentes principales

En los análisis de cambio de nicho de las especies es muy habitual comparar en un análisis multidimensional como los PCA para analizar muchas variables ambientales a la vez. En el curso hemos visto tanto la función prcomp() y FactoMineR. Ahora vamos a utilizar una alternativa con el paquete ggfortifypara visualizar con ggplot, mediante la función autoplot().

ggfortify es una de las extensiones de ggplot para visualizar datos derivados de análisis multivariantes, como en este caso.

install.packages("ggfortify")

Para utilizarlo, vamos a realizar un análisis de componentes principales con PCA(). Utilizaremos como marco de datos las tres primeras columnas de nuestro dataset (donde tenemos alojados los datos climáticos), que no debe tener datos vacíos NA.

library(ggfortify)

pca_bio_data<-na.omit(bio_data)
pca<-PCA(pca_bio_data[1:3],graph=FALSE)

La función fviz() nos va a permitir una visualización más potente de nuestros datos derivados del análisis de componentes principales:

fviz_pca_ind(pca,
             label = "none", # hide individual labels
             habillage = as.factor(pca_bio_data$country), # color by groups
             palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE # Concentration ellipses
             )

 

Predicción en un mapa de la capacidad de invasión

Vamos a realizar en esta sección una versión sumamente simplificada de una técnica en biología llamada “modelos de distribución de especies” que es simplemente utilizar un modelo estadístico que toma como variable respuesta o dependiente la presencia o ausencia de una especie y como variables independientes, los datos ambientales que creemos que condicionan la distribución. Y coger ese modelo para proyectarlo sobre un mapa que contenga la información ambiental (en nuestro caso es exclusivamente climática) y tener un mapa de la distribución potencial de la especie, es decir que nos diga la idoneidad de cada zona para que esa especie prospere.

 

Preparando los datos

Vamos a restringir los datos a utilizar únicamente al área geográfica de España, porque aparentemente las preferencias climáticas han cambiado cuando la especie se convirtió en invasora aquí. Nos vamos a quedar solo con las columnas de las coordenadas porque tenemos que hacer más operaciones.

Acacia_Spain<-subset(pca_bio_data,country=="Spain")
Acacia_Spain<-Acacia_Spain[c("decimalLatitude","decimalLongitude")]

Y le vamos a añadir una columna, “presencia”, que simplemente refiere con el número 1 a que en estas coordenadas esta especie está presente en esas coordenadas. Voy a replicar el número 1 tantas veces como filas haya en mi dataset de coordenadas.

Acacia_Spain$presencia<-rep(1, times=nrow(Acacia_Spain))

Para construir un modelo estadístico que me permita hacer esta predicción, no solo necesito decirle dónde está mi especie, sino también donde no está. Vamos por tanto a generar datos aleatorios por toda la península que vamos a usar como ausencias, o dicho de otra manera, presencia=0. Pero antes, tengo que restringir mis mapas climáticos de todo el globo a solo a la península, así que vamos a recortarlos.

Las coordenadas aproximadas de los los límites los tengo ya disponibles en otro mapa de la península que ya tenía recortado. Lo voy a cargar:

limites<-raster("tmedia.tif")

Y con la función crop() del paquete raster vamos entonces a recortar nuestros 3 mapas climáticos:

bio1Spain<-crop(bio1,limites)
bio2Spain<-crop(bio2,limites)
bio12Spain<-crop(bio12,limites)
plot(bio2Spain)

Ahora solo tenemos que valernos de cualquiera de los tres para con la función sampleRandom() generar puntos con coordenadas aleatorias que vamos a utiliar como supuestas ausencias de Acacia en nuestro futuro modelo estadístico. El argumento xy=TRUE nos asegura que nos devuelva las coordenadas de los puntos. Vamos a generar tantos puntos de ausencia como puntos de presencia tuviéramos en nuestros datos extraídos.

ausencias<-sampleRandom(x=bio1Spain,size=nrow(Acacia_Spain),xy=TRUE)
  ausencias<-as.data.frame(ausencias)

Recomponemos el dataset para que tenga la misma estructura que nuestros datos de presencia que tenemos alojados en Acacia_Spain, quitando la columna de los datos de bio1Spain y añadiendo un valor 0 en la nueva columna presencia, y dándole los mismos nombres a las columnas que a Acacia_Spain

ausencias<-ausencias[c(1,2)]
ausencias$presencia<-rep(x=0,times=nrow(ausencias))
colnames(ausencias)<-colnames(Acacia_Spain)

Ahora podemos añadir los datos de supuestas ausencias a los datos de presencia:

Acacia_Spain<-rbind(Acacia_Spain,ausencias)

Y a continuación extraer los datos nuevamente de los tres mapas de España climáticos y a la información extraida añadirle los datos de presencia y ausencia de Acacia.

climaSpain<-stack(bio1Spain,bio2Spain,bio12Spain)
datosclima<-raster::extract(x=climaSpain,y=Acacia_Spain[1:2],df=TRUE)
datosclima<-cbind(datosclima,Acacia_Spain)
 

Construyendo un modelo predictivo

Con el dataframe conseguido ya podríamos realizar un modelo estadístico de cuál es el comportamiento de nuestra especie objetivo en España para construir luego un mapa. No podemos hacer una regresión logística tal cual; el problema en este caso es que no podemos pasar por alto que vamos a usar como variable respuesta unos datos que contienen o bien 0 (ausencia) o 1 (presencia). Por ello tenemos que recurrir a una familia de modelos estadística que es en esencia una extensión de las regresiones que aceptan otro tipo de distribuciones en los datos. Son los modelos generalizados, conocidos como GLMs.

La sintaxis es similar pero hay que especificarle la distribución que sigue nuestra variable respuesta. En el caso de respuestas que son 0/1 o verdadero/false, se trata de una distribución de poisson:

glm_acacia<-glm(presencia~bio_1+bio_2+bio_12,data=datosclima,family=poisson)

La función pander() me va a permitir visualizar los datos de cualquier tipo de modelo:

pander(glm_acacia)
Fitting generalized (poisson/log) linear model: presencia ~ bio_1 + bio_2 + bio_12
  Estimate Std. Error z value Pr(>|z|)
(Intercept) -27.3 509552 -5.358e-05 1
bio_1 2.854e-16 1493 1.911e-19 1
bio_2 -7.137e-16 3056 -2.335e-19 1
bio_12 -4.106e-17 197.3 -2.082e-19 1

Ya tengo mi modelo. Voy a utilizarlo para predecir la presencia o ausencia de una especie en mis mapas. Pero primero tengo que adecuar la información de los tres mapas a un dataframe que contenga las coordenadas. Es decir, convertir mi mapa en un dataframe que me diga pixel a pixel los valores de mis tres variables bioclimáticas. En este caso es fácil.

datosclima_dataframe<-as.data.frame(climaSpain,xy=TRUE)

La función predict() toma cualquier modelo estadístico, y dándole unos datos cuyos nombres se correspondan con las variables independientes o predictoras que se han utilizado para construirlo, nos devuelve una predicción de acuerdo a dicho modelo. Vamos a valernos de ella para hacer nuestra predicción:

prediccion_acacia<-predict(glm_acacia,datosclima_dataframe[3:5])

Ahora vamos a añadir a estos datos las coordenadas correspondientes a cada celda, sobreescribiendo el objeto prediccion_acacia:

prediccion_acacia<-cbind(datosclima_dataframe$x,datosclima_dataframe$y,prediccion_acacia)
colnames(prediccion_acacia)<-c("x","y","prediccion")

La función rasterfromXYZ() nos convertirá nuestro dataframe con las coordenadas y la predicción en un mapa. Vamos a usar la función scale() para centrar ligeramente los valores y que nos sirvan para interpretar mejor el mapa.

prediccion_acacia_mapa<-rasterFromXYZ(prediccion_acacia)
prediccion_acacia_mapa<-scale(prediccion_acacia_mapa)

Aun así, tened en cuenta que esto es solo una aproximación simplificada. Para estas predicciones, podéis usar paquetes como dismo o biomod2, con los cuales GBIF España trabaja en excelentes cursos presenciales.

 

Visualizando el mapa predictivo

Y ya podemos visualizar nuestro mapa de riesgo de invasión:

plot(prediccion_acacia_mapa)

Ya disponemos de nuestro mapa. Las zonas de mayor riesgo son las costas, más templadas, a menos que sean muy áridas, como Almería, o Alicante. Las zonas más continentales (muy secas y frías) o extremadamente frías, que es un dato que tenemos de nuestro mapa climático bio2, como las mesetas, el sistema central o Sierra Nevada, son las áeas de riesgo inexistente de acuerdo a este modelo.

Finalmente ¿podremos visualizar este mapa de otra manera? Por ejemplo, ¿con ggplot? La respuesta es como siempre afirmativa. Hay varias maneras pero vamos a ver una de ellas por su sencillez. En este caso cargando la extensión para análisis espaciales de rasterVis.

install.packages("rasterVis")
library(rasterVis)
 

En vez de usar la función normal ggplot() usamos la función ad hoc gplot() añadiendo la función geom_tile(). Integramos el parámetro estético fill= que equiparamos al valor de nuestro mapa:

gplot(prediccion_acacia_mapa) + geom_tile(aes(fill = value))

 

Programación de las funciones utilizadas

Esta sección es meramente informativa y la vamos a utilizar simplemente para explicar cómo podríamos utilizar lo que hemos hecho hasta ahora para poder automatizar este código para muchas especies.

Compilación de todo el código utilizado.

Vamos a compilar parte de las instrucciones en una sola función, que tome como argumentos el nombre científico de la especie invasora y el país del área nativa. Dentro de la función podéis ver el código para que ggplot cargue. Además metemos “prints” entre instrucciones para saber dónde está la función en cada momento.

mi_gbif <- function(nombre_cientifico, area_nativa) {
print("Buscando ocurrencias de GBIF")   
  
  RGBIF_occ2<-occ_search(scientificName = nombre_cientifico, hasCoordinate=TRUE,return="data",limit=15000,fields=c("decimalLatitude","decimalLongitude","country"))

     
  GBIF_occ<-RGBIF_occ2%>%
     filter(country=="Spain" | country== area_nativa)%>%na.omit()

  
print("Componiendo mapas")  
  
bio1<-raster("Datasets/bio_1.tif")
bio2<-raster("Datasets/bio_2.tif")
bio12<-raster("Datasets/bio_12.tif")
 
bioclim_world <- stack(bio1,bio2,bio12)

bio_data<-GBIF_occ[1:2]%>%raster::extract(bioclim_world,.,df=TRUE)
bio_data<-cbind(bio_data,GBIF_occ)

bio_data$ID<-NULL 

bio_data%>%
  group_by(country)%>%
  summarise(bio1=mean(bio_1,na.rm=TRUE),bio2=mean(bio_2,na.rm=TRUE),bio12=mean(bio_12,na.rm=TRUE))

##Atent@s al código de ggplot: primero lo creamos en un objeto y después, para poder plotearlo, llamaremos a ese objeto dentro de print()
print("Preparando plots")
plot1<-ggplot(bio_data,aes(x=country,y=bio_1))+stat_boxplot(geom="errorbar")+geom_boxplot()+ggtitle("Bio 1")
print(plot1)
plot2<-ggplot(bio_data,aes(x=bio_2,fill=country))+geom_density()+ggtitle("Bio 2")
print(plot2)
plot3<-ggplot(bio_data,aes(x=bio_12,fill=country,alpha=0.8))+geom_density()+ggtitle("Total rainfall")
print(plot3)

print("Preparando tests")
t.test(bio_1~country,data=bio_data)%>%pander()

t.test(bio_2~country,data=bio_data)%>%pander()

t.test(bio_12~country,data=bio_data)%>%pander()


pca_bio_data<-na.omit(bio_data)
pca<-PCA(pca_bio_data[1:3],graph=FALSE)

fviz_data<-fviz_pca_ind(pca,
             label = "none", # hide individual labels
             habillage = as.factor(pca_bio_data$country), # color by groups
             palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE # Concentration ellipses
             )
print(fviz_data)


}

vamos a probar la función con otra especie:

mi_gbif("Oenothera biennis","United States")

##Utilización de un bucle:

Para poder hacer varias especies podemos hacer un bucle:

latin<-c("Robinia pseudoacacia","Oenothera biennis","Eucaliptus globulus","Carpobrotus edulis")
nombre<-c("Falsa acacia","Oenotera","Eucalipto","Cangrejo rojo")
origen<-c("United States of America","United States of America","Australia","South Africa")
especies<-cbind(latin, nombre, origen)%>%
  data.frame()
colnames(especies)<-c("Especie","Nombre común","Origen")

Generamos para correr el bucle, dos objetos temporales a partir de los datos que nos dice el data.frame creado.

for (i in 1:nrow(especies)){
  especie<-especies[i,1]
  nativa<-especies[i,3]
  print(paste("Extrayendo datos de la especie",especie))
  mi_gbif(especie,nativa)
}