knitr::opts_chunk$set(
    echo = TRUE,
    message = TRUE,
    warning = TRUE
)

1 Introducción

1.1 Justificación del curso

Cursos de R online y libros hay a patadas (Ignasi Bartomeus).
Equivocarnos cuanto más mejor: PICAD CÓDIGO!
GBIF como eje central de la formación para la biodiversidad.

1.2 Objetivos

Entender las ventajas de usar R (o otros lenguages de programación)
Que descubrais que con R se puede hacer casi todo (desde Sudokus hasta pedir pizza)
Saber suficiente R para poder “googlear” lo que necesiteis aprender/resolver a partir de ahora.
Sentar las bases para un uso eficiente de R:
Ahorrar tiempo escribiendo código
Trabajar de manera coherente con los datos
Trabajar de manera reproducible sobre todo con vosotros mismos
No perder el tiempo ni escribiendo ni arrancando el programa.
Ir de los cimientos hasta donde lleguemos, para no perder tiempo después.

1.3 Lo básico

1.3.1 Qué es R

Enlace a R CRAN

Programa nacido en 1993 a partir del código S.
Un programa que es intérprete, lenguaje y entorno de programación a la vez.
10.000 paquetes.
Orientado a objeto = el objeto con el que trabajas define qué función le aplicas.
Carga los datos en la RAM.
Otros programas relacionados: SAS, Python, Julia, Rubi

The 9 Best Languages For Crunching Data

Popularity of data science software

Programas de los que parte: S, Fortran, C++.
Interés de Microsoft…
… y el tidyverse. R como lenguaje que necesita otra fachada.

Hadley Wickham

1.3.2 Por qué R

Gran capacidad para la manipulación de los datos y acabar produciendo la información necesaria sobre los datos.
Comunidad de construcción de utilidades diversas a través de paquetes:
GIS
Estadística
Big data
Web scraping
Genética
Comunidad de apoyo a la programación.
Reproducibilidad.
Conexión con otros lenguajes
Edición directa de figuras

1.3.3 Limitaciones

RAM demandante, cada vez menos.
Curva de aprendizaje lenta.
Inconsistencia en el código.
Calidad de los paquetes muy variable y no garantista.
Documentación fragmentada.
Los errores de programación no siempre son fáciles de encontrar.
¿Qué pasa con los grandes datasets?
No es lento, pero no es el más rápido.

1.4 Recursos de aprendizaje y análisis

StackOverflow
How do I ask a good question?
Google (e.g. error message + r)
Listas de correo…!
España:
Grupo de ecoinformática de la AEET
Redes R locales: Madrid, Sevilla..
Datacamp
Libro de recetas
R for data science
Advanced R
Cheatsheets
R-Bloggers

1.5 Recursos para la biodiversidad: una lista

Genómica, expresión genética, NGS: Bioconductor

ADE4: Ecología y análisis multivariante, comunidades

Vegan: Ecología y análisis multivariante,comunidades

Análisis filogenético: Paquetes geiger, ape, ggtree,adegenet….polysat….

Modelización de nicho: biomod2, dismo, ENMtools

Mapas: shapefiles, raster, RGDAL, maptools, rasterVis (3D)…

Trabajo en paralelo: snow, parallel, foreach

Tidyverse

Análisis estadístico: lmerTest, learnBayes…….

2 Introducción a R y RStudio

2.1 Conceptos básicos del entorno R

R consta de una consola y una ventana que podemos abrir para escribir y ejecutar los scripts, líneas de texto con las órdenes del programa que podemos guardar, reproducir…

Lostres primeros conceptos básicos:

Cada línea de código que se envía a la consola o se escribe directamente es una instrucción, a menos que haya algún elemento de la sintaxis que deje al programa esperando.
Para crear un nuevo objeto usaremos el comando operador <-, no =, por varias razones. Por ejemplo: mi_edad<-36
Para llamar a una función de R se escribe la función y los argumentos o parámetros que queremos, dentro de los paréntesis: sqrt(mi_edad).

packages<-c("rbenchmark","rgbif","ALA4R","pander","parallel","pander","tidyverse","readr")
sapply(packages,require,character.only=T)

## Loading required package: rbenchmark

## Loading required package: rgbif

## Loading required package: ALA4R

## Loading required package: pander

## Loading required package: parallel

## Loading required package: tidyverse

## -- Attaching packages -------------------------------------------------------------------- tidyverse 1.2.1 --

## v ggplot2 2.2.1     v purrr   0.2.4
## v tibble  1.4.2     v dplyr   0.7.4
## v tidyr   0.8.0     v stringr 1.3.0
## v readr   1.1.1     v forcats 0.3.0

## -- Conflicts ----------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

## rbenchmark      rgbif      ALA4R     pander   parallel     pander 
##       TRUE       TRUE       TRUE       TRUE       TRUE       TRUE 
##  tidyverse      readr 
##       TRUE       TRUE

2.2 Uso de GBIF en un entorno R

En este ejercicio vamos a utilizar R para hacer mineria de datos de e ir conociendo el programa Ya os han presentado la gran cantidad de datos de la plataforma. Para ir familiarizándonos con R vamos a realizar un ejercicio de cómo sacar partido a la información de GBIF de manera eficiente.

2.2.1 Gestión básica de carpetas, archivos y directorios

Partimos de la creación de nuestro propio directorio desde las funcionalidades de RStudio. Pero necesitamos primero familizarizarnos con algunos comandos para saber dónde estamos y qué hay en nuestra carpeta: getwd() nos da la ruta al directorio;setwd() nos permite cambiar el directorio. Y dir(), nos permite saber qué hay dentro.

getwd()

dir()

R devuelve una lista de archivos que están en la carpeta, entre comillas porque son nombres, cadenas de texto. Ahora podríamos cargar el archivo “especies.txt”. Lo cargamos con la función read.table().(especies<-read.table("especies.txt")). Un truco para ver tu objeto al mismo tiempo que lo asignas es meterlo entre parénteris. Eso no solo crea el objeto sino que lo devuelve también en pantalla. ¿Qué pasa?

R devuelve un error, no parece tener 5 elementos. Volvemos a cargarlo, pero esta vez especificando el tipo de carácter que nos separa las columnas.

(especies<-read.table("especies.txt",sep=";"))

Ahora el problema es que la primera fila contiene los nombres de las varaibles pero no han sido codificados como tal.

(especies<-read.table("especies.txt",sep=";",header=TRUE))

Ahora sí. Qué tipo de variable es especies?

class(especies)

¿Y algunas de las variables?

class(especies$name)
class(especies[2])
class(especies[,2])

class(especies[3:5])
class(especies[,5])
dim(especies)
nrow(especies)
ncol(especies)
seq_along(especies)
class(seq_along(especies))
class(1)
class("1")

2.2.2 Buscando una especie con R en GBIF

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. Esto es lo que vamos a sacar de este ejercicio. Vamos a usar dos paquetes:

ALA4 R, Atlas of Living Australia R, es un paquete para descargar los datos del Atlas Australiano, como parte de los datos de GBIF.

rgbif para sacar los datos de distribución de nuestra especie objetivo en la base de datos de gbif. Lo primero que necesitamos hacer es instalar la librería con el comando install.packages(), que instalará ese paquete en el directorio del programa R por defecto y todos los paquetes que debemos instalar que necesitamos. También lo podemos hacer desde Rstudio y es una opción quizás más recomendable en la mayoría de los casos.

library(ALA4R)

#manera 1 de seleccionar los datos
my_especie<-subset(especies,invader ==  TRUE)
my_especie<-subset(my_especie,phylum != "Animal")
my_especie<-subset(my_especie,country == "Australia")
my_especie<-my_especie$name
my_especie<-droplevels(my_especie)

#manera2
my_especie<-subset(especies,invader == TRUE & phylum != "Animal" & country == "Australia")
my_especie<-my_especie$name
my_especie<-droplevels(my_especie)

#manera3

my_especie<-especies[4,1]
my_especie
my_especie<-droplevels(my_especie)

Ya tenemos la especie seleccionada.

ALA4R nos exige que demos la razón por la cual queremos descargar los datos. Para ello, vamos a ver cuáles son esas razones. De ellas, elegiremos la 3.

ala_reasons()
ala_reasons()%>%class()

ALA_occ<-occurrences(taxon=my_especie,
            download_reason_id=3)

class(ALA_occ)
ALA_occ$data%>%names()
ALA_occ$data%>%class()
ALA_occ$data%>%dim()
ALA_occ<-subset(ALA_occ$data,scientificName=="Acacia dealbata")
ALA_occ%>%dim()
ALA_occ<-ALA_occ[22:23]
colnames(ALA_occ)<-c("y","x")


#plot(aus,main=paste("Distribución de",my_especie[1]))

Se puede añadir la capa de puntos debajo con: points(x=ALA_occ$x,y=ALA_occ$y)

2.2.2.1 RGBIF

Dentro de este paquete vamos a usar la función occ_search() y sus argumentos.

library(rgbif)
RGBIF_occ0<-occ_search(scientificName = "Acacia dealbata")
 str(RGBIF_occ0)
 RGBIF_occ1<-occ_search(scientificName = "Acacia dealbata", hasCoordinate=TRUE,return="data")
 str(RGBIF_occ1)
 dim(RGBIF_occ1)
 class(RGBIF_occ1)
 RGBIF_occ1
 
 RGBIF_occ2<-occ_search(scientificName = "Acacia dealbata", hasCoordinate=TRUE,return="data",limit=15000)
 names(RGBIF_occ2)
 RGBIF_occ2[3]
 RGBIF_occ2[4]
 RGBIF_occ2[48]
 RGBIF_occ2["country"]
 
 GBIF_occ<-RGBIF_occ2%>%select(decimalLatitude,decimalLongitude,country)
 ggplot(GBIF_occ,aes(x=decimalLongitude,y=decimalLatitude))+geom_point()
 plot(aus)
 points(x=GBIF_occ$decimalLongitude,GBIF_occ$decimalLatitude)
 points(x=ALA_occ$x,y=ALA_occ$y,col="red")
 
 GBIF_occ<-RGBIF_occ2%>%
   select(decimalLatitude,decimalLongitude,country)%>%
   filter(country=="Spain" | country== "Australia")
 
 GBIF_occ<-subset(RGBIF_occ2,country== "Spain" | country== "Australia")
 GBIF_occ<-GBIF_occ[c(3,4,48)]

#bajando las capas de worldclim

dir.create("WorldClim_data", showWarnings = F)
 #curent bioclim
 download.file(url = "http://biogeo.ucdavis.edu/data/climate/worldclim/1_4/grid/cur/bio_10m_esri.zip",
               destfile = "WorldClim_data/current_bioclim_10min.zip",
               method = "auto")
 unzip( zipfile = "WorldClim_data/current_bioclim_10min.zip",
        exdir = "WorldClim_data/current",
        overwrite = T)
 list.files("WorldClim_data/current/bio/")
 bioclim_world <- stack(list.files("WorldClim_data/current/bio", pattern ="bio_", full.names=T), RAT = FALSE)
 plot(bioclim_world)
 class(bioclim_world)
 bioclim_world$bio_1
 plot(bioclim_world$bio_12)
 points(x=GBIF_occ$decimalLongitude,y=GBIF_occ$decimalLatitude,col=GBIF_occ$country)no funciona porque no es un factor
 points(x=GBIF_occ$decimalLongitude,y=GBIF_occ$decimalLatitude,col=as.factor(GBIF_occ$country)) convertimos a factor

#extracción de los datos:
 bio_data<-GBIF_occ[1:2]%>%raster::extract(bioclim_world,.,df=TRUE)#coordenadas invertidas
 head(bio_data)
 bio_data<-GBIF_occ[c(2,1)]%>%raster::extract(bioclim_world,.,df=TRUE)#coordenadas arregladas
 bio_data<-cbind(bio_data,GBIF_occ)

 algunas estadísticas
 summary(bio_data)
 colMeans(bio_data[2:ncol(bio_data)-1])eliminamos el ID y las columnas no numéricas
 bio_data$ID<-NULL
 ?aggregate
 aggregate(bio_data,by=list(bio_data$country),FUN=mean)sin usar fórmula
 boxplot(bio_1~country,data=bio_data)
 y si tuviéramos que hacerlo para todos los casos..?
 t.test(bio_1~country,data=bio_data)
 hist(bio_data$bio_1,main="Histograma")

 install.packages("factoextra")
 library(factoextra)
 library(FactoMineR)
 pca_adealbata<-PCA(bio_data[2:19])
 fviz_pca(pca_adealbata,col.ind=bio_data$country)

3 Visualización de datos con R

GGPLOT2

3.1 Conceptos básicos de la visualización con R

GGPlot(grammar of graphics) es una aproximación a la visualización de datos implementada en el paquete ggplot2 que forma parte del tidyverse. Tanto el tidyverse en general como ggplot tienen filosofías de programación algo más novedosas y que tienen como intención hacer un flujo de trabajo más fluido.

Vamos a explorar primero el concepto de datos “limpios”:

3.1.1 Filosofía tidy data

(ghg_tidy<-read_csv("GHG.csv"))
(ghg_tidy<-read_csv2("GHG.csv"))#usamos tidyverse para subir los datos. Tendríamos que utilizar parse_double para mejorarlo...o algo parecido:type convert

(ghg<-read.csv("GHG.csv",sep=";",header=TRUE))
ghg_tidy<-read_csv2("GHG.csv")%>%
  type_convert()

ghg_tidy<-ghg_tidy%>%gather(3:10,key="year",value="ghg")
ghg_tidy
filter(ghg_tidy,Country=="Algeria")
ghg_tidy$year<-parse_date(ghg_tidy$year,format="%Y")  
ghg_tidy

3.1.2 Tidy data aplicada a estadística

lm_ghg<-lm(ghg~Continent+year,data=ghg_tidy)
summary(lm_ghg)
anova(lm_ghg)

Obviamos el continente para intentar ver el efecto del tiempo:

lm_ghg2<-lm(ghg~year,data=ghg_tidy)
summary(lm_ghg2)
anova(lm_ghg2)

Volvemos a incorporar el Continente pero como factor aleatorio en un modelo mixto. Para ello es mejor considerar el contintente como el factor que verdaderamente es.

ghg_tidy$Continent<-parse_factor(ghg_tidy$Continent,levels=ghg_tidy$Continent)

library(lmerTest)
lmer_ghg<-lmerTest::lmer(ghg~year+(1|Continent),data=ghg_tidy)
summary(lmer_ghg)

¿Cómo podemos plantearnos hasta dónde queremos estirar los datos?

3.2 Creación de figuras de calidad con ggplot

grammar of graphics

data(iris)
head(iris)

##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

plot(iris$Sepal.Length,iris$Petal.Length)

plot(iris$Sepal.Length,iris$Petal.Length,col=as.factor(iris$Species))

iris_plot<-ggplot(iris,aes(x=Sepal.Length,y=Petal.Length))
iris_plot+geom_point()
iris_plot+geom_point(aes(col=Species))
iris_plot+geom_point(aes(col=Species,size=Sepal.Width))

#manejando el tamaño según un parámetro
iris_plot+geom_point(aes(shape=Species,size=Sepal.Width))

#escala continua de color

iris_plot+geom_point(aes(shape=Species,col=Sepal.Width))+
  scale_colour_gradient(low = "green", high = "red")#cambiamos en esta línea la escala por defecto

#con más parámetros
iris_plot+geom_point(aes(col=Species))+
  xlab("Longitud del sépalo")+
  ylab("Longitud del pétalo")+
  theme_bw()+#atajo al tema blanco y negro
  ggtitle("Relación morfométrica en Iris")+
  geom_smooth()+
  geom_smooth(method="lm")


#volvemos a GHG como ejemplo de mapeo estadístico
ggplot(ghg_tidy,aes(x=year,y=ghg))+geom_point()
ggplot(ghg_tidy,aes(x=year,y=ghg))+geom_boxplot()
ggplot(ghg_tidy,aes(x=year,y=ghg,group=Continent))+geom_boxplot()
ggplot(ghg_tidy,aes(x=Continent,y=ghg,fill=Continent))+
  geom_boxplot()+
  facet_grid(.~year)

ghg_NAEU<-filter(ghg_tidy,Continent=="North America" | Continent == "Europe")
ggplot(ghg_NAEU,aes(x=Continent,y=ghg,fill=Continent))+
  geom_boxplot()+
  facet_grid(.~year)

ggplot(ghg_NAEU,aes(x=year,y=ghg,group=Continent))+
  geom_point()+
  geom_smooth()

ggplot(ghg_NAEU,aes(x=year,y=ghg,col=Continent))+
  geom_point()+
  geom_smooth()


best_in_class <- ghg_NAEU %>%
  group_by(Continent) %>%
  filter(row_number(desc(ghg)) == 1)

ggplot(ghg_NAEU,aes(x=year,y=ghg,col=Continent))+
  geom_point()+
  geom_smooth()+
  geom_point(size = 3, shape = 1, data = best_in_class)

library(ggrepel)

ggplot(ghg_NAEU,aes(x=year,y=ghg,col=Continent))+
  geom_point()+
  geom_smooth()+
  geom_point(size = 3, shape = 1, data = best_in_class)+
  ggrepel::geom_label_repel(aes(label = Country), data = best_in_class)

#ejemplos con RGBIF
head(bio_data)
ggplot(bio_data,aes(x=bio_1,col=country))+geom_histogram()
ggplot(bio_data,aes(x=bio_1,fill=country))+geom_histogram(bins=45)
ggplot(bio_data,aes(x=bio_1,col=country))+geom_density()
ggplot(bio_data,aes(x=bio_1,fill=country))+geom_density(aes(alpha=2/3))

ggplot(bio_data,aes(x=decimalLongitude, y=decimalLatitude))+geom_point(aes(col=country,size=bio_1))
ggplot(bio_data,aes(x=decimalLongitude, y=decimalLatitude))+geom_point(aes(shape=country,size=bio_1))
class(aus)
ggplot(as.data.frame(aus,xy=TRUE),aes(x,y))+geom_raster(z)

#quitando la leyenda
ggplot(bio_data,aes(y=bio_1,x=country))+stat_boxplot(geom="errorbar")+
  geom_boxplot(aes(fill=country))+
  scale_fill_discrete(guide=FALSE)#esta línea de código quita la leyenda
  
#cambiar el orden de las apariciones
class(bio_data$country) #nos aseguramos de que es un factor: no lo es
bio_data$country<-factor(bio_data$country,levels=c("Spain","Australia"))

#Mapping variable values to colors
ggplot(bio_data,aes(y=bio_1,x=country))+stat_boxplot(geom="errorbar")+
  geom_boxplot(aes(fill=country))+
  #scale_fill_discrete(guide=FALSE)+
  scale_fill_manual(values=c("#CC6666", "#9999CC"))

#diagramas de barras: emisiones per capita promedio por continente en 2015

ghg_sum<-ghg_tidy%>%filter(year=="2015-01-01" & Continent != "NA")%>% #filtramos y solo nos quedamos los de 2015
  group_by(Continent)%>%
  summarize(ghg=mean(ghg))

ggplot(ghg_sum,aes(x=Continent,y=ghg))+geom_bar(stat="identity")+
  theme_classic()

ghg_tidy%>%filter(year=="2015-01-01" & Continent == "Asia")%>%  #filtramos y solo nos quedamos los de 2015 de Asia
  #arrange(desc(ghg)) #para usar tras el primer intento
  ggplot(aes(x=Country,y=ghg))+
  geom_bar(stat="identity")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
  
ghg_tidy%>%filter(year=="2015-01-01" & Continent == "Asia")%>%  #filtramos y solo nos quedamos los de 2015 de Asia
  arrange(asc(ghg)) #para usar tras el primer intento

subset(ghg_tidy,Country=="China")%>%
  ggplot(aes(x=year,y=ghg))+
  geom_line()

subset(ghg_tidy,Country=="Germany" | Country=="Spain and Andorra")%>%
  ggplot(aes(x=year,y=ghg))+
  geom_line(aes(col=Country))+ 
  scale_colour_discrete(name  ="País",labels=c("Alemania", "España"))

¿Cómo podemos añadir mapas? Este es un ejemplo de funcionalidad de ggplot y la familia de extensiones:

#Sobre coordenadas

library(ggmap)

mapImageData1 <- get_map(location = c(lon = -3.6030016179, lat = 40.348425),
                         color = "color",
                         source = "google",
                         maptype = "satellite",
                         zoom = 6)
ggmap(mapImageData1)+
  geom_point(data=bio_data,aes(x=decimalLongitude,y=decimalLatitude,col="yellow"))

library(ggmap)
aqui<-geocode("Jardin Botanico,Madrid,Spain")

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Jardin%20Botanico,Madrid,Spain&sensor=false

aqui

##         lon      lat
## 1 -3.691127 40.41111

ggmap(get_map(aqui),zoom=4)+geom_point(data=aqui,aes(lon,lat),size=4)

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=40.411106,-3.691127&zoom=10&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

aqui<-geocode(“Jardin Botanico,Madrid,Spain”) aqui ggmap(get_map(aqui),zoom=4)+geom_point(data=aqui,aes(lon,lat),size=4)

3.2.1 Recursos para optimizar ggplot

Referencia oficial

R cookbook para gráficos

Consejos

Chuleta-cheatsheet

Color

3.2.2 Paquetes adicionales

Chuleta-cheatsheet

GGTREE para filogenias

Visualización de mapas raster

4 Programación con R

4.1 Introducción a la programación eficiente

La regla de oro propuesta es automatizar mediante funciones o programación todo aquello que haya que escribir más de 2 veces, por una cuestión de eficiencia: se pueden cometer errores copiando y pegando y consume tiempo. Para ello existen funciones pensadas como atajo y programar bucles (loops) que hacen la misma operación reiterada. Es necesario conocerlos, pero no siempre son necesarios. Ya hay funciones que hacen aquello por nosotros.

Las dos principales opciones que veremos son for y while

for va a recorrer los datos de la manera que le digamos:

for (i in 1:10){
  print (i*10)
}

especies

for (s in especies$name){
  message(paste("La especie es",s),appendLF=TRUE)
  }

for (s in especies$name){
  print(paste("La especie es",s))
}

for (s in especies$name){
  cat(paste("La especie es",s))
}

#jugando con los índices
for (i in 1:nrow(especies)){
  especie<-especies[i,1]
  ncomun<-especies[i,3]
  print(paste(especie,":",ncomun))
}

#recorriendo columnas

for (i in seq_along(bio_data[1:21])){#seq_along previene números extraños; pero hay datos no numéricos, de ahí el 1
  class(bio_data[i])%>%print()
  print(mean(bio_data[,i]))
}
#algo más eficiente:
colMeans(bio_data[1:21])
colMeans(bio_data[1:21])%>%class()
#nos hemos dejado la manera de poder almacenar nuestros resultados:

mis_medias<-c() #aquí vamos a tener un problema de engorde de vectores
for (i in seq_along(bio_data[1:21])){#seq_along previene números extraños; pero hay datos no numéricos, de ahí el 1
  #print(mean(bio_data[,i])) #ya no nos interesa visualizar sino almacenar
  media<-mean(bio_data[,i])
  mis_medias[i]<-media
}
mis_medias

#sin problemas de engorde de vectores:
mis_medias<-c(rep(NA,ncol(bio_data[1:19]))) #anticipamos la longitud del vector
for (i in seq_along(bio_data[1:21])){#seq_along previene números extraños; pero hay datos no numéricos, de ahí el 1
  #print(mean(bio_data[,i])) #ya no nos interesa visualizar sino almacenar
  media<-mean(bio_data[,i])
  mis_medias[i]<-media
}

#hacemos un dataframe para almacenar más variables
mis_medias<-data.frame(variable=character(),media=double()) #o bien sin engorde de vectores
mis_medias<-matrix(NA,ncol=2,nrow=ncol(bio_data)-1)%>%as.data.frame()
head(mis_medias)
colnames(mis_medias)<-c("variable","media")
for (i in seq_along(bio_data[1:21])){#seq_along previene números extraños; pero hay datos no numéricos, de ahí el 1
  #print(mean(bio_data[,i])) #ya no nos interesa visualizar sino almacenar
  mis_medias[i,1]<-colnames(bio_data[i]) #coloco el nombre de la variable en su celda
  mis_medias[i,2]<-mean(bio_data[,i]) #calculo la media y la escribo en su celda. #ojo, tengo que especificar con la coma que es un vector
}
head(mis_medias)

#while:primer ejemplo

i<-0
while (i<5){
  print(i)
  #i+1 si dejamos esto tal cual es un bucle infinito
  i<-i+1
}

#tirada de dados. Si sacamos 6 dejamos de tirar
dado<-0
while (dado!=6){
  #dado<-runif(n=1,min=1,max=6)#distribución aleatoria entre 1 y 6
  dado<-sample(1:6,size=1)#distribución aleatoria entre 1 y 6
  message(paste("Tu resultado es",dado),appendLF=TRUE)
  
}

PAQUETE foreach

#sapply, lapply y apply

apply(bio_data[1:19],2,FUN=mean)

packages<-c("rgbif","ALA4R","kableExtra","parallel","kableExtra","tidyverse")
sapply(packages,require,character.only=T)

lapply(packages,require,character.only=T)

for(i in especies$name){
  print(paste("Searching",i))
  occ_search(scientificName = i)%>%print()
}

sapply(especies$name,occ_search)

#condicionales if else ifelse

for (i in 1:nrow(especies)){
  if (especies[i,5] == TRUE){
    print (paste(especies[i,1],"es invasora"))
  } else {
    print(paste(especies[i,1],"es autóctona"))
  }
}


for (i in 1:nrow(especies)){
  if (especies[i,5] == TRUE){
    print (paste(especies[i,1],"es invasora"))
  } else {
    print(paste(especies[i,1],"es autóctona"))
  }
}

#introducir break y repeat
tiradas<-0
repeat{#repeat permite hacerlo hasta el infintio
  dado1<-sample(1:6,size=1)
  dado2<-sample(1:6,size=1)
  tiradas<-tiradas+1
  if(dado1==5 & dado1==dado2){
    message("Has sacado doble 5")
    print(tiradas)
    print(paste("Probabilidad",1-tiradas/36))
    break #nos saca del bucle.
  } else {
    message(paste0("dado1",dado1,"dado2", dado2))
    message(tiradas)
  }
}

4.2 Ejercicio guiado

Vamos a hacer un ejercicio entre tod@s para resolver una minería de datos de rgbif y hacer un análisis (muy simplificado) de la variabilidad climática de algunas especies. Vamos a ver primero un ejemplo:

mis_datos<-list()
for (i in 1:nrow(especies)){
  if (especies[i,5] == TRUE){
    print (paste(especies[i,1],"es invasora"))
    temp_occ<-occ_search(scientificName = toString(especies[i,1]))#el toString es una ñapa para factores
    mis_datos<-list(mis_datos,temp_occ)
  } else {
    print(paste(especies[i,1],"es autóctona"))
    temp_occ<-occ_search(scientificName = toString(especies[i,1]),country = "ES")
    mis_datos<-list(mis_datos,temp_occ)
    }TRUE
    
}

Si intentamos sacar de la lista mis_datos lo que hemos extraído de GBIF, nos encontramos con una estructura demasiado compleja. No nos sirve para responder a nuestra pregunta.

Buscar en bucle en gbif las coordenadas de las especies. - Si son autóctonas, visualizar sus coordenadas - Si son invasoras, contrastar con un boxplot los datos de dos variables bioclimáticas

for (i in 1:nrow(especies)){
  if (especies[i,5] == TRUE){
    print (paste(especies[i,1],"es invasora"))
    message("Searching occurrences")
    temp_occ<-occ_search(scientificName = toString(especies[i,1]),limit=2000)
    temp_occ<-data.frame(temp_occ$data)
    temp_occ<-cbind(temp_occ["decimalLongitude"],temp_occ["decimalLatitude"],temp_occ["country"])#invertimos el orden de la latitud y la longitud (3 y 4)
    biodata<-temp_occ[c(1,2)]
    message(paste("Extracting data from raster for species",especies[i,1]))
    biodata<-raster::extract(x=bioclim_world,y=biodata,df=TRUE)
    biodata<-cbind(biodata,temp_occ)
    biodata<-biodata[-c(1,24:26)]
    biodata_plot_boxplot12<-ggplot(biodata,aes(x=country,y=bio_12))+
      stat_boxplot(geom="errorbar")+geom_boxplot()+ggtitle(paste("Range of precipitation of",especies[i,1]))
    print(biodata_plot_boxplot12)
    biodata_plot_boxplot10<-ggplot(biodata,aes(x=country,y=bio_10))+
      stat_boxplot(geom="errorbar")+geom_boxplot()+ggtitle(paste("Average temperature of",especies[i,1]))
    print(biodata_plot_boxplot10)
  } else {
    print(paste(especies[i,1],"es autóctona"))
    temp_occ<-occ_search(scientificName = toString(especies[i,1]))
    temp_occ<-data.frame(temp_occ$data)
    temp_occ<-cbind(temp_occ["decimalLongitude"],temp_occ["decimalLatitude"],temp_occ["country"])#invertimos el orden de la latitud y la longitud (3 y 4)
    biodata<-temp_occ[c(1,2)]
    message(paste("Extracting data from raster for species",especies[i,1]))
    biodata<-raster::extract(bioclim_world,biodata,df=TRUE)
    biodata<-cbind(biodata,temp_occ)
    biodata<-biodata[-c(1,24:26)]
    biodata_plot<-ggplot(biodata,aes(x=decimalLongitude,y=decimalLatitude))+
      geom_point(aes(col=country))+ggtitle(especies[i,1])
    print(biodata_plot)
  }
  
}

Otro ejercicio: tirar 1000 veces dos dados e imprimir en pantalla la proporción de dobles 5 en relación a la teoría

contador<-0
  for (n in 1:1000){
    dice1<-sample(1:6,size=1)
    dice2<-sample(1:6,size=1)
    if(dice1==5 & dice1==dice2){
      contador=contador+1
      print("acierto")
    }
  }
  
  message(paste(contador/1000,1/36))

nsimulations=10000
  contador=0
  for (n in 1:nsimulations){
    dice1<-sample(1:6,size=1)
    dice2<-sample(1:6,size=1)
    if(dice1==5 & dice1==dice2){
      contador=contador+1
      print(paste("acierto, contador:",contador))
    }
  }

¿Y si se nos pide hacerlo con 10, 100, 10000….?

nsimulations=10000
  contador=0
  for (n in 1:nsimulations){
    dice1<-sample(1:6,size=1)
    dice2<-sample(1:6,size=1)
    if(dice1==5 & dice1==dice2){
      contador=contador+1
      print(paste("acierto, contador:",contador))
    }
  }

4.3 Creando nuestras propias funciones con R

¿Y si se nos pide hacer la misma operación pero también puede variar el número de caras del dado? Otra manera de optimizar el código es realizarlo a través de una función a la que podamos variar los parámetros que queramos.

#coge el siguiente código y extraer la función. ¿qué ha pasado? ¿qué hace la función?  
  contador=0
  for (n in 1:nsimulations){
    dice1<-sample(1:6,size=1)
    dice2<-sample(1:6,size=1)
    if(dice1==5 & dice1==dice2){
      contador=contador+1
      print(paste("acierto, contador:",contador))
    }
  }


#cuál es la diferencia con el código anterior?
  
  contador=0
  for (n in 1:nsimulations){
    dice1<-sample(1:caras,size=1)
    dice2<-sample(1:caras,size=1)
    if(dice1==5 & dice1==dice2){
      contador=contador+1
      print(paste("acierto, contador:",contador))
    }
  }
  
  dado(10,6)
  
#vamos a arreglar la función: valores por defecto, dar un retorno decente.  
  
  dado<-function(nsimulations=100, caras=6) {
    contador=0
    for (n in 1:nsimulations){
      dice1<-sample(1:caras,size=1)
      dice2<-sample(1:caras,size=1)
      if(dice1==5 & dice1==dice2){
        contador=contador+1
        #print(paste("acierto, contador:",contador)) #lo quitamos para que no moleste tanto print
      }
    }
    return(contador/nsimulations)
  }
  
  
dado()
mi_proporcion<-dado()
mi_proporcion  

simulaciones<-c(1,10,100,1000,10000,100000,100000)
sapply(simulaciones,dado,caras=6)

for (s in simulaciones){
  print(s)
  dado(s,caras=6)%>%print()
}

¿Qué manera de comparar simulaciones es más rápida? Lo podemos realizar con el paquete microbenchmark, para comparar lo que tardan las distintas expresiones.

install.packages("microbenchmark")
library(microbenchmark)

microbenchmark(for (s in simulaciones){
  print(s)
  dado(s,caras=6)%>%print()
}, sapply(simulaciones,dado,caras=6),times=10)

install.packages(“microbenchmark”) library(microbenchmark)

microbenchmark(for (s in simulaciones){ print(s) dado(s,caras=6)%>%print() }, sapply(simulaciones,dado,caras=6),times=10)

Vamos con otro ejemplo: hacer nuestra propia función de la media. ¿Es más eficiente?

#otra comparación: media() vs. mean()
media<-function(columna){
  suma<-sum(columna)
  promedio<-suma/length(columna)
  return(promedio)
}
media(c(4,6))
datos_random<-rnorm(100000,3,5)

#y ahora los comparamos: ocurre algo completamente inesperado para mí.
microbenchmark(media(datos_random),mean(datos_random),times=30)
microbenchmark(media(datos_random),mean(datos_random),times=30)%>%boxplot()


media(datos_random)
mean(datos_random)

5 Publicación de datos con R

5.1 Publicación con RMarkdown

Flujo de Markdown

my_data<-data.frame(matrix(ncol=1000,nrow=10000))
for(i in 1:1000){
 my_data[i]<-rnorm(10000,i,i/5)
 
  
}
resp<-my_data[1]
pred<-my_data[4:6]

pander(aov(X1~X2+X3+X5,data=my_data))

Analysis of Variance Model
	Df	Sum Sq	Mean Sq	F value	Pr(>F)
X2	1	0.02236	0.02236	0.5588	0.4548
X3	1	0.01245	0.01245	0.3111	0.577
X5	1	0.008269	0.008269	0.2067	0.6494
Residuals	9996	399.9	0.04001	NA	NA

5.2 Github

Acceso a Github

5.3 Shiny

Shiny

Curso R GBIF

Alejandro González

March 24, 2018