Prueba T para una media en R
Comenzamos la serie de artículos con una entrada sobre una de las técnicas estadísticas más utilizadas en las ciencias sociales y de la salud: la prueba T de Student. Existen algunas variaciones de la prueba que se aplican en distintas situaciones, y veremos cómo podemos aplicarlas todas con R.
Empezaremos por la prueba T para una sola media, que se utiliza para evaluar si una media obtenida en una muestra es lo suficientemente distinta de un valor determinado que llamamos \(\mu\). Se puede utilizar ésta técnica para comparar los valores que hemos obtenido con un punto de corte determinado.
Imaginemos que hacemos una prueba para un examen de acceso a la universidad en un centro educativo, y desde el gobierno se da un premio a los centros cuya nota media es mayor que 8. Si vemos que los alumnos del centro sacan una media de 8.2, ¿Cómo de seguros podemos estar de que el día del examen volverán a obtener una media superior a 8?
Ignorando factores como el estrés o que el examen sea más difícil, éste es el tipo de preguntas que se suelen responder mediante ésta técnica. Se llama prueba T para una media porque únicamente tenemos la media de un grupo. En cambio, en la prueba T tradicional tenemos dos medias, y la pregunta que tratamos de responder es si ambas medias son iguales o distintas. Eso lo veremos en otra ocasión.
Veamos cómo aplicar ésta técnica con R. El código es el siguiente:
No obstante, no podemos aplicar la función sin más, debemos ver si se cumplen los supuestos. Para entenderlo mejor vamos a plantear el siguiente problema. R tiene una serie de funciones que nos permiten generar números al azar que siguen una distribución determinada. Una de las más utilizadas es la función rnorm(). Que nos permite generar números que siguen la distribución normal.
Por ejemplo rnorm(100) me genera un vector que contiene 100 números que se distribuyen normalmente. Además podemos utilizar otros parámetros: rnorm(n = 50, mean = 100, sd = 15) genera 50 números que siguen una distribución normal con media 100 y desviación típica 15.
Aunque le indiquemos que queremos que nuestros datos tengan media 100, en cada “muestra”, el valor de la media será distinto:
set.seed(123) #Fijamos la semilla aleatoria
n_muestras <- 1000 #Numero de muestras
medias <- rep(NA, n_muestras)
for(i in 1:n_muestras){
medias[i] = mean(rnorm(n = 30, mean = 100, sd = 15))
}
hist(medias, breaks = 20)
Vemos que hemos obtenido valores bastante dispares para las medias. El mínimo es 92.11, y el máximo ha sido 108.98. ¿Cómo podemos saber si éstos valores son razonablemente parecidos al valor teórico que debería salir (100)?
Podemos utilizar ésta técnica para responder a la pregunta. Nos centramos en una “muestra” concreta. Pero antes debemos comprobar el cumplimiento de los supuestos. Para ésta prueba únicamente hay uno, el de normalidad.
muestra <- rnorm(n = 30, mean = 100, sd = 15)
#La prueba T requiere que los datos sean normales, comprobamos si se cumple el supuesto
shapiro.test(muestra)##
## Shapiro-Wilk normality test
##
## data: muestra
## W = 0.95856, p-value = 0.2845El p-valor es mayor que 0.05, por lo que mantenemos la hipótesis nula (los datos se distribuyen normalmente). En éste caso es fácil que los datos se distribuyan normalmente, ya que los hemos generado con una función que utiliza ésta distribución para generar los números (ojo usando ésta función podemos obtener muestras que no son normales, especialmente con valores pequeños de n).
Casi todas las técnicas estadísticas tienen algún supuesto. Antes de aplicarlas debemos comprobar si se cumplen
Si el p-valor fuese menor de 0.05 (o el alfa que seleccionemos), entonces no podemos concluir que los datos son normales y, dado que no se cumple el supuesto, deberíamos aplicar la alternativa no paramétrica: el test de Wilcoxon (ver la función desarrollada al final para ver cómo se aplica el test de Wilcoxon en R).
Como vemos que los datos cumplen los supuestos, podemos aplicar la prueba T. Recuperamos el comando del principio.
##
## One Sample t-test
##
## data: muestra
## t = 0.71503, df = 29, p-value = 0.4803
## alternative hypothesis: true mean is not equal to 100
## 95 percent confidence interval:
## 97.13691 105.94112
## sample estimates:
## mean of x
## 101.539El p-valor es mayor que 0.05, por lo que mantenemos la hipótesis nula (H0). Por defecto, la función utiliza significación bilateral, lo que implica que la hipótesis nula afirma que el valor de la media poblacional (real) es igual a \(\mu\), mientras que la alternativa (H1) dice que el valor de la media en nuestra muestra es distinto (mayor o menor) que \(\mu\), en éste caso fijado en 100. Al mantener H0, concluímos que el valor obtenido no es significativamente distinto del valor de \(\mu\).
Podemos ver cómo van variando los p-valores en función del valor del parámetro \(\mu\). Para nuestra muestra, la media vale 101.54.
valores_mu <- seq(90, 110, 0.1)
p_valores <- rep(NA, length(valores_mu))
for(i in 1:length(valores_mu)){
p_valores[i] = t.test(muestra, mu = valores_mu[i])$p.value
}
plot(valores_mu, p_valores,
pch = 19, col = "violet")
abline(h = 0.05, col = "red", lwd = 2)
Observamos que el p-valor adopta el valor máximo (1) cuando especificamos un valor de \(\mu\) cercano al valor obtenido en la muestra, y va decreciendo a ambos lados a medida que nos alejamos del mismo. La línea roja representa el punto 0.05, lo que significa que todos los valores de \(\mu\) cuyo p-valor está por debajo de la misma son poco compatibles con los datos de nuestra muestra.
Dicho de otro modo, si generamos otra muestra es poco probable (pero no imposible) que su media sea 90 ó 110. Los puntos en los que la línea roja corta a la morada son los límites del intervalo de confianza. El intervalo de confianza comprende aquellos valores que son más compatibles con la hipótesis nula. Si en nuestra muestra obtenemos un valor que no está incluído en el intervalo, debemos rechazar H0.
Para responder a la pregunta que planteábamos sobre los centros educativos deberíamos utilizar la significación unilateral derecha (queremos ver si el valor es mayor que 8). Para ello, podemos usar el parámetro alternative de la función t.test(). Acepta los valores "two.sided" (por defecto) para utilizar la significación bilateral (H0 es que el valor de la media es igual al dado), "greater" para la unilateral izquierda (H0 está formulada en términos de menor o igual) o "less" para la unilateral derecha (H0 está formulada en términos de mayor o igual).
Podemos ver el aspecto que tiene el gráfico anterior al cambiar el tipo de significación:
temp <- round(mean(muestra), 0)
valores_mu <- seq(temp - 10, temp + 10, 0.05)
p_valores_bilateral <- rep(NA, length(valores_mu))
p_valores_derecha <- rep(NA, length(valores_mu))
p_valores_izquierda <- rep(NA, length(valores_mu))
for(i in 1:length(valores_mu)){
p_valores_bilateral[i] = t.test(muestra, mu = valores_mu[i])$p.value
p_valores_izquierda[i] = t.test(muestra, mu = valores_mu[i],
alternative = "greater")$p.value
p_valores_derecha[i] = t.test(muestra, mu = valores_mu[i],
alternative = "less")$p.value
}
plot(valores_mu, p_valores_bilateral,
pch = 19, col = "violet", main = "Sig. bilateral y unilaterales",
xlab = "Valores de mu", ylab = "P valores de los contrastes")
points(valores_mu, p_valores_derecha,
pch = 19, col = "blue")
points(valores_mu, p_valores_izquierda,
pch = 19, col = "red")
abline(h = 0.05, col = "red", lwd = 2)
En violeta tenemos los p-valores con significación bilateral, en azul con la significación derecha y en rojo con la significación izquierda. Observamos que si combinásemos la derecha y la izquierda (cosa que no debemos hacer en ningún caso), el “intervalo de confianza” resultante es más estrecho que el bilateral. Al utilizar la significación unilateral los contrastes tienen mayor potencia.
Vemos que para ambos tipos de significación unilateral la curva es sigmoidea.
