Fecha de publicación

12 de diciembre de 2024

 

Objetivo del manual

  • 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:

Código
# instalar/cargar paquetes

sketchy::load_packages(
  c("ggplot2", 
    "viridis"
    )
  )

 


1 Cómo simular datos

1.1 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():

Código
# simular variable uniforme
unif_var <- runif(n = 100, min = 0, max = 10)

 

El resultado es un vector numérico de longitud 100 (n = 100):

Código
# imprimir variable
unif_var
  [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:

Código
# 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():

Código
# 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:

Código
v1 <- rnorm(n = 5, mean = 10, sd = 3)

v1
[1] 15.4269  8.5450 15.1387  8.9203  9.3693
Código
round(x = v1, digits = 0)
[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

 

1.2 Generación de variables categóricas

La forma más sencilla de generar variables categóricas es utilizar el vector de ejemplo letters' (oLETTERS’) 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:

Código
rep(letters[1:2], each = 4)
[1] "a" "a" "a" "a" "b" "b" "b" "b"

 

También podemos replicar este patrón utilizando el argumento ‘times’. Este código replica el vector anterior 2 veces:

Código
rep(letters[1:2], each = 4, times = 2)
 [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:

Código
# correr rbinom
binom_var <- rbinom(n = 50, size = 1, prob = 0.5)

binom_var
 [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
Código
# convertir a factor
categ_var <- factor(binom_var, labels = c("a", "b"))

categ_var
 [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

 

1.3 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’:

Código
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:

Código
# tomar muestra
sample(x = letters, size = 10)
 [1] "i" "r" "q" "a" "c" "m" "y" "z" "u" "v"

 

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:

Código
sample(x = letters, size = 30)
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:

Código
sample(x = letters, size = 30, 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"

 

1.4 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():

Código
# replicar
repl_rnorm <- replicate(n = 3, expr = rnorm(2), simplify = FALSE)

# ver clase
class(repl_rnorm)
[1] "list"
Código
# imprimir
repl_rnorm
[[1]]
[1]  1.19184 -0.33992

[[2]]
[1]  0.78911 -0.63213

[[3]]
[1] -1.49312 -0.13441

 

1.5 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:

Código
# definir semilla
set.seed(10)

# crear variable uniforme
runif(n = 2)
[1] 0.50748 0.30677

2 Crear juegos de datos

2.1 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():

Código
# 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
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:

Código
# 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
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.


3 Cómo utilizar datos simulados para entender el comportamiento de las herramientas estadísticas

 

3.1 Prueba de concepto: el Teorema del Límite Central

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:

Código
# 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() 

 

Podemos tomar muestras aleatorias usando sample() así:

Código
sample(x = unif_pop, size = 30)
 [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():

Código
# 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:

Código
# ver distribucion/ histograma
ggplot(data = data.frame(samples), mapping = aes(x = samples)) + geom_histogram() 

 

… asi como el promedio:

Código
mean(samples)
[1] 5.0212

 

Como era de esperar, las muestras siguen una distribución normal con una media cercana a la media de la población, que es:

Código
mean(unif_pop)
[1] 5.0527

 

Probemos con una distribución más compleja. Por ejemplo, una distribución bimodal:

Código
# 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() 

Código
# replicar muestreo
samples <- replicate(200, mean(sample(bimod_pop, 10)))

# ver histograma
ggplot(data = data.frame(samples), mapping = aes(x = samples)) + geom_histogram() 

Código
# ver promedios
mean(samples)
[1] 15.231
Código
mean(bimod_pop)
[1] 15.088

 

Ejercicio 2

 

  • 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

3.2 Referencias


Información de la sesión

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