penguins <- read.csv("penguins_size.csv")

DISTRIBUCIONES DE PROBABILIDAD

Teniendo en cuenta la variable aleatoria de que se trate, podemos obtener probabilidades puntuales, acumuladas, cuantiles de la variable de acuerdo a su distribución, etc. Se deben tener en cuenta los parámetros de cada distribución y las características de las variables en estudio. Además, podemos aplicar los conceptos adquiridos en la creación de gráficos para representar las funciones asociadas a las variables aleatorias en estudio.

El paquete stats de R (que se instala por defecto al instalar R, y se carga en memoria siempre que iniciamos sesión) implementa numerosas funciones para la realización de cálculos asociados a distintas distribuciones de probabilidad. Entre las utilizadas más comunmente podemos citar:

Distribuciones Discretas Distribuciones Continuas

Distribución	Nombre en R	Distribución	Nombre en R
Binomial	binom	Uniforme	unif
Poisson	pois	Normal	norm
-	-	t Student	t
-	-	F Fisher	F
-	-	Chi-Cuadrado	chisq

Para cualquier consulta utilizar help(Distributions) Que nos muestra el listado de distribuciones de probabilidad disponibles en el paquete stats. Otras distribuciones están disponibles en otros paquetes.

Para cada distribución, R dispone de cuatro funciones. Se puede acceder a cada una de ellas simplemente precediendo el nombre de la distribución que figura en la tabla anterior por la letra que se indica a continuación:

d: función de densidad o de probabilidad. p: función de distribución q: función para el cálculo de cuantiles. r: función para simular datos con dicha distribución.

Así, por ejemplo, para la distribución normal, la función de densidad se obtiene como dnorm(), la función de distribución como pnorm(), los cuantiles se calculan mediante qnorm() y se pueden generar valores aleatorios con distribución normal mediante rnorm(). Puede consultarse la ayuda, help(dnorm) para conocer la sintaxis específica de estas funciones.

Ejemplo: La variable aleatoria X tiene distribución binomial con n=7 y probabilidad de éxito de 0.8. Calcular:

la probabilidad de X=3

n=7
p=0.8
x=3
a = pbinom(x,n,p)
print("La probabilidad de que X=3 es:")

[1] "La probabilidad de que X=3 es:"

[1] 0.033344

la probabilidad de que X sea a lo sumo 5

x=5
b = dbinom(x,n,p)
print("La probabilidad de que X sea a lo sumo 5 es:"); b

[1] "La probabilidad de que X sea a lo sumo 5 es:"

[1] 0.2752512

la probabilidad de que X sea al menos 3

x = 3
c = 1-dbinom(x-1,n,p)
print("La probabilidad de que X sea al menos 3 es:"); c

[1] "La probabilidad de que X sea al menos 3 es:"

[1] 0.9956992

el valor de X para un probabilidad acumulada de 0.5

q=0.5
d = qbinom (q,n,p)
print("Existe una probabilidad de 0,5 de que X sea menor o igual a "); d

[1] "Existe una probabilidad de 0,5 de que X sea menor o igual a "

[1] 6

Generar una muestra aleatoria de tamaño 30 para la variable X~B(n,p)

#set.seed(15)
c=30
m1 = rbinom(c,n,p)
m1

 [1] 6 5 6 5 5 5 5 5 7 6 7 5 5 7 6 5 6 6 5 7 7 4 7 6 5 5 7 5 5 5

Podemos representar fácilmente la función de probabilidad de la distribución binomial:

plot(dbinom(0:7,7,0.8),type="h",xlab="x",ylab="p(x)", sub="Función de Probabilidad X~B(7,0.8)")

También podemos representar su función de distribución:

plot(stepfun(0:7,pbinom(0:8,7,0.8)),xlab="x",ylab="F(x)",sub="Función de distribución X~B(7,0.8)",main="")

# Distribución de la media en el muestreo. Teorema central del límite

R es una herramienta excelente para ayudar a entender el significado de muchos conceptos estadísticos. Un concepto a menudo mal entendido es el de la distribución de un estadístico en el muestreo. Para aclarar este concepto podemos comenzar por realizar un sencillo experimento: tomar una muestra aleatoria simple en una población normal de, por ejemplo, media μ=170 y desviación típica σ=12. Fijamos inicialmente un tamaño de muestra n=46, tomamos la muestra y calculamos su media:

n=46
muestra1=rnorm(n,170,12)
media1=mean(muestra1)
media1

[1] 169.8113

Si repetimos el proceso obtenemos medias distintas:

muestra2=rnorm(n,170,12)
media2=mean(muestra2)
media2

[1] 173.1973

muestra3=rnorm(n,170,12)
media3=mean(muestra3)
media3

[1] 170.7145

Como a priori es imposible predecir en cada muestreo cuál será el valor medio resultante, la media muestral es una variable aleatoria. Las medias muestrales se parecen a la media de la población, μ=170. Cabe preguntarse entonces: ¿tenderá la media muestral a comportarse de esta manera en todas las muestras? ¿Cuánto llega a apartarse de la media de la población?

Para responder a estas preguntas podemos repetir el proceso anterior no tres, sino muchísimas más veces. Una forma muy eficiente de replicar muchas veces el proceso anterior consiste en encapsular el proceso de muestreo en una función que haremos depender del tamaño de la muestra:

mediaMuestral=function(n){
  muestra=rnorm(n,170,12)
  media=mean(muestra)
  return(media)
}

Cada vez que ejecutemos esta función estaremos eligiendo una muestra de tamaño n de esa población N(170,12) y calculando su media:

mediaMuestral(46)

[1] 167.5578

mediaMuestral(46)

[1] 170.9634

mediaMuestral(46)

[1] 166.6454

Para repetir m veces el proceso de extraer una muestra de tamaño n y calcular su media podemos utilizar la función replicate():

m=10000
muchasMedias=replicate(m,mediaMuestral(46))

Las primeras 20 medias obtenidas en esta simulación son:

muchasMedias[1:20]

 [1] 169.4820 170.0007 171.4020 171.3043 171.8336 169.0188 170.0623 168.7319
 [9] 167.4776 168.9937 173.1806 171.4193 168.0661 171.0253 168.7706 167.1751
[17] 171.0986 168.5304 171.2944 171.0604

La media y desviación típica de todas estas medias muestrales son:

mean(muchasMedias)

[1] 169.9965

sd(muchasMedias)

[1] 1.756664

Por último, representamos gráficamente la distribución de frecuencias de estas medias muestrales mediante un histograma, al que le superponemos una densidad normal:

hist(muchasMedias,xlab="Media muestral", ylab="Frecuencia", col="lightcyan",
     xlim=c(160,180),freq=FALSE,ylim=c(0,0.75),
     main="Histograma de las medias muestrales observadas\n en 10000 muestras de tamaño 46")
curve(dnorm(x,170,sd(muchasMedias)),xlim=c(160,180),col="blue",lwd=2,add=TRUE)

A estas alturas, queda claro que la media muestral es una variable aleatoria, y que su distribución de probabilidad es muy parecida a la normal; además la media de todas las medias muestrales es casi idéntica a la media de la población.

¿Cuál es el efecto de aumentar el tamaño de la muestra? Podemos evaluarlo repitiendo el proceso anterior para n=50, n=100 y n=500 (en todos los casos presentaremos los gráficos en la misma escala para facilitar la comparación):

muchasMedias50=replicate(m,mediaMuestral(50))
muchasMedias100=replicate(m,mediaMuestral(100))
muchasMedias500=replicate(m,mediaMuestral(500))
mean(muchasMedias50); sd(muchasMedias50)

[1] 170.0075

[1] 1.719387

mean(muchasMedias100); sd(muchasMedias100)

[1] 170.0222

[1] 1.19973

mean(muchasMedias500); sd(muchasMedias500)

[1] 170.0052

[1] 0.5406369

hist(muchasMedias50,xlab="Media muestral", ylab="Frecuencia", col="lightcyan",
     xlim=c(160,180),freq=FALSE,ylim=c(0,0.75),
     main="Histograma de las medias muestrales observadas\n en 10000 muestras de tamaño 50")
curve(dnorm(x,170,sd(muchasMedias50)),xlim=c(160,180),col="blue",lwd=2,add=TRUE)

hist(muchasMedias100,xlab="Media muestral", ylab="Frecuencia", col="lightcyan",
     xlim=c(160,180),freq=FALSE,ylim=c(0,0.75),
     main="Histograma de las medias muestrales observadas\n en 10000 muestras de tamaño 100")
curve(dnorm(x,170,sd(muchasMedias100)),xlim=c(160,180),col="blue",lwd=2,add=TRUE)

hist(muchasMedias500,xlab="Media muestral", ylab="Frecuencia",
     col="lightcyan",
     xlim=c(160,180),freq=FALSE,ylim=c(0,0.75),
     main="Histograma de las medias muestrales observadas\n en 10000 muestras de tamaño 500")
curve(dnorm(x,170,sd(muchasMedias500)),xlim=c(160,180),col="blue",lwd=2,add=TRUE,n=200)

Resulta evidente en las gráficas, que el incremento del tamaño de la muestra tiene como consecuencia que las posibles medias muestrales se concentran más en torno a su media; en otras palabras, cuanto más grande sea la muestra, más probable es que la media de la muestra esté muy cerca de la media de la población. Puede comprobarse además que la desviación típica de las medias muestrales es en todos los casos un valor muy parecido a σ/√n=12/√n.

INFERENCIA ESTADÍSTICA

Tanto para estimaciones por intervalos como para pruebas de hipótesis, podemos crear un algoritmo para realizar estos cálculos rápidamente. El procedimiento es semejante para ambas técnicas, ya que el estadígrafo que se utiliza es el mismo.

Supuestos de la prueba T de Student

Para esta prueba es necesario considerar el nivel de medición de las variables, para lo cual se requiere de una variable cuantitativa y una variable cualitativa de escala nominal (dos grupos).

Las condiciones para calcular intervalos de confianza o aplicar un test basados en la distribución T-Student son:

Independencia: Las observaciones tienen que ser independientes unas de las otras. Para ello el muestreo tiene que ser aleatoria y el tamaño de la muestra inferior al 10% de la población.
Normalidad: Las poblaciones que se comparan tienn que distribuirse de forma normal. En caso de cierta asimetría, los T-test son considerablemente robustos cuando el tamaño de muestra es mayor o igual a 30.
Homogeneidad de varianzas: La varianza de ambas poblaciones comparadas debe ser igual. En caso de no cumplirse con esta condición se puede usar una prueba Welch Two Sample t-test. Esta corrección se incorpora a través de los grados de libertad permitiendo compensar la diferencia de varianzas (con el inconveniente de que se pierde precisión).

Para trabajar tomamos de la base /penguins/ la variable longitud de culmen (mm).

Prueba unilateral

1° Paso: Planteo de hipótesis \(H_c\):“El largo de culmen en los pingüinos adelia es mayor a 38 mm” \(H_0: \mu\ = 38mm\) \(H_1: \mu\ > 38mm\)

2° Paso: Tenemos que plantear el nivel de significancia, Por defecto el valor de \(\alpha\ \) es de 0,05. 3° Paso: elegir el estadígrafo de prueba. Dado que a tamaño de muestra grande T-Student se acerca a una distribución normal con media cero y varianza uno, trabajaremos con éste estadígrafo. establecer la regla de decisión y el cálculo del estadígrafo. Luego tomar la decisión de aceptar o rechazar la \(H_0\) o no.

4° Paso: Regla de decisión:

Para determinar el valor crítico necesitamos el nivel de significancia y/o probabilidad acumulada además de los grados de libertad

adelia = subset(penguins,species=="Adelie")
alpha = 0.05
gl = length(adelia$culmen_length_mm)-1
t_c=qt(alpha,gl,lower.tail = FALSE)
t_c

[1] 1.655007

Si no se plantea el nivel de significancia se realiza con un \(\alpha\ = 0.05\) por defecto.

Para realizar el cálculo deestadígrafo de prueba:

t.test(adelia$culmen_length_mm ,mu=38,alternative="greater")


    One Sample t-test

data:  adelia$culmen_length_mm
t = 3.6513, df = 150, p-value = 0.0001799
alternative hypothesis: true mean is greater than 38
95 percent confidence interval:
 38.43266      Inf
sample estimates:
mean of x 
 38.79139

Si es lateral izquierda, el argumento *alternative# es “less”. Por defecto es bilateral.

Prueba bilateral

\(H_c\):“El largo de culmen en los pingüinos adelia es diferente a 38 mm” \(H_0: \mu\ = 38mm\) \(H_1: \mu\ \neq 38mm\)

# Por test
alpha=0.05

t_conf=qt(alpha/2,length(adelia$culmen_length_mm)-1,lower.tail = FALSE)
t_conf

[1] 1.975799

t.test(adelia$culmen_length_mm,mu=38)


    One Sample t-test

data:  adelia$culmen_length_mm
t = 3.6513, df = 150, p-value = 0.0003598
alternative hypothesis: true mean is not equal to 38
95 percent confidence interval:
 38.36312 39.21966
sample estimates:
mean of x 
 38.79139

El intervalo de confianza sale con el test bilateral o puede realizarse de forma manual.

# De manera manual
Media=mean(adelia$culmen_length_mm)
DE=sd(adelia$culmen_length_mm)
alpha=0.05
LI=Media-qt(alpha/2,length(adelia$culmen_length_mm)-1,lower.tail = FALSE)*(DE/sqrt(length(adelia$culmen_length_mm)))
LS=Media+qt(alpha/2,length(adelia$culmen_length_mm)-1,lower.tail = FALSE)*(DE/sqrt(length(adelia$culmen_length_mm)))                                                       
LI ; LS

[1] NA

[1] NA

Diferencia de medias de dos poblaciones normales independientes

Si nos hacemos la pregunta:¿Existen diferencias en el largo de culmen entre los pingüinos adelia de la isla Brisco y Dream?

El primer paso es realizar una exploración de los datos y verificar los supuestos de la prueba; la prueba t requiere que los datos tengan una distribución normal y con varianzas iguales.

Gráfico de histogramas para visualizar la distribución de los datos entre ambas islas:

library(tidyverse)

── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
✔ ggplot2 3.3.6      ✔ purrr   0.3.4 
✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.2.1      ✔ stringr 1.4.1 
✔ readr   2.1.2      ✔ forcats 0.5.2 
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(car)

Loading required package: carData

Attaching package: 'car'

The following object is masked from 'package:dplyr':

    recode

The following object is masked from 'package:purrr':

    some

library(rstatix)


Attaching package: 'rstatix'

The following object is masked from 'package:stats':

    filter

adelia_islas <- adelia %>%
  filter(species == "Adelie", island %in% c("Dream", "Biscoe"))

ggplot(adelia_islas, aes(x = culmen_length_mm, fill = island, col = island))+
  geom_histogram(alpha = 0.3,bins = 9,binwidth = 2)+
  labs(y = "frecuencia absoluta",
       x = "longitud del culmen (mm)",
       title="Historgrama", 
       fill = "isla", 
       col = "isla")+
  facet_grid(island ~.)

Tambien es posible realizar un QQ-Plot (quantil-quantil) para evaluar visualmente si la distribución de los datos de ambas islas se asemejan a una distribución normal

ggplot(adelia_islas, aes(sample = culmen_length_mm, col = island))+
  stat_qq()+
  stat_qq_line()+
  facet_grid(.~ island)

De acuerdo a los gráficos, es evidente que ambas poblaciones tienen distribuciones muy parecida a la normal. Para validar estos resutlados, realizaremos una prueba de Shapiro_wilks para evaluar la significancia de la normalidad.

RECORDAR

La hipótesis nula es la que está bajo estudio y coincide conla hipótesis científica. \(H_0\) : “La variable aleatoria longitud de culmen de los pingüinos adelia de la isla Dream sigue una distribución normal” \(H_1\) : “La variable aleatoria longitud de culmen de los pingüinos adelia de la isla Dream sigue otra distribución”

adelia_islas_shap<- adelia_islas %>%
  group_by(island) %>%
  shapiro_test(culmen_length_mm) %>%
  print()

# A tibble: 2 × 4
  island variable         statistic     p
  <chr>  <chr>                <dbl> <dbl>
1 Biscoe culmen_length_mm     0.977 0.520
2 Dream  culmen_length_mm     0.984 0.640

Tanto el histograma como el qqplot dan soporte a que los datos tienen una distribución aproximadmente normal. Para poder darle validez, se realiza la prueba de shapiro, la cual resulta que las muestras aportan suficiente evidencia con una significancia del 5%, de que la longitud de culmen (mm) de los pingüinos adelia de ambas islas tienen una distribución normal

El siguiente paso es evaluar el supuesto de homogeneidad de varianzas. El histograma nos da un indicio visual que la dispersión de los datos es similar, pero es posible evaluar la significancia con una prueba de Levene implementada en el paquete car:

La prueba de levene plantea en la hipótesis nulas que las varianzas son iguales

leveneTest(culmen_length_mm ~ island, data = adelia_islas)

Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  1  0.0054 0.9415
      98

La prueba nos dice que las varianzas son homogeneas entre los grupos

El siguiente paso para la exploración de los datos es obtener el promedio de los estadígrafos de la longitud del culmen:

adelia_islas_res = adelia_islas %>%
  group_by(island) %>%
  summarise(media = mean(culmen_length_mm),
            DE = sd(culmen_length_mm),
            N = n())
adelia_islas_res

# A tibble: 2 × 4
  island media    DE     N
  <chr>  <dbl> <dbl> <int>
1 Biscoe  39.0  2.48    44
2 Dream   38.5  2.47    56

De manera alternativa, este procedimiento se puede realizar con ayuda del paquete rstatix

adelia_islas_res2 = adelia_islas %>%
  group_by(island) %>%
  get_summary_stats(culmen_length_mm)
adelia_islas_res2

# A tibble: 2 × 14
  island variable       n   min   max median    q1    q3   iqr   mad  mean    sd
  <chr>  <chr>      <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Biscoe culmen_le…    44  34.5  45.6   38.7  37.7  40.7  3.02  2.74  39.0  2.48
2 Dream  culmen_le…    56  32.1  44.1   38.6  36.8  40.4  3.62  2.74  38.5  2.46
# … with 2 more variables: se <dbl>, ci <dbl>

Visualizamos mediante un boxplot los datos:

ggplot(adelia_islas, aes(x = island, y = culmen_length_mm, fill = island))+
  geom_boxplot(outlier.shape = "")+
  geom_point(position = position_jitter(0.1), col = "grey50")+
  theme_light()+
  labs(x = "Isla",
       y = "longitud del culmen (mm)",
       title = "Longitud del culmen por isla en el pingüino Adelie",
       fill = "Isla", col = "Isla")

Finalmente, hacemos la prueba t de student. Esta puede realizarse siguiendo el comando base en R t.test()

t.test(culmen_length_mm ~ island, data = adelia_islas, var.equal = TRUE )


    Two Sample t-test

data:  culmen_length_mm by island
t = 0.95016, df = 98, p-value = 0.3444
alternative hypothesis: true difference in means between group Biscoe and group Dream is not equal to 0
95 percent confidence interval:
 -0.5151265  1.4615551
sample estimates:
mean in group Biscoe  mean in group Dream 
            38.97500             38.50179

En esta prueba, la hipótesis nula es que la media de las dos poblaciones es igual, lo que implica que la hipótesis alternativa es que la diferencia de las medias no es igual a cero.

En la prueba de hiótesis, se desea aceptar o rechazar la hipotesis nula con ciertos intervalos de confianza. Dado que probamos la diferencias entre ambas medias, el intervalo de confianza especifica el intervalo de los valores de esta diferencia.

Alternativametne podemos usar la prueba implementada en rstatixs:

adelia_islas_ttest <- adelia_islas %>%
  t_test(culmen_length_mm ~ island, var.equal = TRUE) %>%
  print()

# A tibble: 1 × 8
  .y.              group1 group2    n1    n2 statistic    df     p
* <chr>            <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl>
1 culmen_length_mm Biscoe Dream     44    56     0.950    98 0.344

Comparación de medias poblacionales dependientes (pareadas)

Dos medias son dependientes o pareadas cuando proceden de grupos o muestras dependientes, esto es, cuando existe una relación entre las observaciones de las muestras. Este escenario ocurre a menudo cuando los resultados se generan a partir de los mismos individuos bajo condiciones distintas.

Las pruebas pareadas tienen la ventaja frente a los independientesde que se puede controlar mejor la variación no sistemática, ya que se bloquean al estar examinando los mismos individuos dos veces, no dos grupos de individuos distintos.

Se realizan los mismos análisis descriptivo e inferencial sobre los supuestos que vimos anteriormente. Lo que debemos cambiar es la sintaxis de la prueba

# t_test(variable ~ grupo, paired = TRUE, var.equal = TRUE)

Otras pruebas las veremos más adelante.

Probabilidad e Inferencia

E. Martinez y E. Maure

2022-09-26