Apuntes de R

Apunte 05: Distribuciones e inferencias

Distribuciones de probabilidad

Las distribuciones de probabilidad son definidas en estadística en términos de funciones matemáticas que describen las probabilidades asociadas a cada posible resultado en el espacio muestral.

Dado que describe la probabilidad de todos los posibles resultados, la suma de las probabilidades en una distribución de probabilidad es siempre igual a 1. Otra forma de expresar esto es que la suma de las probabilidades en una distribución de probablidad discreta o la suma del área bajo la curva en una distribución de probabilidad continua considera el 100% de los posibles casos.

Existen distintos tipos de distribuciones de probabilidad.

Distribuciones discretas

Podemos generar distribuciones de probabilidad para eventos discretos, como la distribución binomial, que nos permite describir eventos con dos posibles resultados (p.e. lanzar una moneda):

dbinom(x, size, prob)

Donde: x = vector de cantidades, size = número de repeticiones, prob = probabilidad de éxito en cada repetición

exitos <- 0:20
plot(exitos, dbinom(exitos, size=20, prob=.3),type='h')

O una distribución uniforme, la cual describe eventos con múltiples resultados en donde todos tienen la misma probabilidad (p.e. lanzar un dado justo):

x <- rep(1:6,each = 20)
barplot(prop.table(table(x)))

O una distribución de Poisson, la que nos permite modelar la probabilidad de que ciertos eventos ocurran o no en cierto intervalo, por ejemplo, errores de tipeo en una página de un libro o accidentes de tránsito en un periodo de tiempo:

exitos <- 0:20
plot(exitos, dpois(exitos, lambda=3), type='h')

Distribuciones continuas

Las distribuciones continuas modelan variables que pueden tomar un rango infinito de valores dentro de un intervalo.

Una de ellas es la distribución exponencial, la que nos permite modelar el tiempo que falta para la ocurrencia de un evento, por ejemplo un terremoto.

x <- seq(0,20,0.01)
plot(x, dexp(x,rate = 2), lwd=2, col="black", type = "l")

O la distribución de Weibull, que es usada para modelar cuánto tiempo pasará antes de que falle un artefacto o componente.

x <- seq(0,5,0.01)
plot(x, dweibull(x,shape = 2), lwd=2, col="black", type = "l")

Existe también la distribución t, una distribución de probabilidad que se utiliza principalmente en inferencia estadística cuando el tamaño de la muestra es pequeño o cuando la desviación estándar de la población es desconocida.

dt(x, df)

Donde: x = vector de cantidades, df = grados de libertad

x <- seq(-4,4,0.01)
plot(x, dt(x,df = 5), lwd=2, col="black", type = "l")

Otro ejemplo es la distribución continua usada en estadística es chi-cuadrado, que se emplea para analizar la variabilidad en muestras de datos y realizar pruebas de hipótesis.

x <- seq(0,20,0.01)
plot(x, dchisq(x,df = 5), lwd=2, col="black", type = "l")

Similarmente, la distribución F es una distribución continua que se utiliza principalmente en análisis de varianza (ANOVA) y pruebas de regresión.

x <- seq(0,10,0.01)
plot(x, df(x,df1 = 5, df2 = 10), lwd=2, col="black", type = "l")

Y, por supuesto la distribución normal, también es una distribución continua que se utiliza ampliamente en estadística, y se caracteriza por ser simétrica en forma de campana.

x <- seq(-5,5,0.01)
plot(x, dnorm(x), lwd=2, col="black", type = "l")

La distribución normal (o las distribuciones normales)

La distribución normal, como se mencionó antes, es continua y modela muchos fenómenos naturales. Una de sus características más importantes es que, gracias al teorema del límite central, cuando se suman un gran número de variables aleatorias independientes, la distribución de esa suma se aproxima a una distribución normal.

Otro punto central a recordar es que la forma y posición de una distribución normal variará de acuerdo a dos parámetros: la media y la desviación estándar.

Inferencias desde la muestra

Para hacer inferencias a respecto a una población a partir de una muestra, podemos apoyarnos en las propiedad del teorema del límite central.

Esto se consigue usando lo que sabemos respecto a la población (media = \(\mu\) y desviación estándar = \(\sigma\)) y lo que sabemos sobre la muestra (media = \(\overline{x}\), desviación estándar = \(s\), tamaño muestra = \(n\)) a traves del error estandar:

\[ E.E. = \frac{\sigma}{\sqrt{n}} \]

Dada la naturaleza de la distribución normal de la muestra y la manera en que se distribuyen las proporciones bajo su curva, podemos utilizar este conocimiento para hacer predicciones sobre qué proporción de muestras contendrán el verdadero promedio poblacional cuando construimos un intervalo de confianza. En otras palabras, podemos utilizar la distribución normal para determinar la probabilidad de que nuestro intervalo contenga el verdadero valor del parámetro de interés.

Para ilustrar esto vamos a instalar una librería llamada TeachingDemos que incluye funciones para ilustrar esto.

#install.packages("TeachingDemos")

library(TeachingDemos)

Esta librería cuenta con la función ci.examp, que nos permite generar simulaciones para ver un número de muestras generadas al azar y ver si estas contienen o no al parámetro poblacional.

Para que siempre podamos reproducir estos resultados vamos a fijar la semilla que usamos para generar resultados (pseudo) aleatorios en R usando el comando set.seed y dándole un número. Si usas el comando set.seed con el mismo valor que otra persona en tu computador, vas a reproducir exactamente la misma muestra al azar que esa persona, y por ende generarán el mismo gráfico.

Después de fijar la semilla, se usa el comando ci.examp (confidence interval example) para generar un ejemplo de 20 muestras al azar y ver cuantas de ellas capturan el promedio poblacional. El comando ci.examp tiene entre sus argumentos n: tamaño muestral, conf.level:grado de confianza) y reps:número de muestras que queremos visualizar.

set.seed(2023)
ci.examp(reps = 20, n = 25, conf.level = 0.95)

Como podemos observar, 19 de las muestras capturan el promedio poblacional (100, en este ejemplo) y solo 1 muestra, la número 20, no lo captura. Esto es exactamente lo que se espera dado el grado de confianza de 95%.

Si disminuyera el grado de confianza al 50%, generando por ende intervalos más pequeños, más muestras fallarán en capturar el promedio.

set.seed(2023)
ci.examp(reps = 20, n = 25, conf.level = 0.5)

En este caso, obtenemos que solo 10 de las 20 muestras capturaron el promedio poblacional.

Si, en cambio, el nivel de confianza aumentara al 99%…

set.seed(2024)
ci.examp(reps = 100, n = 25, conf.level = 0.99)

Entonces, solo 1 de las 100 muestras falla en capturar el promedio poblacional. Sin embargo, podemos ver que nuestros intervalos de confianza son bastante grandes, abarcando un rango de aproximadamente 10 puntos.

¿Cómo podriamos reducir el tamaño de estos intervalos?

Es posible reducir el tamaño de los intervalos si incrementamos el otro componente del error estándar, a saber, el tamaño muestral:

set.seed(2024)
ci.examp(reps = 100, n = 100, conf.level = 0.99)

En este caso, hemos cuadruplicado nuestro tamaño muestral, lo que ha resultado en reducir los intervalos de confianza aproximadamente a la mitad, cubriendo cerca de 5 puntos. Esta relación no es un accidente, dado que en la formula del error estándar:

\[ E.E. = \frac{\sigma}{\sqrt{n}} \]

La desviación estándar poblacional se divide por la raíz cuadrada del tamaño muestral, por lo que necesitamos un aumento de 4 veces el tamaño muestral para reducir a la mitad el error estándar. Por consiguiente, si ahora queremos reducir los intervalos de confianza a 2.5 puntos, entonces vamos a necesitar una muestra de 400 personas:

set.seed(2025)
ci.examp(reps = 100, n = 400, conf.level = 0.99)

Es importante recordar esto para entender que la reducción del error estándar en base al aumento de muestras tiene un rendimiento decreciente.

El ejemplo a continuación es para un caso en que la desviación estándar poblacional es igual a 10:

curve(10/sqrt(x),1,1024, axes = FALSE, xlab = "Tamaño muestral", ylab = "Error estándar")
axis(1)
axis(2, las = 1)
title("Relacion entre error estándar y tamaño muestral")

points(1,10)
points(4,5)
points(16,2.5)
points(64,1.25) 
points(256,.625)

text(1,10, "Error estandar 10 con 1 persona", pos = 4)
text(4,5, "Error estandar 5 con 4 personas", pos = 4)
text(16,2.5, "Error estandar 2.5 con 16 personas", pos = 4)
text(64,1.25, "Error estandar 1.25 con 64 personas", pos = 4)
text(256,.625, "Error estandar .625 con 256 personas", pos = 4)

Calculando el intervalo de confianza de un promedio muestral (con \(\sigma\) conocido)

El un intervalo de confianza se estima identificando un margen de error en torno a un estimador puntual. Concretamente, se obtiene un extremo del intervalo restando el margen de error al estimador puntual y el otro extremo sumando el margen de error al estimador puntual.

\[ \text{Límite inferior} = \text{Estimador puntual} - \text{Valor Z para un grado de confianza} \times \text{Error Estándar} \]

\[ \text{Límite superior} = \text{Estimador puntual} + \text{Valor Z para un grado de confianza} \times \text{Error Estándar} \]

O escrito de forma más sintética:

\[ \text{Estimador puntual} \pm \text{Valor Z para un grado de confianza} \times \text{Error Estándar} \]

Si se nos dice que a partir de una población tiene un promedio \(\mu = 125\) y una desviación estándar \(\sigma = 20\) se obtiene una muestra de 25 personas con una media \(\overline{x} = 100\) y una desviación estándar \(s = 18\)… ¿Cómo se calcularía un intervalo de confianza al 95% del promedio poblacional?

# Datos relevantes
sigma = 20
n = 25
x_barra <- 100

# Calcular el valor Z para un 95%
# Estoy pidiendo el valor absoluto de 5% divido por 2
Z <- abs( qnorm((1 - .95)/2) )

# Calcular error estándar
error_est <- sigma/sqrt(n)

lim_inferior <- x_barra - ( Z * error_est )
lim_superior <- x_barra + ( Z * error_est )

# Intervalo de confianza al 95% redondeado
round(c(lim_inferior, lim_superior),2)

[1]  92.16 107.84

¿Y si ahora queremos calcular un intervalo de confianza al 99%?

# Datos relevantes
sigma = 20
n = 25
x_barra <- 100

# Calcular el valor Z para un 99%
# Estoy pidiendo el valor absoluto de 5% divido por 2
Z <- abs( qnorm((1 - .99)/2) )

# Calcular error estándar
error_est <- sigma/sqrt(n)

lim_inferior <- x_barra - ( Z * error_est )
lim_superior <- x_barra + ( Z * error_est )

# Intervalo de confianza al 95% redondeado
round(c(lim_inferior, lim_superior),2)

[1]  89.7 110.3

Calculando el intervalo de confianza de una proporción

Si en vez de un promedio de valores intervalares o de razón queremos calcular el promedio de una proporción, la formula que usamos varía solo en la forma en la que calculamos el error estándar, y podemos ver que la estructura del intervalo de confianza es igual.

\[ \text{Error Estándar de una proporción} = \sqrt \frac{p \times (1-p)}{n} \]

Con este error estándar, aplicamos una fórmula con la misma estructura para calcular un intervalo de confianza.

\[ \text{Proporción} \pm \text{Valor Z } \times \text{Error Estándar de una proporción} \]

Errores estándar en proporciones

De hecho, la formula es estrictamente la misma en ambos casos:

\[ EE = \sqrt \frac{\text{varianza}}{\text{tamaño muestral}} \]

Ocurre que en una distribución normal la varianza es un valor independiente del promedio: \(\sigma^2\), pero en una distribución binomial (que, como se explicó al inicio, modela experimentos con dos resultados), la varianza se define como dependiente del promedio: \(p \times (1 - p)\).

Veamos un ejemplo: Si se nos dice que la población de estudiantes universitarios en la PUC se obtiene una muestra de 50 personas, de las cuales 40 dicen estar interesados en vacunarse ¿Cómo se calcularía un intervalo de confianza al 95% de la proporción poblacional?

# Datos relevantes
n = 50
p_barra <- 40/50

# Calcular el valor Z para un 95%
# Estoy pidiendo el valor absoluto de 5% divido por 2
Z <- abs( qnorm((1 - .95)/2) )

# Calcular error estándar basado en p*(1-p)
error_est <- sqrt( (p_barra * (1 - p_barra)) / n )

lim_inferior <- p_barra - ( Z * error_est )
lim_superior <- p_barra + ( Z * error_est )

# Intervalo de confianza al 95% redondeado
round(c(lim_inferior, lim_superior),2)

[1] 0.69 0.91

Contrastando intervalos de confianza

Hasta ahora hemos visto como calcular un intervalo de confianza para una muestra, lo cual nos servirá de base para hacer inferencias respecto al promedio o proporción poblacional.

Sin embargo, hay otra posible aplicación de los intervalos de confianza: comparar dos muestras distintas y preguntarnos si es posible que ambas emergieran de la misma población.

Para ilustrar esta idea vamos a extender el ejemplo anterior, sin embargo, las ideas que vamos a ilustrar se aplican a todos los intervalos de confianza, no solo a los intervalos de confianza de proporciones.

Universidad 1: de la población de estudiantes universitarios en la universidad 1 se obtiene una muestra de 50 personas, de las cuales 40 dicen estar interesados en vacunarse.

Universidad 2: de la población de estudiantes universitarios en la universidad 2 se obtiene una muestra de 55 personas, de las cuales 35 dicen estar interesados en vacunarse.

¿Creemos que hay diferencias en las proporciones de personas dispuestas a vacunarse en estas dos universidades? Claramente las muestras son distintas, pero eso no es lo que nos interesa saber, en cambio, lo que realmente queremos saber es si las proporciones poblacionales son distintas.

Veamos los intervalos de confianza para cada universidad.

# Datos relevantes
n_1 = 50
p_barra_1 <- 40/50

n_2 = 55
p_barra_2 <- 35/55


# Calcular el valor Z para un 95%
# Estoy pidiendo el valor absoluto de 5% divido por 2
Z <- abs( qnorm((1 - .95)/2) )


# Calcular error estándar basado en p*(1-p)
error_est_1 <- sqrt( (p_barra_1 * (1 - p_barra_1)) / n )

lim_inferior_1 <- p_barra_1 - ( Z * error_est_1 )
lim_superior_1 <- p_barra_1 + ( Z * error_est_1 )

# Intervalo de confianza al 95% redondeado
round(c(lim_inferior_1, lim_superior_1),2)

[1] 0.69 0.91

# Calcular error estándar basado en p*(1-p)
error_est_2 <- sqrt( (p_barra_2 * (1 - p_barra_2)) / n )

lim_inferior_2 <- p_barra_2 - ( Z * error_est_2 )
lim_superior_2 <- p_barra_2 + ( Z * error_est_2 )

# Intervalo de confianza al 95% redondeado
round(c(lim_inferior_2, lim_superior_2),2)

[1] 0.50 0.77

Veamos las dos proporciones y sus intevalos de confianza gráficamente:

plot(x = c(0 , 1), y =c(.5, 2.5), type = "n"
  , axes = FALSE, ylab = "", xlab = "")

points(x = p_barra_1, y = 1)
points(x = p_barra_2, y = 2)
lines(x = c(lim_inferior_1, lim_superior_1), y = c(1,1), lwd  = 2)
lines(x = c(lim_inferior_2, lim_superior_2), y = c(2,2), lwd  = 2)
axis(1)
axis(2, at = c(1,2), labels = c("Univ 1","Univ 2"), las = 1)

Si bien las proporciones son distintas, podemos ver que los intervalos de confianza no están separados (se solapan o se tocan), lo que indica que si usamos un 95% de confianza, es posible que las muestras de ambas universidades hayan surgido de una misma población (con una proporción poblacional común).

Comparemos estas dos universidades con un una tercera universidad.

Universidad 3: de la población de estudiantes universitarios en la universidad 3 se obtiene una muestra de 45 personas, de las cuales 10 dicen estar interesados en vacunarse.

# Datos relevantes
n_3 = 45
p_barra_3 <- 10/45

# Calcular el valor Z para un 95%
# Estoy pidiendo el valor absoluto de 5% divido por 2
Z <- abs( qnorm((1 - .95)/2) )


# Calcular error estándar basado en p*(1-p)
error_est_3 <- sqrt( (p_barra_3 * (1 - p_barra_3)) / n )

lim_inferior_3 <- p_barra_3 - ( Z * error_est_3 )
lim_superior_3 <- p_barra_3 + ( Z * error_est_3 )

# Intervalo de confianza al 95% redondeado
round(c(lim_inferior_3, lim_superior_3),2)

[1] 0.11 0.34

Y comparemos gráficamente con las otras dos universidades:

plot(x = c(0 , 1), y =c(.5, 3.5), type = "n"
  , axes = FALSE, ylab = "", xlab = "")

points(x = p_barra_1, y = 1)
points(x = p_barra_2, y = 2)
points(x = p_barra_3, y = 3)
lines(x = c(lim_inferior_1, lim_superior_1), y = c(1,1), lwd  = 2)
lines(x = c(lim_inferior_2, lim_superior_2), y = c(2,2), lwd  = 2)
lines(x = c(lim_inferior_3, lim_superior_3), y = c(3,3), lwd  = 2)
axis(1)
axis(2, at = c(1,2,3), labels = c("Univ 1","Univ 2","Univ 3"), las = 1)

El intervalo de confianza de esta tercera universidad claramente está separado de los intervalos de confianza de las dos universidades anteriores. En este caso, con un 95% de confianza podemos inferir que esta tercera universidad no proviene de una misma población que las otras dos proporciones (no comparte la proporción poblacional con las dos primeras universidades).

Inferencias con una prueba Z paso a paso

Partamos con una distribución poblacional inicial, con promedio 30 y desviación estándar 4. Vamos a generar nuestra población con 100 mil casos simulando datos con el comando rnorm indicando con los argumentos la media y desviación estándar.

set.seed(2024)

poblacion.simulada <- rnorm(100000, mean = 30, sd = 4)

hist(poblacion.simulada)

round(mean(poblacion.simulada),2)

[1] 30

round(sd(poblacion.simulada),2)

[1] 4

Ahora podemos sacar muestras de nuestra población:

# Esto es una muestra de 10 casos
muestra1 <- sample(poblacion.simulada, 10)

# Este es el promedio de esta muestra
mean(muestra1)

[1] 29.17965

# Con este comando podemos calcular el promedio para una nueva muestra 
mean(sample(poblacion.simulada, 10))

[1] 30.76872

# Con este comando podemos obtener 20 muestras de nuestra población
replicate(20, mean(sample(poblacion.simulada, 10)))

 [1] 31.39672 29.53590 27.89599 28.29489 28.52923 31.29083 30.41673 31.15253
 [9] 28.64638 31.12427 29.58631 31.08422 30.09349 30.70357 29.88720 30.28285
[17] 31.54082 28.03595 31.62611 32.77431

Ahora que podemos sacar muestras, podemos crear una distribución muestral.

# Y ahora voy a generar 20 mil muestras para construir una distribución muestral
distrib.muestral <- replicate(20000, mean(sample(poblacion.simulada, 10)))

hist(distrib.muestral)

Y podemos compararlas:

# Uso mfrow para crear un lienzo con 1 fila y 2 columnas
par(mfrow=c(1,2))

# Población original
hist(poblacion.simulada, xlim = c(0,60))

# Distribución muestral con muestras de 10 casos
hist(distrib.muestral, xlim = c(0,60))

Veamos las propiedades de nuestra distribución muestral.

# Sabemos que debe tener una media igual a la población
round(mean(distrib.muestral),2)

[1] 30

# Sabemos que debe tener una desv. est. igual a sigma/raizcuadrada(n)

# Error estandár teórico
round(4/sqrt(10),2)

[1] 1.26

# Error estandár de nuestra distribución muestral
round(sd(distrib.muestral),2)

[1] 1.26

Podemos ver además que nuestra distribución muestral se comporta, tal como esperamos, como una distribución normal (media 30, EE = 1.26).

# Uso mfrow para crear un lienzo con 1 fila y 2 columnas
par(mfrow=c(1,2))

# Densidad de nuestra distribución muestral
plot(density(distrib.muestral))

# Densidad de nuestra distribución muestral y una distribución normal
plot(density(distrib.muestral))
curve(dnorm(x, mean=30, sd=1.26), add=TRUE, col="red", lwd = 4)

En clase vimos que podemos calcular el estadígrafo Z para una muestra en particular para caracterizar qué tan lejos está esa muestra del promedio de la distribución muestral. Imaginemos que tenemos una muestra de n = 10 (que corresponde a la distribución muestral que hemos estado usando) con un promedio de 33.

Veamos dánde está en nuestra distribución muestral:

plot(density(distrib.muestral), type = "n")
curve(dnorm(x, mean=30, sd=1.26), add=TRUE, col="red", lwd = 4)

arrows(33,.1, 33, 0)

Podemos calcular el estadígrafo Z para esta muestra de promedio 33:

\[ Z = \frac{\overline{x} - \mu }{E.E.} = \frac{33 - 30}{1.26} = 2.38 \]