Código
# instalar/cargar paquetes
::load_packages(
sketchyc("ggplot2",
"viridis"
) )
12 de diciembre de 2024
Aprender las principales herramientas de simulación de datos en R
Comprender la utlidad de usar datos simulados para entender el comportamiento de las herramientas estadísticas
Paquetes a utilizar en este manual:
La estadística nos permite inferir patrones en los datos. Solemos utilizar conjuntos de datos reales para enseñar estadística. Sin embargo, puede ser circular entender el funcionamiento interno de una herramienta estadística probando su capacidad para inferir un patrón que no estamos seguros de encontrar en los datos (y no tenemos idea del mecanismo que produjo ese patrón). Las simulaciones nos permiten crear escenarios controlados en los que conocemos con seguridad los patrones presentes en los datos y los procesos subyacentes que los han generado.
R ofrece algunas funciones básicas para la simulación de datos. Las más utilizadas son las funciones generadoras de números aleatorios. Los nombres de estas funciones comienzan con r (r____()
). Por ejemplo, runif()
:
El resultado es un vector numérico de longitud 100 (n = 100
):
[1] 9.889093 3.977455 1.156978 0.697487 2.437494 7.920104 3.400624 9.720625
[9] 1.658555 4.591037 1.717481 2.314771 7.728119 0.963015 4.534478 0.847007
[17] 5.606659 0.087046 9.857371 3.165848 6.394489 2.952232 9.967037 9.060213
[25] 9.887391 0.656457 6.270388 4.904750 9.710244 3.622208 6.799935 2.637199
[33] 1.857143 1.851432 3.792967 8.470244 4.980761 7.905856 8.384639 4.569039
[41] 7.994758 3.819431 7.597012 4.367756 9.042177 3.195349 0.825691 8.162891
[49] 8.984762 9.664964 5.730689 7.200795 7.740586 6.277608 7.229893 3.868313
[57] 1.627908 1.872283 3.912495 2.739012 1.919177 5.043918 7.638404 6.936689
[65] 5.440542 6.590872 4.687284 4.818055 3.370636 4.245263 2.870151 6.011915
[73] 8.407423 6.208370 1.345516 5.677224 4.434263 4.379754 6.236172 9.326533
[81] 8.884926 8.785406 2.421769 7.414538 3.876563 0.789517 0.948356 7.621427
[89] 3.478940 4.167667 3.440162 0.084109 9.115750 1.822054 7.228034 5.719633
[97] 5.400364 3.549474 8.240918 1.861368
Podemos explorar el resultado graficando un histograma:
Muestra una distribución uniforme que va de 0 a 10.
También podemos simular números aleatorios procedentes de una distribución normal utilizando rnorm()
:
Tenga en cuenta que todas las funciones generadoras de números aleatorios tienen el argumento ‘n’, que determina la longitud del vector generado (es decir, el número de números aleatorios), además de algunos argumentos adicionales relacionados con parámetros específicos de la distribución.
Las variables continuas (es decir, los vectores numéricos) pueden convertirse en variables discretas (es decir, números enteros) simplemente redondeándolas:
[1] 15.4269 8.5450 15.1387 8.9203 9.3693
[1] 15 9 15 9 9
Ejercicio 1
¿Qué hacen las funciones rbinom()
y rexp()
?
Ejecútela y haga histogramas de sus resultados
¿Qué hacen los argumentos ‘mean’ y ‘sd’ en rnorm()
? Juegue con diferentes valores y comprueba el histograma para hacerse una idea de su efecto en la simulación
La forma más sencilla de generar variables categóricas es utilizar el vector de ejemplo letters' (o
LETTERS’) para asignar niveles de categoría. Podemos hacerlo utilizando la función rep()
. Por ejemplo, el siguiente código crea un vector categórico (caracteres) con dos niveles, cada uno con 4 observaciones:
También podemos replicar este patrón utilizando el argumento ‘times’. Este código replica el vector anterior 2 veces:
[1] "a" "a" "a" "a" "b" "b" "b" "b" "a" "a" "a" "a" "b" "b" "b" "b"
Otra opción es simular una variable a partir de una distribución binomial y luego convertirla en un factor:
[1] 1 0 0 1 1 0 0 0 0 0 0 0 1 1 1 0 0 1 0 1 0 0 0 0 0 0 0 1 0 1 1 1 0 0 0 1 0 0
[39] 1 1 0 0 0 0 0 1 0 1 1 1
[1] b a a b b a a a a a a a b b b a a b a b a a a a a a a b a b b b a a a b a a
[39] b b a a a a a b a b b b
Levels: a b
La otra herramienta importante de R para jugar con datos simulados es sample()
. Esta función permite tomar muestras de tamaños específicos de vectores. Por ejemplo, tomemos el ejemplo del vector ‘letters’:
[1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p" "q" "r" "s"
[20] "t" "u" "v" "w" "x" "y" "z"
Podemos tomar una muestra de este vector como es:
El argumento ‘size’ nos permite determinar el tamaño de la muestra. Tenga en cuenta que obtendremos un error si el tamaño es mayor que el propio vector:
Error in sample.int(length(x), size, replace, prob): cannot take a sample larger than the population when 'replace = FALSE'
Esto sólo puede hacerse cuando el muestreo es con reemplazo (replacement). El muestreo con reemplazo puede aplicarse estableciendo el argumento replace = TRUE
:
[1] "j" "i" "h" "i" "k" "u" "n" "i" "w" "k" "c" "t" "m" "p" "i" "h" "d" "w" "c"
[20] "q" "j" "e" "e" "j" "k" "c" "l" "t" "r" "h"
A menudo, las simulaciones deben repetirse varias veces para descartar resultados espurios debidos al azar o simplemente para probar diferentes parámetros. Las funciones de simulación de datos mencionadas anteriormente pueden ejecutarse varias veces (por ejemplo, iteradas) utilizando la función replicate()
:
[1] "list"
[[1]]
[1] 1.19184 -0.33992
[[2]]
[1] 0.78911 -0.63213
[[3]]
[1] -1.49312 -0.13441
El último truco que necesitamos para ejecutar simulaciones en R es la capacidad de reproducir una simulación (es decir, obtener exactamente los mismos datos y resultados simulados). Esto puede ser útil para que otros investigadores puedan ejecutar nuestros análisis exactamente de la misma manera. Esto puede hacerse fácilmente con la función set.seed()
. Pruebe a ejecutar el siguiente código. Debería obtener la misma salida:
Ahora que sabemos cómo simular variables continuas y categóricas. Podemos juntarlas para crear conjuntos de datos simulados. Esto se puede hacer utilizando la función data.frame()
:
grupo | tamano |
---|---|
a | 4.8158 |
a | 3.6287 |
a | 4.4008 |
b | 5.2946 |
b | 5.3898 |
b | 3.7919 |
Por supuesto, podríamos añadir más variables a este cuadro de datos:
grupo | individuo | tamano | peso |
---|---|---|---|
a | A | 4.6363 | 109.874 |
a | B | 3.3733 | 107.414 |
a | C | 4.7435 | 100.893 |
b | D | 6.1018 | 90.451 |
b | E | 5.7558 | 98.049 |
b | F | 4.7618 | 109.255 |
Y eso es un juego de datos simulados en su forma más básica. Se parece mucho al tipo de datos con los que trabajamos en biología.
El Teorema del Límite Central afirma que, si tomamos muestras aleatorias de una población, los promedios de esas muestras seguirán una distribución normal, aunque la población no esté distribuida normalmente. Además, la distribución normal resultante debe tener un promedio cercano al promedio de la población. El teorema es un concepto clave para la estadística inferencial, ya que implica que los métodos estadísticos que funcionan para las distribuciones normales pueden ser aplicables a muchos problemas que implican otros tipos de distribuciones. No obstante, el objetivo aquí es sólo mostrar cómo se pueden utilizar las simulaciones para entender el comportamiento de los métodos estadísticos.
Para comprobar si esas afirmaciones básicas sobre el Teorema del Límite Central son ciertas, podemos utilizar datos simulados en R. Vamos a simular una población de 1000 observaciones con una distribución uniforme:
Podemos tomar muestras aleatorias usando sample()
así:
[1] 9.28420 1.02626 2.57517 3.32485 6.89990 2.29404 0.33737 8.21366 3.30364
[10] 8.03793 2.59174 7.81770 5.65426 0.63831 2.83470 4.20434 4.76330 4.42193
[19] 6.97830 7.92625 0.68121 3.52323 6.51103 5.38289 7.97210 1.80062 4.21282
[28] 3.33866 8.91223 4.71163
Este proceso puede ser replicado varias veces con replicate()
:
El código anterior toma 100 muestras con 30 valores cada una. Ahora podemos comprobar la distribución de las muestras:
… asi como el promedio:
Como era de esperar, las muestras siguen una distribución normal con una media cercana a la media de la población, que es:
Probemos con una distribución más compleja. Por ejemplo, una distribución bimodal:
Ejercicio 2
Intenta explorar el Teorema del Límite Central como en el caso anterior, pero esta vez utilizando:
rexp()
)rlnorm()
)
R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: America/Costa_Rica
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] viridis_0.6.5 viridisLite_0.4.2 ggplot2_3.5.1 knitr_1.48
loaded via a namespace (and not attached):
[1] gtable_0.3.5 jsonlite_1.8.9 dplyr_1.1.4
[4] compiler_4.3.2 crayon_1.5.3 tidyselect_1.2.1
[7] stringr_1.5.1 gridExtra_2.3 scales_1.3.0
[10] yaml_2.3.10 fastmap_1.2.0 R6_2.5.1
[13] labeling_0.4.3 generics_0.1.3 htmlwidgets_1.6.4
[16] tibble_3.2.1 munsell_0.5.1 xaringanExtra_0.8.0
[19] pillar_1.9.0 rlang_1.1.4 utf8_1.2.4
[22] stringi_1.8.4 xfun_0.47 cli_3.6.3
[25] withr_3.0.1 magrittr_2.0.3 digest_0.6.37
[28] grid_4.3.2 rstudioapi_0.16.0 remotes_2.5.0
[31] packrat_0.9.2 lifecycle_1.0.4 vctrs_0.6.5
[34] evaluate_1.0.0 glue_1.8.0 farver_2.1.2
[37] fansi_1.0.6 colorspace_2.1-1 sketchy_1.0.3
[40] rmarkdown_2.28 tools_4.3.2 pkgconfig_2.0.3
[43] htmltools_0.5.8.1
---
title: <font size="7"><b>Simulación de datos</b></font>
---
```{r,echo=FALSE,message=FALSE}
options("digits"=5)
options("digits.secs"=3)
```
```{r, echo = FALSE, message=FALSE}
library(knitr)
library(ggplot2)
library(viridis)
# ggplot settings
geom_histogram <- function(...) ggplot2::geom_histogram(..., fill = viridis(10, alpha = 0.5)[8], show.legend = FALSE, bins = 20, color = "black")
geom_smooth <- function(...) ggplot2::geom_smooth(..., color = viridis(10, alpha = 0.5)[8])
geom_boxplot <- function(...) ggplot2::geom_boxplot(..., fill = viridis(10, alpha = 0.5)[7])
geom_pointrange <- function(...) ggplot2::geom_pointrange(..., show.legend = FALSE, color = viridis(10, alpha = 0.5)[7], size = 2)
plot_model <- function(...) sjPlot::plot_model(xy_mod, type = "diag", colors = viridis(10, alpha = 0.5)[7])
theme_set(theme_classic(base_size = 20))
```
::: {.alert .alert-info}
# Objetivo del manual {.unnumbered .unlisted}
- Aprender las principales herramientas de simulación de datos en R
- Comprender la utlidad de usar datos simulados para entender el comportamiento de las herramientas estadísticas
:::
------------------------------------------------------------------------
Paquetes a utilizar en este manual:
```{r}
# instalar/cargar paquetes
sketchy::load_packages(
c("ggplot2",
"viridis"
)
)
```
------------------------------------------------------------------------
# Cómo simular datos
## Generación de números aleatorios en R
La estadística nos permite inferir patrones en los datos. Solemos utilizar conjuntos de datos reales para enseñar estadística. Sin embargo, puede ser circular entender el funcionamiento interno de una herramienta estadística probando su capacidad para inferir un patrón que no estamos seguros de encontrar en los datos (y no tenemos idea del mecanismo que produjo ese patrón). **Las simulaciones nos permiten crear escenarios controlados en los que conocemos con seguridad los patrones** presentes en los datos y los procesos subyacentes que los han generado.
R ofrece algunas funciones básicas para la simulación de datos. Las más utilizadas son las funciones generadoras de números aleatorios. Los nombres de estas funciones comienzan con *r* (`r____()`). Por ejemplo, `runif()`:
```{r, echo=FALSE}
set.seed(7)
```
```{r}
# simular variable uniforme
unif_var <- runif(n = 100, min = 0, max = 10)
```
El resultado es un vector numérico de longitud 100 (`n = 100`):
```{r}
# imprimir variable
unif_var
```
Podemos explorar el resultado graficando un histograma:
```{r}
# crear histograma
ggplot(data = data.frame(unif_var), mapping = aes(x = unif_var)) + geom_histogram()
```
Muestra una distribución uniforme que va de 0 a 10.
También podemos simular números aleatorios procedentes de una distribución normal utilizando `rnorm()`:
```{r}
# crear una variable normal
norm_var <- rnorm(n = 1000, mean = 2, sd = 1)
# graficar histograma
ggplot(data = data.frame(norm_var), mapping = aes(x = norm_var)) + geom_histogram()
```
Tenga en cuenta que todas las funciones generadoras de números aleatorios tienen el argumento 'n', que determina la longitud del vector generado (es decir, el número de números aleatorios), además de algunos argumentos adicionales relacionados con parámetros específicos de la distribución.
Las variables continuas (es decir, los vectores numéricos) pueden convertirse en variables discretas (es decir, números enteros) simplemente redondeándolas:
```{r}
v1 <- rnorm(n = 5, mean = 10, sd = 3)
v1
round(x = v1, digits = 0)
```
::: {.alert .alert-info}
<font size="5">Ejercicio 1</font>
- ¿Qué hacen las funciones `rbinom()` y `rexp()`?
- Ejecútela y haga histogramas de sus resultados
- ¿Qué hacen los argumentos 'mean' y 'sd' en `rnorm()`? Juegue con diferentes valores y comprueba el histograma para hacerse una idea de su efecto en la simulación
:::
## Generación de variables categóricas
La forma más sencilla de generar variables categóricas es utilizar el vector de ejemplo `letters' (o`LETTERS') para asignar niveles de categoría. Podemos hacerlo utilizando la función `rep()`. Por ejemplo, el siguiente código crea un vector categórico (caracteres) con dos niveles, cada uno con 4 observaciones:
```{r}
rep(letters[1:2], each = 4)
```
También podemos replicar este patrón utilizando el argumento 'times'. Este código replica el vector anterior 2 veces:
```{r}
rep(letters[1:2], each = 4, times = 2)
```
Otra opción es simular una variable a partir de una distribución binomial y luego convertirla en un factor:
```{r}
# correr rbinom
binom_var <- rbinom(n = 50, size = 1, prob = 0.5)
binom_var
```
```{r}
# convertir a factor
categ_var <- factor(binom_var, labels = c("a", "b"))
categ_var
```
## Muestreo aleatorio
La otra herramienta importante de R para jugar con datos simulados es `sample()`. Esta función permite tomar muestras de tamaños específicos de vectores. Por ejemplo, tomemos el ejemplo del vector 'letters':
```{r}
letters
```
Podemos tomar una muestra de este vector como es:
```{r}
# tomar muestra
sample(x = letters, size = 10)
```
El argumento 'size' nos permite determinar el tamaño de la muestra. Tenga en cuenta que obtendremos un error si el tamaño es mayor que el propio vector:
```{r, error=TRUE}
sample(x = letters, size = 30)
```
Esto sólo puede hacerse cuando el muestreo es con reemplazo (replacement). El muestreo con reemplazo puede aplicarse estableciendo el argumento `replace = TRUE`:
```{r}
sample(x = letters, size = 30, replace = TRUE)
```
## Iterar un proceso
A menudo, las simulaciones deben repetirse varias veces para descartar resultados espurios debidos al azar o simplemente para probar diferentes parámetros. Las funciones de simulación de datos mencionadas anteriormente pueden ejecutarse varias veces (por ejemplo, iteradas) utilizando la función `replicate()`:
```{r}
# replicar
repl_rnorm <- replicate(n = 3, expr = rnorm(2), simplify = FALSE)
# ver clase
class(repl_rnorm)
# imprimir
repl_rnorm
```
## Hacer que las simulaciones sean reproducibles
El último truco que necesitamos para ejecutar simulaciones en R es la capacidad de reproducir una simulación (es decir, obtener exactamente los mismos datos y resultados simulados). Esto puede ser útil para que otros investigadores puedan ejecutar nuestros análisis exactamente de la misma manera. Esto puede hacerse fácilmente con la función `set.seed()`. Pruebe a ejecutar el siguiente código. Debería obtener la misma salida:
```{r}
# definir semilla
set.seed(10)
# crear variable uniforme
runif(n = 2)
```
------------------------------------------------------------------------
# Crear juegos de datos
## Juegos de datos con variables numéricas y categóricas
Ahora que sabemos cómo simular variables continuas y categóricas. Podemos juntarlas para crear conjuntos de datos simulados. Esto se puede hacer utilizando la función `data.frame()`:
```{r}
# crear variable categorica
grupo <- rep(letters[1:2], each = 3)
# crear variable continuaa
tamano <- rnorm(n = 6, mean = 5, sd = 1)
# poner juntas en un data frame
df <- data.frame(grupo, tamano)
# imprimir
df
```
Por supuesto, podríamos añadir más variables a este cuadro de datos:
```{r}
# crear variable categorica
grupo <- rep(letters[1:2], each = 3)
individuo <- rep(LETTERS[1:6])
# crear variables continuas
tamano <- rnorm(n = 6, mean = 5, sd = 1)
peso <- rnorm(n = 6, mean = 100, sd = 10)
# poner todo en un data frame
df <- data.frame(grupo, individuo, tamano, peso)
# imprimir
df
```
Y eso es un juego de datos simulados en su forma más básica. Se parece mucho al tipo de datos con los que trabajamos en biología.
------------------------------------------------------------------------
# Cómo utilizar datos simulados para entender el comportamiento de las herramientas estadísticas
## Prueba de concepto: *el Teorema del Límite Central*
El [Teorema del Límite Central](https://en.wikipedia.org/wiki/Central_limit_theorem) afirma que, si tomamos muestras aleatorias de una población, los promedios de esas muestras seguirán una distribución normal, aunque la población no esté distribuida normalmente. Además, la distribución normal resultante debe tener un promedio cercano al promedio de la población. El teorema es un concepto clave para la estadística inferencial, ya que implica que los métodos estadísticos que funcionan para las distribuciones normales pueden ser aplicables a muchos problemas que implican otros tipos de distribuciones. No obstante, el objetivo aquí es sólo mostrar cómo se pueden utilizar las simulaciones para entender el comportamiento de los métodos estadísticos.
Para comprobar si esas afirmaciones básicas sobre el Teorema del Límite Central son ciertas, podemos utilizar datos simulados en R. Vamos a simular una población de 1000 observaciones con una distribución uniforme:
```{r, eval=FALSE}
# simular popublacion uniforme
unif_pop <- runif(1000, min = 0, max = 10)
# ver histograma
ggplot(data = data.frame(unif_pop), mapping = aes(x = unif_pop)) + geom_histogram()
```
```{r, echo=FALSE}
set.seed(10)
# simulate uniform population
unif_pop <- runif(1000, min = 0, max = 10)
# ver distribucion
ggplot(data = data.frame(unif_pop), mapping = aes(x = unif_pop)) + geom_histogram()
```
Podemos tomar muestras aleatorias usando `sample()` así:
```{r}
sample(x = unif_pop, size = 30)
```
Este proceso puede ser replicado varias veces con `replicate()`:
```{r}
# replicar
samples <- replicate(n = 100, expr = mean(sample(x = unif_pop, size = 30)))
```
El código anterior toma 100 muestras con 30 valores cada una. Ahora podemos comprobar la distribución de las muestras:
```{r, eval=FALSE}
# ver distribucion/ histograma
ggplot(data = data.frame(samples), mapping = aes(x = samples)) + geom_histogram()
```
```{r, echo=FALSE}
# ver distribucion
ggplot(data = data.frame(samples), mapping = aes(x = samples)) + geom_histogram()
```
... asi como el promedio:
```{r}
mean(samples)
```
Como era de esperar, las muestras siguen una distribución normal con una media cercana a la media de la población, que es:
```{r}
mean(unif_pop)
```
Probemos con una distribución más compleja. Por ejemplo, una distribución bimodal:
```{r, eval=FALSE}
# usar semilla
set.seed(123)
# simular variables
norm1 <- rnorm(n = 1000, mean = 10, sd = 3)
norm2 <- rnorm(n = 1000, mean = 20, sd = 3)
# juntar en una sola variable
bimod_pop <- c(norm1, norm2)
# ver histograma
ggplot(data = data.frame(bimod_pop), mapping = aes(x = bimod_pop)) + geom_histogram()
```
```{r, echo=FALSE}
# definir semilla
set.seed(123)
norm1 <- rnorm(n = 1000, mean = 10, sd = 3)
norm2 <- rnorm(n = 1000, mean = 20, sd = 3)
bimod_pop <- c(norm1, norm2)
# ver distribucion
ggplot(data = data.frame(bimod_pop), mapping = aes(x = bimod_pop)) + geom_histogram()
```
```{r, eval=FALSE}
# replicar muestreo
samples <- replicate(200, mean(sample(bimod_pop, 10)))
# ver histograma
ggplot(data = data.frame(samples), mapping = aes(x = samples)) + geom_histogram()
```
```{r, echo=FALSE}
samples <- replicate(200, mean(sample(bimod_pop, 10)))
# ver distribucion
ggplot(data = data.frame(samples), mapping = aes(x = samples)) + geom_histogram()
```
```{r}
# ver promedios
mean(samples)
mean(bimod_pop)
```
::: {.alert .alert-info}
<font size="5">Ejercicio 2</font>
- Intenta explorar el Teorema del Límite Central como en el caso anterior, pero esta vez utilizando:
1. Una distribución exponencial (`rexp()`)
2. Una distribución log-normal (`rlnorm()`)
- Para cada distribución: grafique un histograma y compare los promedios de la población y de las muestras
:::
## Referencias
- [R's rbinom -- Simulate Binomial or Bernoulli trials](https://www.programmingr.com/examples/neat-tricks/sample-r-function/r-rbinom/)
- [R's rnorm -- selecting values from a normal distribution](https://www.programmingr.com/examples/neat-tricks/sample-r-function/r-rnorm/)
- [R's exp -- Simulating Exponential Distributions](https://www.programmingr.com/examples/neat-tricks/sample-r-function/rexp/)
- [Simulating data in R](https://aosmith.rbind.io/2018/08/29/getting-started-simulating-data/)
------------------------------------------------------------------------
# Información de la sesión {.unnumbered .unlisted}
```{r session info, echo=F}
sessionInfo()
```