8 Variables

R pre-introductorio enseña estadística. No mucha, pero es algo.

8.1 Cualidades y cantidades

Exploremos el conjunto de datos penguins (viene en el paquete palmerpenguins):

## Explore data
head(penguins)
#> # A tibble: 6 × 8
#>   species island bill_…¹ bill_…² flipp…³ body_…⁴ sex    year
#>   <fct>   <fct>    <dbl>   <dbl>   <int>   <int> <fct> <int>
#> 1 Adelie  Torge…    39.1    18.7     181    3750 male   2007
#> 2 Adelie  Torge…    39.5    17.4     186    3800 fema…  2007
#> 3 Adelie  Torge…    40.3    18       195    3250 fema…  2007
#> 4 Adelie  Torge…    NA      NA        NA      NA <NA>   2007
#> 5 Adelie  Torge…    36.7    19.3     193    3450 fema…  2007
#> 6 Adelie  Torge…    39.3    20.6     190    3650 male   2007
#> # … with abbreviated variable names ¹​bill_length_mm,
#> #   ²​bill_depth_mm, ³​flipper_length_mm, ⁴​body_mass_g

Ok, tiene ocho variables. Pongamos atención a cómo son las variables:

## Explore data
customGlimpse(penguins)
#>            col_name col_index col_class
#> 1           species         1    factor
#> 2            island         2    factor
#> 3    bill_length_mm         3   numeric
#> 4     bill_depth_mm         4   numeric
#> 5 flipper_length_mm         5   integer
#> 6       body_mass_g         6   integer
#> 7               sex         7    factor
#> 8              year         8   integer

En estadística, qué podemos hacerle a una variable depende enteramente de cuál sea su clase.

En este conjunto de datos tenemos unas variables numéricas (numeric, integer) y otras que son factores (factor). ¿A qué se debe esta diferencia? Pues se debe a la diferencia elemental entre variables cuantitativas (e.g., la edad, la estatura) y variables cualitativas (e.g., el género).

Revisemos un par de esas variables cuantitativas:

penguins |>
  select(bill_depth_mm, year) |>
  head()
#> # A tibble: 6 × 2
#>   bill_depth_mm  year
#>           <dbl> <int>
#> 1          18.7  2007
#> 2          17.4  2007
#> 3          18    2007
#> 4          NA    2007
#> 5          19.3  2007
#> 6          20.6  2007

R dice que bill_depth_mm es double y que year es integer; en efecto, arriba podemos observar que year tiene números enteros mientras bill_depth_mm presenta decimales.

Ya aprendimos algo nuevo: R tiene una clase específica para los números decimales: double. Si bien R guarda los enteros y los decimales en objetos diferentes, la mayoría del tiempo las mismas funciones se pueden usar en ambas variables porque, de hecho, R posee una clase que agrupa las variables double e integer: numeric.

Dado que bill_depth_mm y year son variables numéricas, podemos hacer con ellas lo que normalmente haríamos con un vector compuesto por números, por ejemplo:

## Compute
range(penguins$bill_depth_mm)
#> [1] NA NA

## My bad! Try again
range(penguins$bill_depth_mm, na.rm = TRUE)
#> [1] 13.1 21.5

De la función range() obtenemos el valor mínimo y el máximo de una variable cuantitativa, que individualmente se obtendrían con las funciones min() y max():

## Compute
range(penguins$year, na.rm = TRUE)
#> [1] 2007 2009

## Compute
min(penguins$year)
#> [1] 2007

## Compute
max(penguins$year)
#> [1] 2009

¿Qué pasa si le aplico range() a una variable cualitativa?

## Explore variable
car::some(penguins[, 2]) |> head(10)
#> # A tibble: 10 × 1
#>    island   
#>    <fct>    
#>  1 Torgersen
#>  2 Dream    
#>  3 Biscoe   
#>  4 Dream    
#>  5 Biscoe   
#>  6 Biscoe   
#>  7 Dream    
#>  8 Dream    
#>  9 Dream    
#> 10 Dream

## Compute
range(penguins$species, na.rm = TRUE)
#> Error in Summary.factor(structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, : 'range' not meaningful for factors

Pasa que nos da un error. (Aprovecharé para recalcar que esta es de las pocas veces en las que R arroja un mensaje de error que sí se entiende: “‘range’ not meaningful for factors”).

Ok, ¿entonces qué podemos analizar cuando se trata de una variable cualitativa?

Puesto que las variables cualitativas agrupan información en categorías (por eso se les llama categóricas), podemos analizarlas contando cuántas observaciones hay para cada uno de sus niveles:

## Count
table(penguins$species)
#> 
#>    Adelie Chinstrap    Gentoo 
#>       152        68       124

En R, la clase de las variables cualitativas es factor y la función table() cuenta cuántas observaciones corresponden a cada nivel:

## Inspect
class(penguins$sex)
#> [1] "factor"

## Count
table(penguins$sex)
#> 
#> female   male 
#>    165    168

El potencial de agrupar variables en categorías no lo poseen las variables cuantitativas. Si no me creen, observen qué pasa cuando le aplico la función table() a una variable numérica:

## Count
table(penguins$bill_length_mm)
#> 
#> 32.1 33.1 33.5   34 34.1 34.4 34.5 34.6   35 35.1 35.2 35.3 
#>    1    1    1    1    1    1    1    2    2    1    1    1 
#> 35.5 35.6 35.7 35.9   36 36.2 36.3 36.4 36.5 36.6 36.7 36.8 
#>    2    1    3    2    4    3    1    2    2    2    2    1 
#> 36.9   37 37.2 37.3 37.5 37.6 37.7 37.8 37.9 38.1 38.2 38.3 
#>    1    2    2    3    2    3    3    5    2    4    2    1 
#> 38.5 38.6 38.7 38.8 38.9   39 39.1 39.2 39.3 39.5 39.6 39.7 
#>    1    3    1    3    2    3    1    3    1    3    5    4 
#> 39.8 40.1 40.2 40.3 40.5 40.6 40.7 40.8 40.9   41 41.1 41.3 
#>    1    1    3    2    2    4    1    2    4    1    7    2 
#> 41.4 41.5 41.6 41.7 41.8   42 42.1 42.2 42.3 42.4 42.5 42.6 
#>    2    2    1    1    1    3    1    2    1    1    3    1 
#> 42.7 42.8 42.9 43.1 43.2 43.3 43.4 43.5 43.6 43.8   44 44.1 
#>    2    2    2    1    4    2    1    3    1    1    1    2 
#> 44.4 44.5 44.9   45 45.1 45.2 45.3 45.4 45.5 45.6 45.7 45.8 
#>    1    3    2    1    3    6    2    2    5    2    3    3 
#> 45.9   46 46.1 46.2 46.3 46.4 46.5 46.6 46.7 46.8 46.9   47 
#>    1    2    3    5    1    4    5    2    2    4    2    1 
#> 47.2 47.3 47.4 47.5 47.6 47.7 47.8 48.1 48.2 48.4 48.5 48.6 
#>    2    2    1    4    2    1    1    2    2    3    3    1 
#> 48.7 48.8   49 49.1 49.2 49.3 49.4 49.5 49.6 49.7 49.8 49.9 
#>    3    1    3    3    2    2    1    3    3    1    3    1 
#>   50 50.1 50.2 50.3 50.4 50.5 50.6 50.7 50.8 50.9   51 51.1 
#>    5    2    3    1    2    5    1    2    4    2    1    2 
#> 51.3 51.4 51.5 51.7 51.9   52 52.1 52.2 52.5 52.7 52.8 53.4 
#>    4    1    2    1    1    3    1    2    1    1    1    1 
#> 53.5 54.2 54.3 55.1 55.8 55.9   58 59.6 
#>    1    1    1    1    1    1    1    1

Evidentemente, el resultado del código anterior es un despropósito y no tiene sentido asumir que 47.5, 47.6, 47.7 y 47.8 son categorías específicas.

Afortunadamente, R nos da alternativas para convertir variables cuantitativas en factores, como la función cut():

## Prepare data
(vector <- c(1:10))
#>  [1]  1  2  3  4  5  6  7  8  9 10

## Cut
cut(vector, 
    breaks = c(0, median(vector), Inf),
    labels = c("low", "high"))
#>  [1] low  low  low  low  low  high high high high high
#> Levels: low high

Por algo el análisis exploratorio de datos es de máxima importancia: debemos conocer nuestras variables antes de trabajarlas de cualquier modo, porque de hecho qué podremos hacer con ellas depende enteramente de cómo son.

La función skim() del paquete skimr es fenomenal para realizar análisis exploratorio de datos:

## Explore data
skim(penguins) |> kable()

skim_type	skim_variable	n_missing	complete_rate	factor.ordered	factor.n_unique	factor.top_counts	numeric.mean	numeric.sd	numeric.p0	numeric.p25	numeric.p50	numeric.p75	numeric.p100	numeric.hist
factor	species	0	1.0000000	FALSE	3	Ade: 152, Gen: 124, Chi: 68	NA	NA	NA	NA	NA	NA	NA	NA
factor	island	0	1.0000000	FALSE	3	Bis: 168, Dre: 124, Tor: 52	NA	NA	NA	NA	NA	NA	NA	NA
factor	sex	11	0.9680233	FALSE	2	mal: 168, fem: 165	NA	NA	NA	NA	NA	NA	NA	NA
numeric	bill_length_mm	2	0.9941860	NA	NA	NA	43.92193	5.4595837	32.1	39.225	44.45	48.5	59.6	▃▇▇▆▁
numeric	bill_depth_mm	2	0.9941860	NA	NA	NA	17.15117	1.9747932	13.1	15.600	17.30	18.7	21.5	▅▅▇▇▂
numeric	flipper_length_mm	2	0.9941860	NA	NA	NA	200.91520	14.0617137	172.0	190.000	197.00	213.0	231.0	▂▇▃▅▂
numeric	body_mass_g	2	0.9941860	NA	NA	NA	4201.75439	801.9545357	2700.0	3550.000	4050.00	4750.0	6300.0	▃▇▆▃▂
numeric	year	0	1.0000000	NA	NA	NA	2008.02907	0.8183559	2007.0	2007.000	2008.00	2009.0	2009.0	▇▁▇▁▇

En realidad ese |> kable() lo necesité por motivos particulares (compilación de los archivos que dan forma a este ebook) que no vale la pena detallar. Ustedes no deberían tener problema en usarla así, sin más:

## Explore data
skim(penguins) 

¡Es imposible no amar esta función! skimr::skim() nos entrega una panorámica del conjunto de datos: número de filas y de columnas, frecuencia de las variables según sean cuantitativas o cualitativas, valores faltantes (NA), si los niveles están ordenados o no, cuántas observaciones hay por cada nivel (en el caso de las variables cualitativas), y un recuento pormenorizado de las variables cuantitativas que incluye medidas descriptivas como la media, la desviación estandar, la mediana, etcétera. Es tan completa esta función que hasta dibuja un histograma en la última columna.

Otra opción excelente es describe() del paquete psych:

psych::describe(penguins) |> kable()

	vars	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
species*	1	344	1.918605	0.8933198	2.00	1.898551	1.48260	1.0	3.0	2.0	0.1591315	-1.7318773	0.0481646
island*	2	344	1.662791	0.7261940	2.00	1.579710	1.48260	1.0	3.0	2.0	0.6086049	-0.9064333	0.0391538
bill_length_mm	3	342	43.921930	5.4595837	44.45	43.906934	7.04235	32.1	59.6	27.5	0.0526530	-0.8931397	0.2952205
bill_depth_mm	4	342	17.151170	1.9747932	17.30	17.172628	2.22390	13.1	21.5	8.4	-0.1422086	-0.9233523	0.1067846
flipper_length_mm	5	342	200.915205	14.0617137	197.00	200.335766	16.30860	172.0	231.0	59.0	0.3426554	-0.9991866	0.7603704
body_mass_g	6	342	4201.754386	801.9545357	4050.00	4154.014598	889.56000	2700.0	6300.0	3600.0	0.4662117	-0.7395200	43.3647348
sex*	7	333	1.504504	0.5007321	2.00	1.505618	0.00000	1.0	2.0	1.0	-0.0179376	-2.0056743	0.0274400
year	8	344	2008.029070	0.8183559	2008.00	2008.036232	1.48260	2007.0	2009.0	2.0	-0.0532601	-1.5092478	0.0441228

Una contribución fundamental de la estadística tiene que ver con la descripción de nuestras variables. Afortunadamente, describir variables cualitativas está bien fácil: basta con contar cuántas observaciones hay en cada nivel, tal como haremos a continuación con el data frame diamonds:

## Explore data
head(diamonds[, 1:7])
#> # A tibble: 6 × 7
#>   carat cut       color clarity depth table price
#>   <dbl> <ord>     <ord> <ord>   <dbl> <dbl> <int>
#> 1  0.23 Ideal     E     SI2      61.5    55   326
#> 2  0.21 Premium   E     SI1      59.8    61   326
#> 3  0.23 Good      E     VS1      56.9    65   327
#> 4  0.29 Premium   I     VS2      62.4    58   334
#> 5  0.31 Good      J     SI2      63.3    58   335
#> 6  0.24 Very Good J     VVS2     62.8    57   336

## Count
table(diamonds$cut)
#> 
#>      Fair      Good Very Good   Premium     Ideal 
#>      1610      4906     12082     13791     21551

Y listo… Aunque por allá un bar plot no va mal:

## Plot
ggplot(diamonds, aes(x = cut, fill = cut)) +
  geom_bar() +
  scale_fill_brewer(palette = "Reds") 

Desde luego que me quedó espantoso, pero si ustedes le dedican tiempo lo podrán mejorar montones.

8.2 Media

La descripción de variables cuantitativas sí nos da más oficio. Primero, miremos cómo se ve una variable cuantitativa:

## Explore data
head(penguins[, 6], 20)
#> # A tibble: 20 × 1
#>    body_mass_g
#>          <int>
#>  1        3750
#>  2        3800
#>  3        3250
#>  4          NA
#>  5        3450
#>  6        3650
#>  7        3625
#>  8        4675
#>  9        3475
#> 10        4250
#> 11        3300
#> 12        3700
#> 13        3200
#> 14        3800
#> 15        4400
#> 16        3700
#> 17        3450
#> 18        4500
#> 19        3325
#> 20        4200

Arriba tomamos un pedacito, apenas unas pocas filas. Ahora vamos a verla toda:

## Pull
penguins$body_mass_g
#>   [1] 3750 3800 3250   NA 3450 3650 3625 4675 3475 4250 3300
#>  [12] 3700 3200 3800 4400 3700 3450 4500 3325 4200 3400 3600
#>  [23] 3800 3950 3800 3800 3550 3200 3150 3950 3250 3900 3300
#>  [34] 3900 3325 4150 3950 3550 3300 4650 3150 3900 3100 4400
#>  [45] 3000 4600 3425 2975 3450 4150 3500 4300 3450 4050 2900
#>  [56] 3700 3550 3800 2850 3750 3150 4400 3600 4050 2850 3950
#>  [67] 3350 4100 3050 4450 3600 3900 3550 4150 3700 4250 3700
#>  [78] 3900 3550 4000 3200 4700 3800 4200 3350 3550 3800 3500
#>  [89] 3950 3600 3550 4300 3400 4450 3300 4300 3700 4350 2900
#> [100] 4100 3725 4725 3075 4250 2925 3550 3750 3900 3175 4775
#> [111] 3825 4600 3200 4275 3900 4075 2900 3775 3350 3325 3150
#> [122] 3500 3450 3875 3050 4000 3275 4300 3050 4000 3325 3500
#> [133] 3500 4475 3425 3900 3175 3975 3400 4250 3400 3475 3050
#> [144] 3725 3000 3650 4250 3475 3450 3750 3700 4000 4500 5700
#> [155] 4450 5700 5400 4550 4800 5200 4400 5150 4650 5550 4650
#> [166] 5850 4200 5850 4150 6300 4800 5350 5700 5000 4400 5050
#> [177] 5000 5100 4100 5650 4600 5550 5250 4700 5050 6050 5150
#> [188] 5400 4950 5250 4350 5350 3950 5700 4300 4750 5550 4900
#> [199] 4200 5400 5100 5300 4850 5300 4400 5000 4900 5050 4300
#> [210] 5000 4450 5550 4200 5300 4400 5650 4700 5700 4650 5800
#> [221] 4700 5550 4750 5000 5100 5200 4700 5800 4600 6000 4750
#> [232] 5950 4625 5450 4725 5350 4750 5600 4600 5300 4875 5550
#> [243] 4950 5400 4750 5650 4850 5200 4925 4875 4625 5250 4850
#> [254] 5600 4975 5500 4725 5500 4700 5500 4575 5500 5000 5950
#> [265] 4650 5500 4375 5850 4875 6000 4925   NA 4850 5750 5200
#> [276] 5400 3500 3900 3650 3525 3725 3950 3250 3750 4150 3700
#> [287] 3800 3775 3700 4050 3575 4050 3300 3700 3450 4400 3600
#> [298] 3400 2900 3800 3300 4150 3400 3800 3700 4550 3200 4300
#> [309] 3350 4100 3600 3900 3850 4800 2700 4500 3950 3650 3550
#> [320] 3500 3675 4450 3400 4300 3250 3675 3325 3950 3600 4050
#> [331] 3350 3450 3250 4050 3800 3525 3950 3650 3650 4000 3400
#> [342] 3775 4100 3775

Ese chorizo de números ilustra que las variables cuantitativas tienen una distribución de valores. Como vimos antes, esos valores no se pueden agrupar en categorías. Tendremos que describir esa variable de alguna otra forma.

Al describir la distribución de valores de una variable, nuestro primer objetivo es hallar un único valor que sea la mejor representación de toda esa distribución. Ese valor más representativo es la media, es decir, el promedio de todos los valores.

La función mean() nos devuelve el promedio de una variable:

## Compute
mean(penguins$body_mass_g)
#> [1] NA

No funcionó… Y no funcionó por una razón conocida: R no computa promedios si los datos tienen valores ausentes (NA).

La misma función mean() trae consigo un argumento que permite arreglar este problema, si lo configuramos adecuadamente (na.rm = TRUE):

## Compute
mean(penguins$body_mass_g, na.rm = TRUE) 
#> [1] 4201.754

Ya obtuvimos el promedio: el valor que mejor representa al chorizo de números que desplegamos hace un ratito.

¿Qué es la media entonces? Ya dijimos que la media es el promedio de la distribución de valores de una variable. Primero los sumamos todos los valores juntos (numerador) y luego dividimos el resultado por el número de valores (denominador):

## Add values
(num <- sum(penguins$body_mass_g, na.rm = TRUE))
#> [1] 1437000

## Get the number of observations
(den <- length(penguins$body_mass_g))
#> [1] 344

## Compute the mean (divide)
num / den
#> [1] 4177.326

Si habíamos calculado que la media es 4201.754386, ¿por qué no nos dio igual? La respuesta es: porque en la función length() no desconté los NA:

## Verify
mean(penguins$body_mass_g, na.rm = TRUE) == (num / den)
#> [1] FALSE

La función length() no dispone de un argumento para remover los NA, a diferencia de mean() que, como hemos visto, sí lo tiene. Esta es de esas ocasiones en que R nos obliga a dar un vuelco, a buscar otra función. ¿Existe? Sí, sí existe.

Podemos usar la función na.omit() para remover los NA:

## Add values
num <- sum(penguins$body_mass_g, na.rm = TRUE)

## Get the number of values
den <- length(na.omit(penguins$body_mass_g)) # na.omit()

## Compute the mean 
num / den
#> [1] 4201.754

## Compute
mean(penguins$body_mass_g, na.rm = TRUE) 
#> [1] 4201.754

Ahora sí nos dio consistente. Por supuesto que todo este código que crea num y den y luego calcula num / den era innecesario en primer lugar si de por sí R ya nos daba una salida más fácil y directa:

## Compute
mean(penguins$body_mass_g, na.rm = TRUE)
#> [1] 4201.754

Ese ejemplo, sin embargo, ilustra a la perfección cómo podemos capitalizar R para estudiar conceptos estadísticos mediante la descomposición de funciones sencillas como mean() en una serie de comandos que acaba con el mismo resultado pero nos permite examinar el concepto en cuestión paso a paso. Descomponer código es algo que desde luego no realizaríamos en un análisis para el que nos contraten, pero sí podemos intentarlo a la hora de estudiar.

Gracias a la descomposición de código pude probarme a mí mismo que el promedio es, en efecto, la suma de los valores que toma la variable dividido entre el número de observaciones que corresponden a la variable:

## Compute 
mean(penguins$body_mass_g, na.rm = TRUE) 
#> [1] 4201.754

## Ok, let's prove it
sum(penguins$body_mass_g, na.rm = TRUE) / length(na.omit(penguins$body_mass_g))
#> [1] 4201.754

Es muy fácil esto.

8.3 Varianza y desviación estándar

La media es una medida de centralidad: si tenemos que hallar el número que mejor describe a la distribución de una variable, la media es usualmente el mejor candidato (salvo si la variable tienen valores atípicos, en cuyo caso conviene usar la mediana, una medida más robusta).

La media es el valor al que más se parecen, en promedio, los valores en cuestión, el que mejor los describe (siempre que no haya outliers, tal como acabo de advertir).

Ahora digamos que también nos gustaría saber cuán variable es una variable.

Preguntarse cuán-variable-es-una-variable tiene sentido y les voy a demostrar por qué:

## Compute
mean(1, 1, 1, 1, 1, 1, 1)
#> [1] 1

## Compute
mean(1, 2, 3, 4, -4, 0)
#> [1] 1

En estos dos casos la media es la misma, pero en el primer caso todos los valores son iguales (¡no varía la variable!) y en el segundo ninguno se repite.

Resulta que la respuesta a cuán-variable-es-una-variable la conseguimos también mediente un promedio: primero calculamos la distancia entre cada uno de los valores respecto a la media, y luego sacamos el promedio de esas distancias.

Entonces, si los valores son c(1, 2, 3, 4 -1, 0), a cada uno le voy a restar la media, que es igual a mean(1, 2, 3, 4 -1, 0), que es igual a 1, y finalmente calculo el promedio de las distancias.

Así se hace:

## Compute (value - mean)
1 - mean(1, 2, 3, 4, -4, 0)
#> [1] 0
2 - mean(1, 2, 3, 4, -4, 0)
#> [1] 1
3 - mean(1, 2, 3, 4, -4, 0)
#> [1] 2
4 - mean(1, 2, 3, 4, -4, 0)
#> [1] 3
-4 - mean(1, 2, 3, 4, -4, 0)
#> [1] -5
0 - mean(1, 2, 3, 4, -4, 0)
#> [1] -1

Acabo de calcular por separado las distancias de cada valor respecto a la media. Pero lo que yo necesito es un promedio, es decir, tengo que sumarlos todos y dividirlos por cuántos son.

Primero voy a sumar todas las distancias:

## Compute
1 - mean(1, 2, 3, 4, -4, 0) +
  2 - mean(1, 2, 3, 4, -4, 0) +
  3 - mean(1, 2, 3, 4, -4, 0) +
  4 - mean(1, 2, 3, 4, -4, 0) +
  -4 - mean(1, 2, 3, 4, -4, 0) +
  0 - mean(1, 2, 3, 4, -4, 0) 
#> [1] 0

Ok, acá hay un problema: si nada más sumo las distancias de cada valor respecto a la media, siempre voy a obtener cero.

Para que eso no suceda, debo dar un paso intermedio: elevar al cuadrado las distancias antes de sumarlas todas juntas, de modo tal que, al convertirlas en números positivos todas, ya no se van a cancelar entre sí cuando las sume:

## Compute
((1 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
  ((2 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
  ((3 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
  ((4 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
  ((-4 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
  ((0 - mean(1, 2, 3, 4, -4, 0)) ^ 2)
#> [1] 40

Y ahora que ya las logramos sumar, vamos a dividir ese total entre el número de valores:

## Compute
(((1 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
   ((2 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
   ((3 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
   ((4 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
   ((-4 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
   ((0 - mean(1, 2, 3, 4, -4, 0)) ^ 2)) / 
  length(c(1, 2, 3, 4, -4, 0))
#> [1] 6.666667

Y listo, ese que obtuvimos ahí es el promedio del cuadrado de las distancias entre cada uno de los valores respecto a la media, also known as varianza.

Como extrema precaución, vamos a calcular la varianza en R usando la función var(), así verificaremos el resultado que obtuvimos por descomposición de código:

## Compute
var(c(1, 2, 3, 4, -4, 0))
#> [1] 8

Y resulta que no nos dio lo mismo…

El problema esta vez es que el código no es consistente con la fórmula matemática para calcular la varianza: la varianza es el promedio del cuadrado de las distancias entre cada uno de los valores respecto a la media, eso es correcto conceptualmente y sí lo había explicado bien, pero el promedio se calcula tomando en cuanta no el número \(N\) de valores sino \(N - 1\):

## Compute
(((1 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
   ((2 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
   ((3 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
   ((4 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
   ((-4 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
   ((0 - mean(1, 2, 3, 4, -4, 0)) ^ 2)) / 
  (length(c(1, 2, 3, 4, -4, 0)) - 1) ## N - 1
#> [1] 8

Ya ahora sí la pegamos: nuestra descomposición de código para calcular la varianza coincidió con el resultado de la función var().

Recapitulemos: la varianza es una medida de cuán variable es una variable.

Calculemos ahora la varianza de aquel otro caso en el que todos los valores eran iguales:

## Compute
var(c(1, 1, 1, 1, 1, 1))
#> [1] 0

## Ok, let's prove it!
(((1 - mean(1, 1, 1, 1, 1, 1)) ^ 2) +
    ((1 - mean(1, 1, 1, 1, 1, 1)) ^ 2) +
    ((1 - mean(1, 1, 1, 1, 1, 1)) ^ 2) +
    ((1 - mean(1, 1, 1, 1, 1, 1)) ^ 2) +
    ((1 - mean(1, 1, 1, 1, 1, 1)) ^ 2) +
    ((1 - mean(1, 1, 1, 1, 1, 1)) ^ 2)) / 
  (length(c(1, 1, 1, 1, 1, 1)) - 1)
#> [1] 0

No hay ninguna sorpresa aquí: si todos los valores son iguales, evidentemente, no hay variación alguna. Una varianza igual a cero significa que la variable no varía (es más, no es una variable sino una constante)… Este es un caso extremo, que he utilizado tan sólo para demostrar que sí tiene sentido preguntarse cuán variable es una variable.

Ahora quiero destacar que para calcular la varianza elevamos al cuadrado las distancias entre cada valor respecto a la media, y luego las sumamos, y luego las dividimos entre las \(N - 1\) observaciones. Pero enfoquémonos en el hecho de que las elevamos al cuadrado: ese paso nos permitió llevar a cabo la suma de las distancias (ya vimos que, si no las elevábamos al cuadrado, esas distancias se cancelaban entre sí y resultaban en cero); no obstante, elevarlas al cuadrado también nos privó de la posibilidad de interpretar la varianza en las unidades originales de la variable (que en este caso no existían porque eran puros números que me inventé, números sin significado de fondo, pero obviamente lo usual es que nuestros datos sí tengan significado y, por ende, unidades de medida).

Para, digamos, “regresar” la varianza a las unidades originales de la variable basta con sacarle la raíz cuadrada:

## Compute
sqrt(
  (((1 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
     ((2 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
     ((3 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
     ((4 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
     ((-4 - mean(1, 2, 3, 4, -4, 0)) ^ 2) +
     ((0 - mean(1, 2, 3, 4, -4, 0)) ^ 2)) / 
    (length(c(1, 2, 3, 4, -4, 0)) - 1)
)
#> [1] 2.828427

Plot twist: el resultado de sacarle la raíz cuadrada a la varianza es, de hecho, la desviación estandar.

R calcula la desviación estandar con la función sd():

## Compute
sd(c(1, 2, 3, 4, -4, 0))
#> [1] 2.828427

En síntesis:

• Varianza: desviación cuadrada promedio respecto a la media.

• Desviación estandar: desviación promedio respecto a la media.

Volviendo al data frame de pingüinos, ahora estamos en condiciones de efectuar otros cálculos además de la media:

## Compute the mean
mean(penguins$body_mass_g, na.rm = TRUE) 
#> [1] 4201.754

## Compute the standard deviation
sd(penguins$body_mass_g, na.rm = TRUE) 
#> [1] 801.9545

## Compute the variance
var(penguins$body_mass_g, na.rm = TRUE) 
#> [1] 643131.1

Para que se formen una mejor idea de cómo la varianza arruina las unidades originales de la variable, sólo compárenla con el valor máximo y miren cómo lo rebasa:

## Get the maximum value
max(penguins$body_mass_g, na.rm = TRUE) 
#> [1] 6300

## Compute the variance
var(penguins$body_mass_g, na.rm = TRUE) 
#> [1] 643131.1

La desviación estandar corrige ese comportamiento de la varianza y nos regresa al rango de la variable:

## Get the maximum value
max(penguins$body_mass_g, na.rm = TRUE) 
#> [1] 6300

## Compute the standard deviation
sd(penguins$body_mass_g, na.rm = TRUE) 
#> [1] 801.9545

Hasta el momento hemos trabajado con nada más una variable a la vez. De aquí en adelante vamos a trabajar con dos variables.

Cuando analizamos dos variables juntas nuestro objetivo es determinar si existe relación entre ellas. Para saberlo dispondremos de tres mediciones diferentes: la covarianza, el coeficiente de correlación y la pendiente de la línea de regresión.

Más información, en el próximo capítulo.