Presentación

En este espacio se encuentran los contenidos para el Módulo 1 - Introducción a R y estadística básica del Curso de Estadística ENZOEM, que tiene la siguiente estructura:

El lenguaje R y el entorno R-Studio: funciones y objetos. Materiales
Flujo de trabajo en R: directorio de trabajo, lectura y escritura de datos. Materiales
Visualización de datos en R y ggplot2. Materiales
Distribuciones de probabilidad y estadística básica. Presentaciones

1. El lenguaje R y el entorno R-Studio: funciones y objetos

¿Por qué R?

R es un entorno y lenguaje de programación con enfoque al análisis estadístico. Este tipo de herramientas pueden tener desventajas (estár más pendiente de aprender a utilizar la herramienta que a enterarse de la teoría, una curva de aprendizaje muy pronunciada, etc.). Sin embargo, también hay buenas razones para trabajar con esta herramienta:

R es una herramienta libre y gratuita para todos los sistemas operativos.
R no es solamente un software para análisis estadístico.
Es muy probable que tarde o temprano nos veamos “obligad@s” utilizar esta herramienta.

Descarga e instalación de R y RStudio

En este curso utilizaremos dos softwares que deberemos descargar e instalar:

R lo podemos descargar en https://cran.r-project.org/ para los sistemas operativos Windows, macOS y Linux.
RStudio lo podemos descargar en https://posit.co/download/rstudio-desktop/, también para cualquer sistema operativo

Una vez descargados, RStudio se vinculará con R automáticamente.

Un paseo por R y RStudio

Es importante diferenciar entre R y RStudio. El software o lenguaje que utilizamos para los análisis estadisticos es R, por tanto podemos utilizar la consola de R por sí sola. RStudio es un IDE (del inglés Integrated Development Environment), un entorno de desarrollo integrado, una “carcasa” que envuelve a R y facilita su uso. Es imporante familiarizarse con RStudio, ya que utilizaremos R a través de él. Tómate unos minutos para explorar sus diferentes paneles y personalizar su diseño:

Diferencia entre consola y script
Visualiza tu entorno de trabajo, tu ventana de gráficos (plots), tus librerías (paquetes) cargadas, etc.
También puedes cambiar el color de fondo, aumentar el tamaño de letra, etc.

Las `funciones()`

Las funciones() son la “maquinaria” de R, las que realizan el trabajo. Se pueden identificar porque generalmente van seguidas de unos paréntesis entre los cuales se colocan sus argumentos. Los argumentos de una función son elementos que necesita esa función para ejecutarse y suelen ir separados por comas $,$. Por ejemplo, existe una función que se llama “concatenar” (que en el lenguaje R se escribe c()), que simplemente sirve para unir en el mismo vector una serie de elementos (letras, números, etc.). Vamos a utilizar esa función para unir en un único vector una serie de números.

# La almohadilla se utiliza para incluir un comentario. Todo lo que se encuentre
# después de una almohadilla no se ejecutará en la consola.

# Unir los números 3, 7, 12 y 4 en un único vector
c(3, 7, 12, 4)

## [1]  3  7 12  4

¿Cuál es la función? ¿Cuáles son sus argumentos?

Imaginemos que ahora queremos hallar la media de esos números. Podemos utilizar otra función denominada mean() a la cual le introduciremos los numeros que hemos concatenado anteriormente como único argumento:

# Obtener la media de los números 3, 7, 12 y 4
mean(c(3, 7, 12, 4))

## [1] 6.5

Todas las funciones de R se alamcenan en “bibliotecas” o librerías (habitualmente denominados paquetes). Estos paquetes podemos entenderlos como cajas de herramientas donde se almacenan las herramientas que queremos utilizar. Hay algunos paquetes que vienen instalados y cargados por defecto en R. Sin embargo, otros los tenemos que descargar e instalar por primera vez y luego cargarlos en nuestra sesión cada vez que queramos utilizarlos. Por ejemplo, el paquete lme4 se utiliza frecuentemente para ajustar modelos generales lineales mixtos. Podemos descargarlo, instalarlo y cargarlo en nuestra sesión de R de la siguiente manera:

# Descargamos e instalamos el paquete
install.packages("lme4")
# Cargamos el paquete en nuestra sesión
library(lme4)

Una de las grandes ventajas (o inconvenientes?) de R es que es un software libre, por lo que cualquiera puede desarrollar sus propios paquetes con las herramientas (funciones) que necesite y ponerlo a disposición de la comunidad de usuarios. Si tenéis curiosidad, aquí podéis encontrar un pequeño tutorial sobre como hacerlo.

Los objetos

Los objetos en R son los contenedores donde almacenamos los resultados (outputs) de las funciones. Podemos identificarlos porque suelen aparecer por primera vez precediendo a los caracteres <-, que simbolizan una flecha que señala hacia la izquierda. Cada vez que se quiera crear un objeto se le ha de dar un nombre, el que queramos, aunque suele ser conveniente darle un nombre que tenga sentido. Por ejemplo, vamos a almacenar en un objeto que vamos a llamar “numeros” la concatenación de valores que creamos anteriormente:

# Almacenamos en un objeto llamado "numeros" el resultado de concatenar 3, 7, 12 y 4
numeros <- c(3, 7, 12, 4)

# Ahora podemos "llamar" a "numeros" para ver qué tiene dentro
numeros

## [1]  3  7 12  4

# También podemos utilizar a "numeros" dentro de otra función, por ejemplo mean()
mean(numeros)

## [1] 6.5

# Por último podemos guardar dentro de otro objeto el resultado de la linea anterior

m <- mean(numeros)
m

## [1] 6.5

Todos los objetos en R tienen una clase, que informa sobre el tipo de objeto que es. Por ejemplo, si es un vector de números será numeric, pero si lo que almacena son caracteres su clase será character. Hay muschísimas clases de objetos (¡incluso se pueden crear clases nuevas!). Conocer la clase de nuestros objetos es muy important, puesto que algunas funciones necesitan que sus argumentos sean de una clase específica, y sino no funcionarán. Por ejemplo, no podemos hacer la media de las letras “a”, “b” y “c”, pero sí podremos hacer la media de los números 1, 2 y 3.

# Para averiguar la clase de un objeto usamos la función class()
class(numeros)

## [1] "numeric"

# Vamos a crear un objeto con los caracteres "a", "b" y "c"
letras <- c("a", "b", "c")

# Exploramos la clase de letras
class(letras)

## [1] "character"

# Vamos a intentar hacer la media de letras
mean(letras)

## Warning in mean.default(letras): argument is not numeric or logical: returning
## NA

## [1] NA

2. Flujo de trabajo en R: directorio de trabajo, lectura y escritura de datos

Entorno y directorio de trabajo

Todos los objetos que vayamos creando o cargando en R se pueden visualizar en el panel Environment que suele situarse en la parte superior derecha en RStudio. En ese panel aparece información sobre los objetos, una previsualización o incluso, si el objeto lo permite, podemos hacer click en él y visualizarlos de forma intuitiva en formato “hoja de cálculo”.
También es importante ser conscientes del directorio de trabajo en el que estamos trabajando, esto es, la carpeta en la que se guardarán los ficheros que salgan de R, o la carpeta desde donde se cargarán los ficheros. Para saber nuestro directorio de trabajo actual utilizamos la función getwd(), mienstras que para cambiarla utilizaremos setwd("mi_directorio_de_trabajo_nuevo").

# Exploramos mi directorio de trabajo actual
getwd()

## [1] "/home/javifl/github/web/tutorials"

# Cambiamos el directorio de trabajo
setwd("/home/javifl/github/web/tutorials/estadisticaBasica_files")

Trabajo con `data.frame`

Uno de los objetos que más vamos a utilizar en la práctica es el data.frame, que podría asimilarse a un fichero similar a lo que podríamos visualizar en una hoja de cálculo tipo Excel. Vamos a utilizar la función read.csv leer un archivo de texto plano delimitado por comas (CSV del inglés “comma separated values”) y ensayar algunas cosas con él. Podéis echar un vistazo a este archivo pinchando aquí.

# Guardamos el archivo en un obejto que denominaremos "datos"
datos <- read.csv("https://raw.githubusercontent.com/jabiologo/web/master/tutorials/estadisticaBasica_files/datosPracticas.csv")

class(datos)

## [1] "data.frame"

# Podemos utilizar la función head() para explorar los primeros elementos de datos
head(datos)

##     peso    sex      hab      tri
## 1 419.35  Macho Agricola 47.87247
## 2 465.11 Hembra   Urbano 13.03228
## 3 434.17 Hembra   Urbano 18.05707
## 4 508.56  Macho Agricola 20.87707
## 5 514.75  Macho Agricola 15.29327
## 6 478.53  Macho   Urbano 15.52720

# Podemos utilizar la función str() para ver la estructura interna del data.frame
str(datos)

## 'data.frame':    300 obs. of  4 variables:
##  $ peso: num  419 465 434 509 515 ...
##  $ sex : chr  "Macho" "Hembra" "Hembra" "Macho" ...
##  $ hab : chr  "Agricola" "Urbano" "Urbano" "Agricola" ...
##  $ tri : num  47.9 13 18.1 20.9 15.3 ...

# O incluso podemos pedirle un resumen con summary()
summary(datos)

##       peso           sex                hab                 tri        
##  Min.   :330.0   Length:300         Length:300         Min.   :-4.733  
##  1st Qu.:430.9   Class :character   Class :character   1st Qu.:18.955  
##  Median :454.9   Mode  :character   Mode  :character   Median :25.853  
##  Mean   :453.0                                         Mean   :25.780  
##  3rd Qu.:475.9                                         3rd Qu.:32.421  
##  Max.   :555.1                                         Max.   :52.502

Tenemos varias formas de manejar los datos contenidos en un data.frame. Por ejemplo, si queremos explorar alguna columna podemos utilizar el símbolo del dólar $ para seleccionarla por su nombre. Por ejemplo datos$peso. Si lo que queremos es seleccionar un dato en concreto, podemos utilizar los corchetes [] para acceder a los elementos del data.frame indicando su posición en cuanto a filas y columnas separados por una coma con la forma datos[fila,columna]. Por ejemplo, si queremos seleccionar el dato que está en la fila 3 y columna 4, podemos emplear datos[3,4]. Si dejamos uno de los dos huecos en blanco, obtendremos toda la fila o toda la columna. Por ejemplo, si quremos todos los datos de la fila 7, podemos hacer datos[7, ]. Veamos estos ejemplos.

# Visualizamos los nombres de las columnas de datos
colnames(datos)

## [1] "peso" "sex"  "hab"  "tri"

# Visualizamos los 15 primeros elementos de la columna peso
datos$peso[1:15]

##  [1] 419.35 465.11 434.17 508.56 514.75 478.53 429.03 463.79 445.34 412.27
## [11] 443.78 329.97 414.74 461.40 388.19

# Visualizamos el dato que está en la fila 3 y columna 4
datos[3,4]

## [1] 18.05707

# Visualizamos todos los datos de la fila 7
datos[7,]

##     peso    sex    hab      tri
## 7 429.03 Hembra Urbano 32.48139

Por último, vamos a ver una forma más moderna (y cómoda) de trabajar con los datos a través de tydiverse. Tydiverse es una coleccion de paquetes de R con una sintaxis en teoría más intuitiva que se centran en el manejo de datos. Utilizan el formato de “tuberías” o pipes mediante el símbolo %>%. Veamos un ejemplo.

# Instalamos (en caso necesario) y cargamos las librerías de tidyverse
# install.packages(tidiverse)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     1.4.0
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

# Tomamos el objeto datos, y lo filtramos dejando tan solo los machos
datosMacho <- datos %>% filter(sex == "Macho")

# Visualizamos los primeros elementos del nuevo objeto
head(datosMacho)

##     peso   sex      hab      tri
## 1 419.35 Macho Agricola 47.87247
## 2 508.56 Macho Agricola 20.87707
## 3 514.75 Macho Agricola 15.29327
## 4 478.53 Macho   Urbano 15.52720
## 5 463.79 Macho Matorral 23.83045
## 6 412.27 Macho Agricola 46.89978

# También podemos obtener una muestra aleatoria de un número de filas
muestra11 <- datos %>% sample_n(11)

# Número de filas de muestra11
nrow(muestra11)

## [1] 11

muestra11

##      peso    sex      hab       tri
## 1  482.59  Macho Agricola 26.005255
## 2  461.66 Hembra   Urbano 13.256390
## 3  467.74 Hembra Matorral 21.298767
## 4  484.60  Macho Matorral 19.423977
## 5  414.65 Hembra Agricola 30.917451
## 6  478.27 Hembra   Urbano 16.110716
## 7  525.18  Macho Agricola 10.248180
## 8  478.92  Macho Agricola 33.059242
## 9  388.30 Hembra   Urbano 34.535212
## 10 465.71 Hembra   Urbano  9.599989
## 11 476.96 Hembra Matorral 14.317273

La ayuda `help()` y preguntar a Google (o ahora a ChatGPT con precaución)

Otra de las ventajas de R con respecto a otros lenguajes de programación es que sus paquetes deben de cumplir una serie de estándares para poder estar en el repositorio “oficial” CRAN. Uno de los requerimientos es que las funciones de los diferentes paquetes deben de estar bien documentadas, es decir, debe existir un manual que indique cómo se usa cada una de las funciones del paquete. Ese manual sigue siempre la misma estructura: nombre y descripción de la función, cómo se usa, los argumentos que necesita, el objeto que resulta al aplicar la función (Value), y unos ejemplos de cómo usar la función. Este manual o ayuda puede incluir más apartados, pero los mencionados suelen ser obligatorios. La forma de “llamar a la ayuda” en R es utilizando la función help(). Veamos un ejemplo:

# Vamos a intentar obtener la media del objeto letras
mean(letras)

## Warning in mean.default(letras): argument is not numeric or logical: returning
## NA

## [1] NA

# No entendemos por qué obtenemos este error... vamos a consultar la ayuda para
# entender qué necesita la función mean()
help(mean)

# Como vemos en el apartado Arguments, a la función mean() sólo se le pueden
# proporcionar objetos que sean de la clase numeric, logical, vectors, date,
# date-time y time_interval.
# Dado que la clase de letras es character, ahora entendemos por qué no podemos
# utilizar mean con letras.

Quizás la lección más importante es que normalmente ningún usuario habitual ha asistido a ningún curso especializado de R: la mayoría aprende a base de prueba/error y muchas horas delante del ordendaor con ERROR o WARNINGS en nuestra pantalla. Por eso es esencial no perder la paciencia y preguntar siempre a Google antes de preguntar al compañero que sabe (siemrpe hay alguien que ha tenido la misma duda antes que nosotros y ha dejado un comentario en un foro o blog). De esta forma interiorizaremos mucho mejor el funcionamiento de este lenguaje y aprenderemos mucho más rápido!

Merece la pena echarle un vistazo (o incluso imprimir si lo creemos conveniente) a esta cheatsheet sobre el uso básico de R.

Prácticas 1.1 y 1.2: Flujo de trabajo en R y RStudio

Los fungicidas triazoles son compuestos químicos que se aplican habitualmente en semillas de cultivos para prevenir el crecimiento de hongos patógenos de plantas. Sin embargo, cuando las semillas son consumidas por la fauna silvestre, estos compuestos pueden producir efectos crónicos perjudiciales en su sauld y desarrollo. Queremos estudiar el efecto de los fungicidas triazoles sobre la condición corporal (peso) en perdices rojas Alectoris rufa. Para ello se han capturado un total de 300 perdices en tres hábitats diferentes (semiurbano, agrícola y monte matorralizado) a las que se les ha sexado y extraido muestras de heces para obtener la concentración de fungicidas triazoles (ngramos de compuesto/gramo de heces).
En el los materiales de este modulo se puede encontrar el archivo enzoem1_1-2_datos.xlsx. Este archivo pretende reflejar los típicos problemas a los que nos solemos enfrentar cuando empezamos a analizar nuestros propios datos con R. El objetivo de esta práctica es sencillo: realizar un histograma de los pesos de los animales capturados con la función hist()y un gráfico de cajas que compare el peso entre machos y hembras con la función boxplot().

3. Visualización de datos en R y `ggplot2`

Conceptos básicos de ggplot2

ggplot2 es un paquete de R para crear gráficos de alta calidad. Se basa en una gramática que permite describir y construir gráficos, incorporando capas de información complementarias de una forma coherente.

La gramática de los gráficos se basa en la idea de que cualquier gráfico se puede construir a partir de los siguientes componentes:

Datos (Data): El conjunto de datos que se va a visualizar.
Mapeo: El mapeo define cómo se asignan las partes de los datos a los atributos estéticos de los objetos geométricos. Esto se hace usando la función aes()
Capas (Layers): Las capas son el corazón de cualquier gráfico. Cada capa consiste en:
- Geometría (Geometries): Determina cómo se muestran los datos (puntos, líneas, rectángulos, etc.).
- Transformación estadística (Statistics): Puede calcular nuevas variables a partir de los datos.
- Ajuste de posición: Determina dónde se muestra cada pieza de datos. Las capas se construyen usando funciones geom_*() y stat_*().
Facetas (Facets): Permiten dividir los datos en múltiples paneles basados en una o más variables. Esto se hace usando funciones como facet_wrap() y facet_grid().
Coordenadas (Coordinates): Definen el sistema de coordenadas del gráfico. Por ejemplo, coord_flip() invierte los ejes x e y.
Temas Visuales (Themes): Personalización de la apariencia del gráfico. theme( )

A partir de estos componentes es posible generar una gran variedad de tipos y formatos de gráficos de una forma flexible.

CheatSheet (chuleta código)

Como crear un gráfico básico

Empecemos creando un gráfico de dispersión básico a partir de los datos palmerpenguins. Los datos fueron recogidos y puestos a disposición por la Dra. Kristen Gorman y la Estación LTER Palmer, en la Antártida.

Artwork by @allison_horst

Instalamos y/o cargamos los datos
```
library (palmerpenguins)
```
```
## 
## Adjuntando el paquete: 'palmerpenguins'
```
```
## The following objects are masked from 'package:datasets':
## 
##     penguins, penguins_raw
```
```
data("penguins")
```
Valoramos la calidad de los datos de forma gráfica, identificando la distribución de los NA en la matriz de datos.
```
#install.packages("visdat")
library ("visdat")
visdat::vis_dat(penguins)
```
Listado y descripción de las variables pertenecientes a la BD de PalmerPenguins
- species: Especie de pingüino: Adélie (Adelia), Chinstrap (Barbijo), Gentoo (Papúa).
- island: Isla donde se encontró el pingüino (Biscoe, Dream, Torgersen).
- bill_length_mm: Longitud del pico en milímetros.
- bill_depth_mm: Profundidad del pico en milímetros.
- flipper_length_mm: Longitud de la aleta en milímetros.
- body_mass_g: Masa corporal en gramos.
- sex: Sexo del pingüino (macho, hembra).
  
  Para simplificar los análisis y reducir los potenciales errores o warnings se eliminan los NA

library(tidyr)
penguins1 <- penguins %>% drop_na()
names (penguins1)

## [1] "species"           "island"            "bill_length_mm"   
## [4] "bill_depth_mm"     "flipper_length_mm" "body_mass_g"      
## [7] "sex"               "year"

visdat::vis_dat(penguins1)

Generamos la gráfica base

Representamos la longitud del pico frente a la anchura del pico

-   "ggplot (data = penguins1, aes(x = bill_length_mm, y = bill_depth_mm))": Especifica el conjunto de datos y el mapeo estético.

-   "geom_point()": Añade puntos al gráfico.

ggplot2 permite guardar esta figura como un objeto en una variable. Luego puedes usar este objeto añadiéndole nuevas capas de información, por ejemplo uniendo los puntos con una línea.

p1 + geom_line(color = "firebrick")

Personalización del Gráfico

Se puede personalizar el gráfico añadiendo más capas y modificando la estética.

Añadir Título y Etiquetas de los ejes:

ggplot(data = penguins1, 
       aes(x = bill_length_mm, y = bill_depth_mm))+
 geom_point() +
 labs(title = "Relación entre la longitud y anchura del pico",
       x = "Longitud del pico (mm)",
       y = "Anchura del pico (mm)")

Cambiar Colores y Tamaños:

Modificamos la estética, incorporando un color para cada isla, además el tamaño del punto estará condicionado por el peso del animal. También podríamos usar una forma diferente dependiendo de la especies.
```
ggplot(data = penguins1,
       aes(x = bill_length_mm, y = bill_depth_mm, color = island, size = body_mass_g, shape = species ))+
  geom_point() +
  labs(title="Relación entre la longitud y anchura del pico",
       x = "Longitud del pico (mm)",
       y = "Anchura del pico (mm)")
```
- “color = island”: Colorea los puntos según la isla de origen.
- “size” = body_mass_g: Ajusta el tamaño de los puntos según el peso de los animales.
- “shape” = permite ajustar la forma de los puntos en base a una variable

Facetado

El facetado permite crear múltiples gráficos basados en subgrupos de datos. Por ejemplo representando por separado la relación para cada especie.

ggplot(data = penguins1,
       aes(x = bill_length_mm, y = bill_depth_mm))+
 geom_point() +
 geom_line(color = "blue", linewidth = 0.2)+ 
 facet_wrap(~ species) +
 labs(title="Relación entre la longitud y anchura del pico",
       x = "Longitud del pico (mm)",
       y = "Anchura del pico (mm)")

Temas

Los temas permiten personalizar la apariencia del gráfico. En esta gráfica se han incorporado diferentes aproximaciones (lm y loess) para valorar visualmente el ajuste de la asociación para cada especie.

ggplot(data = penguins1,
       aes(x = bill_length_mm, y = bill_depth_mm))+
 geom_point(color = "#8B8878", size = 2) +
 theme_minimal() +
 facet_wrap(~ species) +
 geom_smooth(aes(colour = "loess"), method = "loess", se = F) +  geom_smooth(aes(colour = "lm"), method = "lm", se = FALSE) +  
 labs(title="Relación entre la longitud y anchura del pico",
       x = "Longitud del pico (mm)",
       y = "Anchura del pico (mm)")

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Cambio del color de los puntos en base a una variable (los colores se generan automáticamente). Usando diferentes colores, se muestran los individuos capturados en cada isla. También se modifica el color del texto de los ejes

ggplot(data = penguins1,
       aes(x= bill_length_mm, y= bill_depth_mm, color=island))+
 geom_point(size = 2) +
 theme_minimal() +
 facet_wrap(~ species) +
 geom_smooth(aes(colour = "loess"), method = "loess", se = F) +  geom_smooth(aes(colour = "lm"), method = "lm", se = FALSE) +  
 labs(title="Relación entre la longitud y anchura del pico",
       x = "Longitud del pico (mm)",
       y = "Anchura del pico (mm)")+
 theme(axis.text = element_text(color = "blue", size = 12),
        axis.text.x = element_text(face = "italic"))

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Se puede seleccionar manualmente los diferentes colores a utilizar para los puntos y las 2 líneas de tendencia generadas en el gráfico.

ggplot(data = penguins1,
       aes(x= bill_length_mm, y= bill_depth_mm, color=island))+
 geom_point(size = 2) +
 scale_color_manual(values = c("darkorange", "darkorchid", "cyan4","red","black"))+
 theme_minimal() +
 facet_wrap(~ species) +
 geom_smooth(aes(colour = "loess"), method = "loess", se = F) +  geom_smooth(aes(colour = "lm"), method = "lm", se = FALSE) +  
 labs(title="Relación entre la longitud y anchura del pico",
       x = "Longitud del pico (mm)",
       y = "Anchura del pico (mm)")+
 theme(axis.text = element_text(color = "blue", size = 12),
        axis.text.x = element_text(face = "italic"))

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Representación de la longitud de las aletas frente al peso de los animales diferenciando las especies por color y forma. También se activa el intervalo de confianza alrededor de la línea de suavizado.

ggplot(data=penguins1,aes(x= flipper_length_mm, y= body_mass_g))+
  geom_point(aes(color = species, 
                 shape = species,
                 size = 1))+
  geom_smooth(aes(group = species,color = ""), method = "lm", se = T) +
  scale_color_manual(values = c("orange","purple","cyan4", "black"))+
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

Colores

Hay diferentes opciones para poder incorporar diferentes colores en glos gráficos

el uso de los colores estándares que tiene R. (leer mas)

Se puede consultar el nombre y la codificación HEX de los colores en Color-Chart eg: “#8B8878” =cornskill4
Alternativamente se pueden utilizar paletas específicas Paletas de colores

Uso de una librería para ajustar una paleta de colores prefijada

library(RColorBrewer)

ggplot(data = penguins1, aes(x = bill_length_mm, y = bill_depth_mm, color = species)) +
 geom_point(size=3) +
 scale_color_brewer(palette = "Set2") +
 theme_minimal()

Explorando otras geometrías

Necesitaremos crear diferentes geometrías (formas) a partir de nuestros datos (por ejemplo, diagramas de barras, histogramas, diagramas de dispersión, diagramas de caja).

Esto se hace añadiendo capas “geoms” al comando inicial ggplot(). Hay muchas funciones ggplot. Cada una de estas funciones empieza por “geom_”, por lo que nos referiremos a ellas genéricamente como geom_XXXX(). Algunas de los geoms más comunes son:

Puntos (e.g. scatter plots) - geom_point()
Gráficos de líneas - geom_line() or geom_path()
Líneas de tendencias - geom_smooth()
Histogramas - geom_histogram()
Gráficos de barras - geom_bar() or geom_col()
Gráficos de caja - geom_boxplot()

Histograma

Histograma representando la longitud del pico por especies, se le otorga cierta transparencia (alpha = 0.4), los colores se establecen de forma manual (scale_fill_manual)

ggplot(data = penguins1, mapping = aes(x=bill_length_mm))+
  geom_histogram (aes(fill = species), alpha = 0.4, position = "identity") + 
  scale_fill_manual(values = c("orange","orchid","cyan2"))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Variaciones del argumento position:

“stack”: Apila los elementos uno sobre otro. Es útil para gráficos de barras apiladas.
“dodge”: Coloca los elementos uno al lado del otro. Es útil para comparar diferentes grupos en gráficos de barras.
“fill”: Similar a stack, pero ajusta las alturas de las barras para que sumen 1. Es útil para mostrar proporciones.
“jitter”: Añade un pequeño desplazamiento aleatorio a los puntos para evitar la superposición. Es útil para gráficos de dispersión con muchos puntos.

Histográma de densidades para las tres especies

ggplot(data = penguins1, mapping = aes(x = bill_length_mm)) +
 geom_density(aes(fill = species), alpha = 0.4, position = "identity") +
 scale_fill_manual(values = c("orange", "orchid", "cyan2")) +
 labs(title = "Proportional density")

Gráfico de barra

ggplot(data = penguins1, aes(x = factor(species))) +
 geom_bar(fill = "skyblue", color = "black") +
 labs(title = "Número de pinguinos por espcie ",
      x = "Especies de pinguinos",
      y = "Frecuencia") +
 theme_minimal()

Contando número de individuos por especie e isla

ggplot(penguins1, aes(x = island, fill = species)) +
  geom_bar(alpha = 0.8) +
  scale_fill_manual(values = c("orange","orchid","cyan4"), 
                    guide = "none") +
  theme_minimal() +
  facet_wrap(~species, ncol = 1) +
  coord_flip()

Contando número de individuos por especie e isla, y diferenciando por sexo

ggplot(penguins1, aes(x = island, fill = interaction(species, sex))) +
 geom_bar(alpha = 0.8) +
 scale_fill_manual(values = c("orange", "orchid", "cyan4", "lightblue", "plum", "darkgreen"), 
         guide = guide_legend(title = "Species and Sex")) +
 theme_minimal() +
 facet_wrap(~species, ncol = 1) +
 coord_flip ()

Gráfico de caja

ggplot(data = penguins1, aes(x = species, y= body_mass_g)) +
 geom_boxplot(fill = "skyblue", color = "black") +
 labs(title = "Peso de los animales por espcie ",
      x = "Especies de pinguinos",
      y = "Peso") +
 theme_minimal()

pb<-ggplot(data = penguins1, aes(x = species, y= body_mass_g)) +
 geom_boxplot(fill = "indianred", color = "black") +
 labs(title = "Peso de los animales por espcie ",
      x = "Especies de pinguinos",
      y = "Peso") +
 theme_minimal()

 pb+ coord_flip()

Combinando diferentes gráficos (Grafíco de Violín y distribución de los puntos)

pb+ coord_flip()+ 
    geom_violin(alpha=0.5)+ 
    geom_jitter(aes(color=species),
              alpha=0.5, 
              position = position_jitter(width = .2))

4. Distribuciones de probabilidad y estadística básica

Tipos de variables

Cuantitativas
- Continuas
- Discretas
Semicuantitativas u ordinales
Cualitativas o nominales

Distribuciones de probabilidad: la Normal o Gaussiana

Las distribuciones de probabilidad son funciones matemáticas que definen el comportamiento de una variable. Pueden entenderse como modelos teóricos que podemos utilizar para describir las mediciones que tomamos durante nuestros experimientos. Las distribuciones se definen mediante un nombre y unos parámetros, que manejan y describen su comportamiento.

Una de las distribuciones de probabilidad más conocida es la distribución Normal o Gaussianas. Esta distribución describe variables continuas que pueden ir desde menos infinito hasta más infinito (rango). Tiene dos parámetros: la media ($\mu$) y la varianza ($\sigma$). Variando los valores que toman esos dos parámetros podemos obtener diferentes formas para esa distribución.

¿Se te ocurre alguna variable que no pueda ser descrita por una normal?
¿Conoces el nombre de alguna otra distribución de probabilidad?

Para saber si una determinada variable sigue una distribución normal podemos utilizar varios test estadísticos. Entre los más utilizados se encuentra el test de Shapiro-Wilk. Lo que hace este test es comparar nuestra variable con una distribución normal teórica, asumiendo que ambas son iguales (nuestra hipótesis nula $H_0$). Si rechazamos la hipótesis nula ($p-valor < 0.05$), querrá decir que las dos distribuciones difieren, por lo que no podremos decir que nuestra variable se comporta como una distribución normal. Veamos un ejemplo:

# Guardamos el archivo en un obejto que denominaremos "datos"
datos <- read.csv("https://raw.githubusercontent.com/jabiologo/web/master/tutorials/estadisticaBasica_files/datosPracticas.csv")

# Podemos utilizar la función str() para ver la estructura interna del data.frame
str(datos)

## 'data.frame':    300 obs. of  4 variables:
##  $ peso: num  419 465 434 509 515 ...
##  $ sex : chr  "Macho" "Hembra" "Hembra" "Macho" ...
##  $ hab : chr  "Agricola" "Urbano" "Urbano" "Agricola" ...
##  $ tri : num  47.9 13 18.1 20.9 15.3 ...

# Primeros elementos de la base de datos
head(datos)

##     peso    sex      hab      tri
## 1 419.35  Macho Agricola 47.87247
## 2 465.11 Hembra   Urbano 13.03228
## 3 434.17 Hembra   Urbano 18.05707
## 4 508.56  Macho Agricola 20.87707
## 5 514.75  Macho Agricola 15.29327
## 6 478.53  Macho   Urbano 15.52720

# Histograma de nuestra variable dependiente (peso)
hist(datos$peso)

# Comprobamos la normalidad de nuestra variable
shapiro.test(datos$peso)

## 
##  Shapiro-Wilk normality test
## 
## data:  datos$peso
## W = 0.9966, p-value = 0.7764

# ¿Qué significa aquí un p-valor > 0.05?
# ¿Cuáles son aquí nuestras hipótesis nula y alternativa?
# Cómo trabaja la función shapiro.test()?
help("shapiro.test")

Comparación entre medias

Uno de los análisis estadísticos más comunes es la comparación de medias entre diferentes grupos: diferencias de talla entre machos y hembras, diferencias en la respuesta de un marcador molecular entre tratamiento y control, etc. Cuando sólo tenemos dos grupos podemos utilizar un T-test (o su “versión” no paramétrica U de Mann-Whitney). Veamos un ejemplo.

# Diferencia entre medias
# ¿Quién pesa más, los machos o las hembras? Filtramos los datos usando tidyverse
datosHe <- datos %>% filter(sex == "Hembra")
datosMa <- datos %>% filter(sex == "Macho")

# Recordamos: 
# H0 (hipótesis nula): no hay diferencias significativas entre gurpos
# H1 (hipótesis alternativa): hay diferencias significativas entre gurpos

# Paramétrico: T-Student
t.test(datosHe$peso, y = datosMa$peso)

## 
##  Welch Two Sample t-test
## 
## data:  datosHe$peso and datosMa$peso
## t = -7.5131, df = 294.22, p-value = 6.974e-13
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -35.46227 -20.74003
## sample estimates:
## mean of x mean of y 
##  438.7512  466.8524

# No paramétrico: U de Mann-Whitney
wilcox.test(datosHe$peso, y = datosMa$peso)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  datosHe$peso and datosMa$peso
## W = 6161, p-value = 1.276e-11
## alternative hypothesis: true location shift is not equal to 0

En este caso vemos que sí que existen diferencias significativas entre el peso de los machos y de las hembras, puesto que en el T-test obtenemos un P-valor muy bajo, 6.974e-13. Además, el estdístico t nos da algo de información adicional… ¿Quién pesa más, los machos o las hembras? ¿Qué le ocurre al estadístico t si cambiamos el orden de esta forma: t.test(datosMa$peso, y = datosHe$peso)?

¿Cuál es mi variable dependiente o respuesta? ¿Cuál es mi variable independiente o predictor? ¿De qué tipo es cada una de nuestras variables?
¿Podrias realizar un gráfico para representar estas diferencias en el peso?

Comparación entre dos variables cuantitativas continuas

Otro de los análisis estadísticos más comunes es la comparación de dos variables cuantitativas continuas: peso y talla, cantidad de un marcador molecular con un índice de condición física, etc. Para realizar esta comparación podemos utiilizar las correlaciones. Veamos un ejemplo.

# ¿Cómo afecta la concentración de fungidida al peso?
# Correlación entre dos variables continuas

# Paramétrico: Correlación de Pearson
cor.test(datos$tri, datos$peso, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  datos$tri and datos$peso
## t = -23.23, df = 298, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8395769 -0.7583332
## sample estimates:
##        cor 
## -0.8026473

# No paramétrico: Correlación de Spearman
cor.test(datos$tri, datos$peso, method = "spearman")

## Warning in cor.test.default(datos$tri, datos$peso, method = "spearman"): Cannot
## compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  datos$tri and datos$peso
## S = 8050549, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.7890307

# Regresión
regresion <- lm(peso ~ tri, data = datos)
summary(regresion)

## 
## Call:
## lm(formula = peso ~ tri, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -53.822 -15.102  -0.568  15.030  60.303 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 526.1267     3.3748  155.90   <2e-16 ***
## tri          -2.8370     0.1221  -23.23   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.05 on 298 degrees of freedom
## Multiple R-squared:  0.6442, Adjusted R-squared:  0.643 
## F-statistic: 539.6 on 1 and 298 DF,  p-value: < 2.2e-16

# Gráfico de dispersión (scatter-plot)
plot(datos$tri, datos$peso)

¿Qué significa que el coeficiente de correlación de Pearson sea negativo (-0.8026473)? ¿Qué relación existe entre este coeficiente y el Multiple R-squared obtenido a en la regresión?
Es importante recordar que esta regresión nos permite construir la fórmula de la recta con la cual podríamos teóricamente predecir los valores de nuestra variable respuesta en función de nuestra variale predictora. De forma general, esta fórmula se podría escribir de la siguiente manera:

\[\begin{equation} \text{Variable respuesta = Intercepto + Coeficiente(Variable_predictora)} \end{equation}\]

Al coeficiente también lo podemos llamar “estimador” y a todos los elementos de esa fórmula se suelen sustitur por diferetes caracteres, aunque suelen significar lo mismo. Aquí tienes algunos ejemplos: \[\begin{equation} Y = a + bX \end{equation}\] \[\begin{equation} Y = \alpha + \beta X \end{equation}\] \[\begin{equation} Y = \beta_0 + \beta_1 X \end{equation}\]

Así, teniendo en cuenta los resultados de summary(regresion), podríamos escribir la siguiente fórmula: \[\begin{equation} \text{Peso} = 526.1267 - 2.837 * (\text{concentración de triazoles}) \end{equation}\]

De esta forma, podríamos predecir el peso de una hipotética perdiz a partir de una concentración de triazol dada.

El concepto de residuo

En el contexto de los modelos lineales, el concepto de residuo hace referencia a “aquella información que no somos capaces de explicar en función de nuestros predictores”. Es la distancia que existe entre el valor observado y el valor predicho. Para poder utilizar los test paramétricos, estos residuos deberían seguir una distribución normal.

summary(regresion)

## 
## Call:
## lm(formula = peso ~ tri, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -53.822 -15.102  -0.568  15.030  60.303 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 526.1267     3.3748  155.90   <2e-16 ***
## tri          -2.8370     0.1221  -23.23   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.05 on 298 degrees of freedom
## Multiple R-squared:  0.6442, Adjusted R-squared:  0.643 
## F-statistic: 539.6 on 1 and 298 DF,  p-value: < 2.2e-16

# Calculamos los residuos
datos$residuos <- residuals(regresion)

# Histograma de los residuos
hist(datos$residuos)

# Comprobamos la normalidad de los residuos
shapiro.test(datos$residuos)

## 
##  Shapiro-Wilk normality test
## 
## data:  datos$residuos
## W = 0.99611, p-value = 0.6711

A continuación realizaremos una comparación entre las medias del peso de de todos los individuos entre los diferentes hábitats. ¿Hay diferencias significativas? ¿Entre qué hábitats?
¿Se te ocurre cómo podríamos analizar a la vez el efecto de los fungicidas, el sexo y el hábitat sobre el peso de las perdidces?

ENZOEM Módulo 1 - Introducción a R y estadística básica