1. El Teorema del Límite Central es uno de los más importantes en la inferencia estadística y habla sobre la convergencia de los estimadores como la proporción muestral a la distribución normal. Algunos autores afirman que esta aproximación es bastante buena a partir del umbral n>30.

a. Realice una simulación en la cual genere una población de N=1000 (Lote) y además que el porcentaje de individuos (plantas) enfermas sea del 50%.

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.8     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.0
## ✔ readr   2.1.2     ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(datos)

## Generación de población 

N <- rbinom(1000,1,0.5)

b. Genere una función que permita obtener una muestra aleatoria de la población y calcule el estimador de la proporción muestral para un tamaño de muestra dado n.

## Definición de función que calcula la media muestral 

sampling = function(n){
  
  x <- sample(N,n)
  return(mean(x))
}      

#Calcular el estimador para una muestra de n = 40 

sampling(40)
## [1] 0.5

c. Repita el escenario anterior (b) 500 veces y analice los resultados en cuanto al comportamiento de los 500 estimadores. ¿Qué tan simétricos son los datos?, ¿Son sesgados y qué pasa en cuanto a variabilidad?

salida = 1:500

for (i in 1:500){
  
  salida[[i]] <- sampling(40)  #Llenado de vector con los 500 estimadores 
}

shapiro.test(salida)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida
## W = 0.98838, p-value = 0.0005236
hist(salida, breaks = 10, main = str_c("Sampling distribution of p where p =  0.5 & n =", 40), xlab = "Sample mean")

qqnorm(salida,main = str_c("Distribución de residuos para p =  0.5 & n =", 40))
qqline(salida, col= 2)

Con respecto a la distribución del estimador de las 500 muestras de tamaño 40 podemos decir que la distribución se centra más claramente en 0.5 a diferencia de el estimador individual en el punto b que fue de 0.525. Por otro lado podemos identificar que la distribución de los 500 estimadores es aproximadamente normal (Simetrica), esta idea la podemos sustentar con el resultado de la prueba shapiro el cual arroja un p value bastante pequeño.

d. Realice los ejercicios completos b y c para tamaños de muestra n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500. Y compare los resultados de los estimadores en cuanto a la normalidad. Investigue y utilice pruebas de bondad y ajuste (shapiro wilks) y métodos gráficos (grafico qq de normalidad).

vector = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)


for (m in 1:10){

size = vector[m]  #Extrae el tamaño n de la iteración actual 

salida = 1:500  #Genera vector donde se almacenaran los 500 estimadores 

for (i in 1:500){
  
  salida[[i]] <- sampling(size)  #Llenado de vector con los 500 estimadores 
}

shapiro.test(salida)

hist(salida, breaks = 10, main = str_c("Sampling distribution of p where p =  0.5 & n =", vector[m]), xlab = "Sample mean")

qqnorm(salida,main = str_c("Distribución de residuos para p =  0.5 & n =", vector[m]))
qqline(salida, col= 2)


}

c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)

salida5 = 1:500
for (i in 1:500){
  salida5[[i]] <- sampling(5)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida5)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida5
## W = 0.92674, p-value = 6.633e-15
salida10 = 1:500
for (i in 1:500){
  salida10[[i]] <- sampling(10)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida10)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida10
## W = 0.9643, p-value = 1.15e-09
salida15 = 1:500
for (i in 1:500){
  salida15[[i]] <- sampling(15)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida15)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida15
## W = 0.97043, p-value = 1.671e-08
salida20 = 1:500
for (i in 1:500){
  salida20[[i]] <- sampling(20)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida20)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida20
## W = 0.97891, p-value = 1.24e-06
salida30 = 1:500
for (i in 1:500){
  salida30[[i]] <- sampling(30)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida30)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida30
## W = 0.98646, p-value = 0.000134
salida50 = 1:500
for (i in 1:500){
  salida50[[i]] <- sampling(50)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida50)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida50
## W = 0.99013, p-value = 0.001954
salida60 = 1:500
for (i in 1:500){
  salida60[[i]] <- sampling(60)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida60)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida60
## W = 0.99073, p-value = 0.003117
salida100 = 1:500
for (i in 1:500){
  salida100[[i]] <- sampling(100)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida100)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida100
## W = 0.99468, p-value = 0.08083
salida200 = 1:500
for (i in 1:500){
  salida200[[i]] <- sampling(200)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida200)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida200
## W = 0.99465, p-value = 0.07861
salida500 = 1:500
for (i in 1:500){
  salida500[[i]] <- sampling(500)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida500)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida500
## W = 0.99561, p-value = 0.1736

De las anteriores simulaciones podemos extraer varios findings: A mediada que el tamaño de la muestra incrementa, la normalidad en la distribución del estimador se va haciendo más evidente. Podríamos decir que desde un análisis gráfico después de n = 30 la distribución se puede considerar normal. A través de la prueba Shapiro se puede evaluar de una forma más técnica la normalidad en estas distribuciones, como era de esperarse a medida que aumenta el n, el valor-p de la prueba shapiro incrementa. Es importante aclarar que el no rechazo de la hipótesis nual significa aceptar normalidad en la distribución. A pesar de que gráficamente se podía hablar de normalidad desde n = 30, usando los p values de las pruebas se empieza a aceptar la hipótesis de normalidad a partir de n = 100 a un nivel de significancia del 5%.

e. Repita toda la simulación (puntos a – d) pero ahora con lotes con 10% y 90% de plantas enfermas. Concluya todo el ejercicio.

Simulación para lotes 10% de probabilidad
## Generación de población 

N <- rbinom(1000,1,0.1)

## Definición de función que calcula la media muestral 

sampling = function(n){
  
  x <- sample(N,n)
  return(mean(x))
}

vector = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)


for (m in 1:10){

size = vector[m]  #Extrae el tamaño n de la iteración actual 

salida = 1:500  #Genera vector donde se almacenaran los 500 estimadores 

for (i in 1:500){
  
  salida[[i]] <- sampling(size)  #Llenado de vector con los 500 estimadores 
}

shapiro.test(salida)

hist(salida, breaks = 10, main = str_c("Sampling distribution of p where p =  0.1 & n =", vector[m]), xlab = "Sample mean")

qqnorm(salida,main = str_c("Distribución de residuos para p =  0.1 & n =", vector[m]))
qqline(salida, col= 2)

}

Simulación para lotes 90% de probabilidad
## Generación de población 

N <- rbinom(1000,1,0.9)

## Definición de función que calcula la media muestral 

sampling = function(n){
  
  x <- sample(N,n)
  return(mean(x))
}

vector = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)


for (m in 1:10){

size = vector[m]  #Extrae el tamaño n de la iteración actual 

salida = 1:500  #Genera vector donde se almacenaran los 500 estimadores 

for (i in 1:500){
  
  salida[[i]] <- sampling(size)  #Llenado de vector con los 500 estimadores 
}

shapiro.test(salida)

hist(salida, breaks = 10, main = str_c("Sampling distribution of p where p =  0.9 & n =", vector[m]), xlab = "Sample mean")

qqnorm(salida,main = str_c("Distribución de residuos para p =  0.9 & n =", vector[m]))
qqline(salida, col= 2)

}

Al igual que como lo vimos en la población con p= 0.5, los resultados con las poblaciones con p = 0.1 y p=0.9 tienden a tener los mismos resultados, a medida que se aumenta el valor p la distribución del estimador tiende a la normalidad, se podría pensar que en los tamaños de muestra pequeños las distribuciones del estimador tienden a estar mas sesgadas en p=0.1 y p=0.9 que en p=0.5, pero este aparente sesgo realmente es ocacionado por la posición de la media dentro del dominio (0,1). Al estar 0.5 mas centrado en el dominio parece tener una normalidad más temprana, pero lo cierto es que la normalidad se presenta a la misma velocidad solo que las campanas están centradas en diferentes valores de p. 

2. La comparación de tratamientos es una práctica fundamental en las ciencias agropecuarias y para esto a nivel estadístico se cuenta con algunas herramientas para apoyar el proceso de toma de decisiones y lograr concluir con algún grado de confianza que los resultados observados en una muestra son representativos y se pueden asociar a los tratamientos y no se deben únicamente al azar. Por medio una simulación validemos algunos de estos resultados.

a. Suponga un escenario en el cual usted aplicó tratamientos diferentes a dos lotes y desea analizar si alguno de los dos presenta un mejor desempeño en el control de una plaga presente en ambos al momento inicial. Para ello utilizará como criterio de desempeño el tratamiento que menor % de plantas enfermas presente después de un tiempo de aplicación (es decir, si se presentan o no diferencias en las proporciones de enfermos P1 y P2). Realice una simulación en la cual genere dos poblaciones de N1=1000 (Lote1) y N2=1500 (Lote2), además asuma que el porcentaje de individuos (plantas) enfermas en ambos lotes sea la misma 10% (es decir, sin diferencias entre los tratamientos).

#Generación de lote 1 

N1 <- rbinom(1000,1,0.1)

#Generación de lote 2

N2 <- rbinom(1500,1,0.1)

b. Genere una función que permita obtener una muestra aleatoria de los lotes y calcule el estimador de la proporción muestral para cada lote (p1 y p2) para un tamaño de muestra dado n1=n2. Calcule la diferencia entre los estimadores p1-p2.

sampling2 = function(n){
  
  x <- sample(N1,n)
  y <- sample(N2,n)
  
  return(mean(x)-mean(y))
}

sampling2(50)
## [1] 0.04

c. Repita el escenario anterior (b) 500 veces y analice los resultados en cuanto al comportamiento de los 500 estimadores (diferencias p1-p2). ¿Qué tan simétricos son los datos?, ¿Son siempre cero las diferencias?

salida = 1:500

for (i in 1:500){
  
  salida[[i]] <- sampling2(40)  #Llenado de vector con los 500 estimadores 
}

shapiro.test(salida)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida
## W = 0.98572, p-value = 8.097e-05
hist(salida, breaks = 10, main = str_c("Sampling distribution of p1-p2 where p1=p2=0.1 & n1=n2 =", 40), xlab = "Sample mean")

De esta simulación podemos extraer que la distribución de diferentes es simétrica y que estas no son siempre cero a pesar de que su distribución si esta centada en el 0.

d. Realice los puntos b y c para tamaños de muestra n1=n2=5, 10, 15, 20, 30, 50, 60, 100, 200, 500. Y compare los resultados de los estimadores (p1-p2) en cuanto a la normalidad. También analice el comportamiento de las diferencias y evalúe. ¿Considera que es más probable concluir que existen diferencias entre los tratamientos con muestras grandes que pequeñas, es decir, cuál considera usted que es el efecto del tamaño de muestra en el caso de la comparación de proporciones?

#Generación de lote 1 

N1 <- rbinom(1000,1,0.1)

#Generación de lote 2

N2 <- rbinom(1500,1,0.1)

sampling2 = function(n){
  
  x <- sample(N1,n)
  y <- sample(N2,n)
  
  return(mean(x)-mean(y))
}


vector = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)


for (m in 1:10){

size = vector[m]  #Extrae el tamaño n de la iteración actual 

salida = 1:500  #Genera vector donde se almacenaran los 500 estimadores 

for (i in 1:500){
  
  salida[[i]] <- sampling2(size)  #Llenado de vector con los 500 estimadores 
}

shapiro.test(salida)

hist(salida, breaks = 10, main = str_c("Sampling distribution of p1-p2 where p1=p2=0.1 & n1=n2 =", vector[m]), xlab = "Sample mean")

}

salida5 = 1:500
for (i in 1:500){
  salida5[[i]] <- sampling2(5)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida5)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida5
## W = 0.8959, p-value < 2.2e-16
salida10 = 1:500
for (i in 1:500){
  salida10[[i]] <- sampling2(10)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida10)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida10
## W = 0.94812, p-value = 3.13e-12
salida15 = 1:500
for (i in 1:500){
  salida15[[i]] <- sampling2(15)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida15)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida15
## W = 0.96432, p-value = 1.158e-09
salida20 = 1:500
for (i in 1:500){
  salida20[[i]] <- sampling2(20)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida20)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida20
## W = 0.97062, p-value = 1.827e-08
salida30 = 1:500
for (i in 1:500){
  salida30[[i]] <- sampling2(30)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida30)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida30
## W = 0.97941, p-value = 1.644e-06
salida50 = 1:500
for (i in 1:500){
  salida50[[i]] <- sampling2(50)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida50)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida50
## W = 0.98743, p-value = 0.0002649
salida60 = 1:500
for (i in 1:500){
  salida60[[i]] <- sampling2(60)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida60)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida60
## W = 0.98982, p-value = 0.001541
salida100 = 1:500
for (i in 1:500){
  salida100[[i]] <- sampling2(100)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida100)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida100
## W = 0.99134, p-value = 0.005068
salida200 = 1:500
for (i in 1:500){
  salida200[[i]] <- sampling2(200)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida200)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida200
## W = 0.99402, p-value = 0.04631
salida500 = 1:500
for (i in 1:500){
  salida500[[i]] <- sampling2(500)  #Llenado de vector con los 500 estimadores 
}
shapiro.test(salida500)
## 
##  Shapiro-Wilk normality test
## 
## data:  salida500
## W = 0.99333, p-value = 0.02612

No, no se considera que sea más probable concluir que existen diferencias entre los tratamientos con muestras grandes que pequeñas, por el contrario, la centralidad de la distribución de las diferencias en 0 se hace más fuerte a medida que el tamaño de la muestra se incrementa.

e. Ahora realice nuevamente los puntos a-d bajo un escenario con dos lotes, pero de proporciones de enfermos diferentes (P1=0.1 y P2=0.15), es decir, el tratamiento del lote 1 si presentó un mejor desempeño reduciendo en un 5% el porcentaje de enfermos. Bajo este nuevo escenario compare la distribución de estas diferencias (p1- p2) con las observadas bajo igualdad de condiciones en los lotes. ¿Qué puede concluir? ¿Existen puntos en los cuales es posible que se observen diferencias de p1- p2 bajo ambos escenarios (escenario 1: sin diferencias entre P1 y P2, escenario 2: diferenciade 5%)?

#Generación de lote 1 

N1 <- rbinom(1000,1,0.1)

#Generación de lote 2

N2 <- rbinom(1500,1,0.15)

sampling2 = function(n){
  
  x <- sample(N1,n)
  y <- sample(N2,n)
  
  return(mean(x)-mean(y))
}

vector = c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)


for (m in 1:10){

size = vector[m]  #Extrae el tamaño n de la iteración actual 

salida = 1:500  #Genera vector donde se almacenaran los 500 estimadores 

for (i in 1:500){
  
  salida[[i]] <- sampling2(size)  #Llenado de vector con los 500 estimadores 
}

shapiro.test(salida)

hist(salida, breaks = 10, main = str_c("Sampling distribution of p1-p2 where p1=0.1 & p2=0.15 & n1=n2 =", vector[m]), xlab = "Sample mean")

}

Al igual que en los puntos anteriores la diferencia real de las proporciones se hace más evidente a medida que el tamaño de la muestra se hace más grande. Incluso podemos evidenciar que se puede concluir que no hay diferencas en las proporciones a valores bajos de n. La centralización de la distribución de diferencias más abajo de cero se empieza a notar a partir de n = 20, por lo que para valores menores se puede llegar a concluir que no hay diferencia significativa en las poblaciones igual que cuando realmente no las hay.

3. Con base a los artículos “Statistical Errors: P values, the gold standard of statistical validity, are not as reliable as many scientists assume” & “Statisticians issue warning on P values: Statement aims to halt missteps in the quest for certainty” escriba un resumen (máximo 2 páginas) sobre ambos artículos e incluya en este sus opiniones en cuanto al uso del valor p como criterio de decisión en inferencia estadística.

A lo largo de la historia del último siglo, el valor p se ha convertido en una herramienta estadística de mucho uso por parte de los investigadores, este fue establecido con un objetivo que con poca frecuencia se le da lugar. En este ensayo analizaremos al valor p desde su concepto, su aplicación en la actualidad y aterrizaremos las ideas en el taller desarrollado previamente. Desarrollaremos este escrito teniendo como marco las simulaciones realizadas en este RMarkdown y los ensayos Statisticians issue warning on P values y STATISTICAL ERRORS P values, the ‘gold standard’ of statistical validity, are not as reliable as many scientists assume.

El valor p se define como la probabilidad de tener un evento especifico o más extremo asumiendo la veracidad de una hipótesis nula previamente ya establecido. Cuando el valor p es menor a un nivel de significancia previamente establecido (normalmente 0.05) ósea que la probabilidad de que suceda tal evento es uno más extremo sea menor a este valor, se rechaza la hipótesis nula y se concluye que el resultado es significativo. Análogamente, cuando el valor p es mayor a este nivel de significancia no se rechaza la hipótesis nula, por lo que no hay evidencia suficiente para negar la validez de esta.

Si un estimador resulta ser significativo, y la prueba apoya la hipótesis planteada, no significa estrictamente que esta hipótesis sea cierta, simplemente que es poco probable que se presente un resultado como el muestreado o uno mas extremo asumiendo que la hipótesis nula es cierta. Y es justamente en este punto donde aparece uno de los errores mas comunes con respecto al uso del valor p. Durante el último siglo muchos investigadores han optado por usar la significancia de un valor p como elemento decisivo para sustentar una hipótesis y asegurar su validez. Esto por su puesto ha generado múltiples inconsistencia y falsos positivos en diferentes campos tal y como se expusieron en los artículos del marco de desarrollo de este ensayo. Un ejemplo claro de esto es el estudio desarrollado por Matt Motyl en donde se concluyó que los extremistas suelen ver el mundo a blanco y negro mientras que los moderados ven más sombras grises, el resultado de su estudio arrojó un valor p de 0.01 pero al usar más datos el p-value cambió a 0.59.

Analicemos uno de los puntos desarrollados en este taller, específicamente el punto 1.d; en este ejercicio se solicitó desarrollar una simulación en la cual se analizara la normalidad de 500 estimadores de la proporción de una muestra que proviene de una población con parámetro p = 0.5 en distintos niveles del tamaño muestral, los valores de n utilizados fueron 10, 15, 20, 30, 50, 60, 100, 200 y 500. El análisis de normalidad se hiso utilizando la prueba Shapiro y los valores p arrojados. El análisis de normalidad se hiso utilizando la prueba Shapiro y podemos evidenciar que incluso en valores más altos de 30 la prueba rechaza la normalidad de los estimadores, aun cuando sabemos que si por la ley de los grandes números estos si son.

Podríamos definir entonces, o como lo decía Fisher, el valor P se debe usar como guía numérica aproximada de la fuerza que una evidencia tendría contra la hipótesis nula. Pero no debe ser un elemento decisivo al momento de definir si una hipótesis es verdadera o no. Ahora, si al ejecutar múltiples pruebas con diferentes muestras, se encuentran valores de p consecutivamente bajos se podría decir que los efectos observados son poco probables que se debieran solo por variabilidad natural.

Podríamos concluir que un valor significativo de p no significa que algo sea cierto, sino que se debe tener en cuenta o se debe considerar para un estudio más profundo, se debería tener como punto de partida para estudiar el fenómeno de interés a fondo y así poder apoyar la hipótesis con diferentes evidencias estadísticas.