Selección de Datos

library("dplyr")
data<-read.csv("femaleMiceWeights.csv")
control<-filter(data,Diet=="chow") %>% select(Bodyweight) %>% unlist
# control<-data$Bodyweight[data$Diet=="chow"] otra forma sin usar paquete
tratamiento<-filter(data,Diet=="hf") %>% select(Bodyweight) %>% unlist
# control<-data$Bodyweight[data$Diet=="hf"] otra forma sin usar paquete
head(control);head(tratamiento)
## Bodyweight1 Bodyweight2 Bodyweight3 Bodyweight4 Bodyweight5 Bodyweight6 
##       21.51       28.14       24.04       23.45       23.68       19.79
## Bodyweight1 Bodyweight2 Bodyweight3 Bodyweight4 Bodyweight5 Bodyweight6 
##       25.71       26.37       22.80       25.34       24.97       28.14

La media para control es 23.81gramos y para tramiento 26.83gramos, La diferencia entre control y la media es 3.02gramos, tiene un mayor peso la dieta con un alto contenido de grasa.

Esta diferencia es realmente debido a la dieta o es debido a otra casa que la casualidad de las dietas. Puede ser que seleccionaramos los ratones mas pesados o los mas ligeros, porque algunos de datos de la dieta control estan por encima de la media, 3 datos (28.14, 28.4, 26.91)

Variables Aleatorias

Para comprender sobre las variables aleatorias, vamos a ocupar toda la población de los datos, esto rara vez lo tenemos para analizar.

poblacion<-unlist(read.csv("femaleControlsPopulation.csv"))
mean(sample(poblacion,12))
## [1] 24.06083
mean(sample(poblacion,12))
## [1] 25.02833
mean(sample(poblacion,12))
## [1] 25.405

Cada vez que que hacemos un muestreo de la poblacion con la función sample(), obtenemos valores diferentes y por lo tanto medias diferentes. Vemos como fluctua una variable aleatoria, porque tenemos acceso a todos los datos de la población, y podemos simular los muestreos.

Así la diferencia de los 3g entre tratamiento y control, podria deberse al azar, para comprobarlo se tiene que probar estadisticamente.

Seleccionamos muestras, en donde este garantizado que no tenga efecto en la dieta, tomamos una muestra aleatoria de la población y comparamos sus medias.

obs<-mean(tratamiento)-mean(control)
controlP <- sample(poblacion,12)
tretamientP <- sample(poblacion,12)
mean(tretamientP) - mean(controlP)
## [1] 0.5033333

Podemos repetir una y otra vez, obteniendo distintas diferencias, estos multiples valores de medias para la Ho, se le llama distribución nula. Son todas las realizaciones posibles en virtud de la Ho.

Si conocemos la distribución nula, podemos saber la proporcion de valores para cualquier intervalo.
Para registrar todas diferencias y obtener la distribución nula creamos un bucle para registrar 10000 valores.

n<-10000
null<-vector("numeric",n)
for (i in 1:n){
      controlP <- sample(poblacion,12)
      tretamientP <- sample(poblacion,12)
      null[i]<-mean(tretamientP) - mean(controlP)
}
hist(null,density = 20)

p-Valor

Podemos saber cual es el valor mas alto de esta distribución nula 5.1716667 que es mucho mayor que 3.0208333. Con el histograma podemos tener una idea que tan probable es ver valores mayores de 3, que es la diferencia entre tratamiento y control.
Podemos calcular la proporción de veces que la distribución nula es mas grande a la observada.

mean(null>=obs)
## [1] 0.0121
mean(abs(null)>=obs)
## [1] 0.0248

Esto nos dice que la distribución Ho es mas grande que media la observada (3.02) en 1.21% veces, en valores absolutos es 2.48%. A esto se le llama p-valor que responde a la pregunta: Cual es la probabilidad de un resultado de la Distribución Nula, es mas grande que lo que observamos cuando la Ho es verdadera.

Al cambiar el número de interacciones o muestras, no hay mucha diferencia de los p-valores, pero aumentado o disminuyendo el tamaño de la muestra esto causa una diferencia mayor y cambia el p-valor.

n<-10000
null5<-vector("numeric",n)
for (i in 1:n){
      sam50 <- sample(poblacion,5)
      null5[i]<-mean(sam50)-mean(poblacion)
}
hist(null5,density = 20)

mean(null5>=obs)
## [1] 0.0294

En este ejemplo al disminuir el tamaño de muestra a 5 ratones por muestra, el p-valor cambia de 0.0248 a 0.0461, o sea, que la distribución Ho es mas grande que la media observada en 4.61% veces, aumento casi el doble.

Distribución

La forma más sencilla de pensar en una distribución es como una descripción compacta de muchos números. Por ejemplo, supongamos que ha medido la altura de todos los hombres en una población.

library(UsingR)
x <- father.son$fheight
round(sample(x,30),1)
##  [1] 66.9 60.8 62.4 70.2 66.2 68.7 64.6 68.5 66.3 67.1 72.5 72.2 68.8 66.4
## [15] 68.8 67.8 68.8 71.2 68.6 67.6 72.4 68.1 72.0 69.3 70.6 67.4 67.4 64.9
## [29] 65.6 70.3

Para definir una una distribución calculamos, para todos los valores posibles de a, la proporción de números en nuestra lista inferiores a a. Usamos la siguiente notación: \[ F(a) \equiv \mbox{Pr}(x \leq a) \]

Esto se llama la función de distribución acumulativa (FDC), permite calcular la proporción en cualquier punto. Cuando el FDC se deriva de los datos, también llamamos la FDC empírica (FDCE). Podemos trazar F(a) frente a una de esta manera:

smallest <- floor( min(x) )
largest <- ceiling( max(x) )
valores <- seq(smallest, largest,len=300)
heightecdf <- ecdf(x)
plot(valores, heightecdf(valores),type="l",xlab="a (Height in inches)",ylab="Pr(x <= a)")

Además, el ECDF en realidad no es tan popular como histogramas, que nos dan la misma información, pero nos muestran que la proporción de los valores en intervalos: \[ \mbox{Pr}(a \leq x \leq b) = F(b) - F(a) \]

hist(x,density = 20,xlab="Alturas en Pulgada",main="Alturas de Hombres Adultos")

Mostrando este histograma es mucho más informativo que mostrar números. Con este sencillo argumento, podemos aproximar el número de individuos en cualquier intervalo dado. Por ejemplo, hay cerca de 70 personas de más de seis pies (72 pulgadas) de altura.

Ejemplo

  • Que proporción de países en 1952, tienen una esperanza de vida menor igual a 40 años.
library(gapminder)
data(gapminder)
head(gapminder)
##       country continent year lifeExp      pop gdpPercap
## 1 Afghanistan      Asia 1952  28.801  8425333  779.4453
## 2 Afghanistan      Asia 1957  30.332  9240934  820.8530
## 3 Afghanistan      Asia 1962  31.997 10267083  853.1007
## 4 Afghanistan      Asia 1967  34.020 11537966  836.1971
## 5 Afghanistan      Asia 1972  36.088 13079460  739.9811
## 6 Afghanistan      Asia 1977  38.438 14880372  786.1134
ano1952<-gapminder[gapminder$year==1952,]
mean(ano1952$lifeExp<=40)
## [1] 0.2887324
  • Que proporción de países en 1952, tienen una esperanza de vida entre 40 y 60 años.
mean(ano1952$lifeExp<=60) - mean(ano1952$lifeExp<=40)
## [1] 0.4647887