Cuando el histograma de una lista de números se aproxima a la distribución normal, se puede utilizar una fórmula matemática para aproximar la proporción de valores o resultados en cualquier intervalo dado:
\[ \mbox{Pr}(a < x < b) = \int_a^b \frac{1}{\sqrt{2\pi\sigma^2}} \exp{\left( \frac{-(x-\mu)^2}{2 \sigma^2} \right)} \, dx \]
\(\mu\) es la media y \(\sigma\) es desviación estandar \[ \mu = \frac{1}{M} \sum\limits_{i=1}^M x_{i} \] \[ \sigma ^{2} = \frac{1}{M} \sum\limits_{i=1}^M \big( x_{i}- \mu \big)^{2} \]
Si X esta normalmente distribuido, entonces es facil calcular una proporcion en un intervalo con i individuos. Z se deduce como distribucion normal estandar, y se dice cuantas deviaciones estandar de distancia esta de la media. Esto es definido como puntuación Z.
La media y desviación estandar pueden describir la propocion de cualquier intervalo, un resumen de los datos. Solo si los datos se distribuyen de manera normal.
Que diferencia se nota cuando aumentas el tamaño de la mustra en una población, como el ejemplo siguiente.
poblacion<-unlist(read.csv("femaleControlsPopulation.csv"))
n <- 1000
null5<- vector("numeric",n)
null50<-vector("numeric",n)
for(i in 1:n){
x5<-sample(poblacion,5)
x50<-sample(poblacion,50)
null5[i] <- mean(x5)
null50[i]<-mean(x50)
}
par(mfrow=c(1,2))
hist(null5,density = 20)
hist(null50,density = 20)
Podemos calcular la proporción especifica en cada intervalo de valores.
La proporción de la muestra de tamaño de 5 (\(\mu\)= 23.89 y \(\sigma\)= 1.55), entre 23 y 25. Con esta ecuación mean(null5 < 25 & null5 > 23)
es 0.524. Este contiene menos de la mita de tos los valores.
Y para la muestra de tamaño de 50 (\(\mu\)= 23.9 y \(\sigma\)= 0.41), entre el mismo intervalo mean(null50 < 25 & null50 > 23)
es 0.982. El intervalo contiene serca de la totalidad de los datos.
Ahora con el mismo intervalo (23-25) en una distribución normal, cual es la proporción con una media de 23.9 y una desviación de 0.43: para esto ocupamos la función pnorm()
como pnorm((25-23.9)/0.43)-pnorm((23-23.9)/0.43)
es pnorm()
como 0.9766.
La inferencia estadística es un enfoque muy potente, que permite responder a una pregunta, con solo mirar una parte de la población, muestra aleatoria. El primer paso en la inferencia estadística es entender la población.
\[Media Población -> \mu _{X} = \frac{1}{m} \sum\limits_{i=1}^m x_{i} \Longleftrightarrow \mu _{Y} = \frac{1}{m} \sum\limits_{i=1}^m y_{i}\] \[Observaciones -> X_{1} ,..., X_{n} \Longleftrightarrow Y_{1} ,..., Y_{n} \] \[Media Muestra -> \overline{X} = \frac{1}{M} \sum\limits_{i=1}^M X_{i} \Longleftrightarrow \overline{Y} = \frac{1}{M} \sum\limits_{i=1}^M Y_{i}\]
Lo que queremos saber es: Que tan cerca esta la media muestreal de la media poblacional. Con ayuda del teorema del límite central podemos contestar.
En el ejemplo del peso del ratón, tenemos dos poblaciones:. Ratones hembra y machos en dietas de control y en dietas altas en grasa, el peso es el resultado de interés. En este ejemplo tenemos toda la población.
library(rafalib)
datos<-read.csv("mice_pheno.csv")
expMice<-na.omit(datos)
head(expMice)
## Sex Diet Bodyweight
## 1 F hf 31.94
## 2 F hf 32.48
## 3 F hf 22.82
## 4 F hf 19.92
## 5 F hf 32.22
## 6 F hf 27.50
machoCont<-expMice$Bodyweight[expMice$Diet=="chow" & expMice$Sex=="M"]
machoGras<-expMice$Bodyweight[expMice$Diet=="hf" & expMice$Sex=="M"]
hembraCont<-expMice$Bodyweight[expMice$Diet=="chow" & expMice$Sex=="F"]
hembraGras<-expMice$Bodyweight[expMice$Diet=="hf" & expMice$Sex=="F"]
par(mfrow=c(2,2))
hist(machoCont,density = 20)
hist(machoGras,density = 20)
hist(hembraCont,density = 20)
hist(hembraGras,density = 20)
Ahora vamos a tomar dos muestra de la población, de un tamaño de 25 ratones.
set.seed(1)
Xm<-sample(machoCont,25)
set.seed(1)
Ym<-sample(machoGras,25)
set.seed(1)
Xh<-sample(hembraCont,25)
set.seed(1)
Yh<-sample(hembraGras,25)
Muestra aleatoria de la poblacion de Machos control y dieta, la media de X y Y son 32.1 y 34.77 respectivamente.
abs((mean(machoCont)-mean(machoGras))-(mean(Ym)-mean(Xm)))
el resultado es 1.212.Muestra aleatoria de la poblacion de Hembras control y dieta, la media de X y Y son 23.17 y 26.28 respectivamente.
abs((mean(hembraCont)-mean(hembraGras))-(mean(Yh)-mean(Xh)))
el resultado es 0.736Estos nos dice que las hembras de las muestras estimadas estan mas sercanas a la diferencia de la población que los machos. Esto se debe a la pequeña varianza de las hembras 3.42 en comparación con los machos 4.42.
La media de la muestra sigue una distribución normal, con la media de la población y la desviación estandar, entonces podemos saber la probabilidad de cualquier intervalo. La desviación estandar es como varia la población con respecto a la media, estre más grande mayor varia y a la inversa.
El TLC es uno de los resultados matemáticos de uso más frecuente en la ciencia. Nos dice que cuando el tamaño de la muestra es grande, la media de \(\overline{Y}\) de una muestra aleatoria sigue una distribución normal centrada en la media de la población \(\mu _{Y}\) y la desviación estandar es igual a la desviación estandar de la población \(\sigma _{Y}\) dividido por la raíz cuadrada del tamaño de la muestra N:
\[ \frac{\bar{Y} - \mu}{\sigma_Y/\sqrt{N}} \]
se aproxima a una distribución normal centrada en 0 y con una desviación estándar 1.
Ahora nos interesa la diferencia de medias entre dos muestras. Debido a que tanto la muestra y poblacion son normales, la diferencia es normal, así, y la varianza (la desviación estándar al cuadrado) es la suma de las dos varianzas. Bajo la hipótesis nula de que no hay diferencia entre las medias de población.
\[ \frac{\bar{Y}-\bar{X}}{\sqrt{\frac{\sigma_X^2}{M} + \frac{\sigma_Y^2}{N}}} \sim N \big(0,1\big) \]
Es aproximada por una distribución normal centrada en 0 (no hay diferencia) y con una desviación estándar \(\sqrt{\sigma_X^2 +\sigma_Y^2}/\sqrt{N}\)
Se aproxima por una distribución normal centrada en 0 y desviación estándar 1. Utilizando esta aproximación hace que el calculo de p-valores sea simple, porque conocemos la proporción de distribución de cualquier valor.
De una distribución normal, para calcular la proporción de datos en 1, 2 o 3 desviaciones estandar se utiliza la funcion pnorm()
.
x<-rnorm(10000)
hist(x,xlab=NULL,ylab=NULL,main=NULL,axes=F)
par(new=T)
plot(density(x),col=2,lwd=3)
Para una desviación estandar pnorm(1)-pnorm(-1)
igual 0.683.
Para dos desviaciones estandar pnorm(2)-pnorm(-2)
igual 0.954.
Para tres desviaciones estandar pnorm(3)-pnorm(-3)
igual 0.997.
Del los datos de expMice
, calcular la proporción de machos control en una desviación estandar, para caulcular la proporcion, calculamos \(Z_{i}=\frac{X_{i}-\overline{X}}{\sigma_{X}}\), luego calculamos la media de z menor igual a 1.
z<-(machoCont-mean(machoCont))/popsd(machoCont)
mean(abs(z)<=1)
## [1] 0.6950673
Con la funcción qqplot()
podemos comparar la distribuccion de nuestros datos con la distribución normal.
qqnorm(z)
abline(0,1,lwd=3,col=2)
Entre mas cerca de la linea (teorica) esten los puntos, los datos se distribuye de manera normal.
Vemos todos las graficas de las cuatro poblaciones hembras y machos, en dieta y control.
par(mfrow=c(2,2))
qqnorm(z,main = "Macho-Control");abline(0,1,lwd=3,col=2)
z1<-(machoGras-mean(machoGras))/popsd(machoGras)
qqnorm(z1,main = "Macho-AltoGrasa");abline(0,1,lwd=3,col=2)
z2<-(hembraCont-mean(hembraCont))/popsd(hembraCont)
qqnorm(z2,main = "Hembra-Control");abline(0,1,lwd=3,col=2)
z3<-(hembraGras-mean(hembraGras))/popsd(hembraGras)
qqnorm(z3,main = "Hembra-AltoGrasa");abline(0,1,lwd=3,col=2)
Si la distribución de una lista de números es aproximadamente normal, entonces si elegimos un número al azar de esta distribución, se siguen una distribución normal. Sin embargo, es importante recordar que aunque indica que algunos datos tiene una distribución no implica necesariamente que esta al azar. Para esto ocupamos la función replicate()
.
MediasMachoCont <- replicate(10000, mean( sample(machoCont, 25)))
par(mfrow=c(1,2))
hist(MediasMachoCont)
qqnorm(MediasMachoCont)
qqline(MediasMachoCont,lwd=3,col=2)
Contamos con los datos de toda la población, y vamos a utilizar los datos del experimento, para observar que tambien el TLC se aproxima a los parametros de la población.
muCont<-mean(hembraCont)
muGras<-mean(hembraGras)
sdCont<-popsd(hembraCont)
sdGras<-popsd(hembraGras)
Creamos una función y un bucle para hacer los gráficos.
Ns <- c(3,12,25,50)
B <- 10000
res <- sapply(Ns,function(n) {
replicate(B,mean(sample(hembraGras,n))-mean(sample(hembraCont,n)))
})
par(mfrow=c(2,2))
for (i in seq(along=Ns)) {
titleavg <- signif(mean(res[,i]),3)
titlesd <- signif(popsd(res[,i]),3)
title <- paste0("N=",Ns[i]," Avg=",titleavg," SD=",titlesd)
qqnorm(res[,i],main=title)
qqline(res[,i],col=2)
}
Las cuatro aproximaciones se parecen a una normal, Debido a que la población misma está relativamente cerca de una distribución normal, los promedios son cerca de lo normal tambien (la suma de las normales es también una normal). > Cuanto más se desvian los puntos de la linea, peor es la aproximación a una normal.
En la práctica, en realidad calculamos una relación: dividimos por la desviación estándar estimada. Aquí es donde el tamaño de la muestra comienza a importar más.
Ns <- c(3,12,25,50)
B <- 10000
computetstat <- function(n) {
y <- sample(hembraGras,n)
x <- sample(hembraCont,n)
(mean(y)-mean(x))/sqrt(var(y)/n+var(x)/n)
}
res <- sapply(Ns,function(n) {
replicate(B,computetstat(n))
})
mypar(2,2)
for (i in seq(along=Ns)) {
qqnorm(res[,i],main=Ns[i])
qqline(res[,i],col=2)
}
Así vemos que para N = 3, la CLT no proporciona una aproximación útil. Para N = 12, hay una ligera desviación en los valores más altos, aunque la aproximación parece útil. Para el 25 y el 50, la aproximación es perfecto. Esta simulación sólo demuestra que N = 12 es lo suficientemente grande en este caso.