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.
  1. 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%.
poblacion = c(rep("enfermo",500),rep("sano" , 500))
  1. 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.
calc_p_gorro=function(n){
muestra=sample(poblacion,size = n)
p_gorro=sum(muestra=="enfermo")/n
return(p_gorro)
}
calc_p_gorro(n = 100)
## [1] 0.51
  1. 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?
require(ggplot2)
posibles_p_gorros = sapply(rep(100,500),calc_p_gorro)

resumen = summary(posibles_p_gorros)
Coeficiente_asimetria_Pearson_CAP = 3*(mean(posibles_p_gorros)-median(posibles_p_gorros))/sd(posibles_p_gorros)
Varianza= var(posibles_p_gorros)
Desviacion = sd(posibles_p_gorros)
Coeficiente_variación_CV = sd(posibles_p_gorros)/mean(posibles_p_gorros)*100

resultados = data.frame(Coeficiente_asimetria_Pearson_CAP, Varianza, Desviacion, Coeficiente_variación_CV)
resumen
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.370   0.460   0.500   0.499   0.530   0.680
resultados
Coeficiente_asimetria_Pearson_CAP Varianza Desviacion Coeficiente_variación_CV
-0.064476 0.0023416 0.0483901 9.698192 
hist(posibles_p_gorros)
abline(v=median(posibles_p_gorros), col="blue", lwd=2)
abline(v=mean(posibles_p_gorros), col="red", lwd=2)
legend(x ="topright",legend = c("media", "mediana"), col = c("red", "blue"),lwd = 4 )

#--------------------------------------ANÁLISIS---------------------------------------------

# Para una muestra de n=100 y 500 repeticiones, se observa que la media y la mediana son muy similares, esto es un síntoma de simetría que se puede observar en el histograma, en donde la media se representa por la linea de color rojo y la mediana con color azul; el coeficiente de asimetría de Pearson CAP indica que la distribución es simétrica ya que CAP ≈ 0. En el histograma se observa que no existe sesgo (o es muy leve dado que el CAP es muy cercano a cero), la media de la muestra es muy similar a la media real , los datos se distribuyen con una alta simetría. A pesar que los datos no tienen sesgo, presentan un coeficiente de variación del 9% lo que significa que la media aritmética es representativa del conjunto de datos. 

#-------------------------------------------------------------------------------------------
  1. 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).
require(ggplot2)
require(ggpubr)

posibles_p_gorros_5 = sapply(rep(5,500),calc_p_gorro)
posibles_p_gorros_10 = sapply(rep(10,500),calc_p_gorro)
posibles_p_gorros_15 = sapply(rep(15,500),calc_p_gorro)
posibles_p_gorros_20 = sapply(rep(20,500),calc_p_gorro)
posibles_p_gorros_30 = sapply(rep(30,500),calc_p_gorro)
posibles_p_gorros_50 = sapply(rep(50,500),calc_p_gorro)
posibles_p_gorros_60 = sapply(rep(60,500),calc_p_gorro)
posibles_p_gorros_100 = sapply(rep(100,500),calc_p_gorro)
posibles_p_gorros_200 = sapply(rep(200,500),calc_p_gorro)
posibles_p_gorros_500 = sapply(rep(500,500),calc_p_gorro)

p_value_5 =shapiro.test(posibles_p_gorros_5)$p.value
p_value_10 =shapiro.test(posibles_p_gorros_10)$p.value
p_value_15 =shapiro.test(posibles_p_gorros_15)$p.value
p_value_20 =shapiro.test(posibles_p_gorros_20)$p.value
p_value_30 =shapiro.test(posibles_p_gorros_30)$p.value
p_value_50 =shapiro.test(posibles_p_gorros_50)$p.value
p_value_60 =shapiro.test(posibles_p_gorros_60)$p.value
p_value_100 =shapiro.test(posibles_p_gorros_100)$p.value
p_value_200 =shapiro.test(posibles_p_gorros_200)$p.value
p_value_500 =shapiro.test(posibles_p_gorros_500)$p.value

vector_p_value = c(p_value_5, p_value_10, p_value_15, p_value_20, p_value_30, p_value_50, p_value_60,p_value_100,p_value_200,p_value_500 )
tamaño_muestra= c(5,10,15,20,30,50,60,100,200,500)
datos_p_value=data.frame(vector_p_value,tamaño_muestra)
g1 = ggplot(datos_p_value, aes(y=vector_p_value, x=tamaño_muestra)) +  geom_line()+ geom_point()
g2 =ggqqplot(posibles_p_gorros_5)
g3 =ggqqplot(posibles_p_gorros_500)
ggarrange(g1,g2, g3,labels = c("A", "B(N=5)", "C(N=500)"),ncol = 3, nrow = 1, widths = c(2, 1.5,1.5))

plot(density(posibles_p_gorros_5),xlim=c(0, 1), ylim=c(0, 26), lwd=3)
lines(density(posibles_p_gorros_10), lwd=3, col='blue')
lines(density(posibles_p_gorros_15), lwd=3, col='red')
lines(density(posibles_p_gorros_20), lwd=3, col='green')
lines(density(posibles_p_gorros_30), lwd=3, col='purple')
lines(density(posibles_p_gorros_50), lwd=3, col='orange')
lines(density(posibles_p_gorros_60), lwd=3, col='yellow')
lines(density(posibles_p_gorros_100), lwd=3, col='gray')
lines(density(posibles_p_gorros_200), lwd=3, col='cyan')
lines(density(posibles_p_gorros_500), lwd=3, col='magenta')
legend(x = "topright",legend = c("N=5","N=10","N=15","N=20","N=30","N=50","N=60","N=100","N=200","N=500"), col = c("black","blue", "red","green", "purple", "orange", "yellow","gray","cyan","magenta"),lwd = 4 )

#--------------------------------------ANÁLISIS---------------------------------------------

#En cuanto a normalidad se puede observar que el p_value incrementa a medida que incrementa el tamaño de muestra como se observa en el gráfico A. Esto quiere decir que entre más pequeño es el p_value, más se aleja de la hipótesis nula (P=0.5). En los gráficos QQ para normalidad (B y C) se compara para muestras de n=5 y n=500 respectivamente, se observa que a medida que aumenta el tamaño de muestras los puntos se distribuyen cerca a la linea, lo que es un síntoma de normalidad. Finalmente, en el diagrama de densidad se evidencia cómo a medida que aumenta el tamaño de muestra, los datos se distribuyen en forma de campana de Gauss.

#-------------------------------------------------------------------------------------------
  1. Repita toda la simulación (puntos a – d) pero ahora con lotes con 10% y 90% de plantas enfermas. Concluya todo el ejercicio.
#-------------------------------------------------------------------------------------------
#Población 10% enfermos, 90% sanos
require(ggplot2)
poblacion_10_90 = c(rep("enfermo",100),rep("sano" , 900))
calc_p_gorro_10_90=function(n){
muestra=sample(poblacion_10_90,size = n)
p_gorro=sum(muestra=="enfermo")/n
return(p_gorro)
}
posibles_p_gorros_10_90 = sapply(rep(100,500),calc_p_gorro_10_90)

resumen = summary(posibles_p_gorros_10_90)
Coeficiente_asimetria_Pearson_CAP = 3*(mean(posibles_p_gorros_10_90)-median(posibles_p_gorros_10_90))/sd(posibles_p_gorros_10_90)
Varianza= var(posibles_p_gorros_10_90)
Desviacion = sd(posibles_p_gorros_10_90)
Coeficiente_variación_CV = sd(posibles_p_gorros_10_90)/mean(posibles_p_gorros_10_90)*100
resultados = data.frame(Coeficiente_asimetria_Pearson_CAP, Varianza, Desviacion, Coeficiente_variación_CV)
resumen
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0200  0.0800  0.1000  0.1006  0.1200  0.1900
resultados
Coeficiente_asimetria_Pearson_CAP Varianza Desviacion Coeficiente_variación_CV
0.0665058 0.0007822 0.0279675 27.79515 
hist(posibles_p_gorros_10_90)
abline(v=median(posibles_p_gorros_10_90), col="blue", lwd=2)
abline(v=mean(posibles_p_gorros_10_90), col="red", lwd=2)
legend(x ="topright",legend = c("media", "mediana"), col = c("red", "blue"),lwd = 4 )

#----------------------------

require(ggpubr)
posibles_p_gorros_5 = sapply(rep(5,500),calc_p_gorro_10_90)
posibles_p_gorros_10 = sapply(rep(10,500),calc_p_gorro_10_90)
posibles_p_gorros_15 = sapply(rep(15,500),calc_p_gorro_10_90)
posibles_p_gorros_20 = sapply(rep(20,500),calc_p_gorro_10_90)
posibles_p_gorros_30 = sapply(rep(30,500),calc_p_gorro_10_90)
posibles_p_gorros_50 = sapply(rep(50,500),calc_p_gorro_10_90)
posibles_p_gorros_60 = sapply(rep(60,500),calc_p_gorro_10_90)
posibles_p_gorros_100 = sapply(rep(100,500),calc_p_gorro_10_90)
posibles_p_gorros_200 = sapply(rep(200,500),calc_p_gorro_10_90)
posibles_p_gorros_500 = sapply(rep(500,500),calc_p_gorro_10_90)

p_value_5 =shapiro.test(posibles_p_gorros_5)$p.value
p_value_10 =shapiro.test(posibles_p_gorros_10)$p.value
p_value_15 =shapiro.test(posibles_p_gorros_15)$p.value
p_value_20 =shapiro.test(posibles_p_gorros_20)$p.value
p_value_30 =shapiro.test(posibles_p_gorros_30)$p.value
p_value_50 =shapiro.test(posibles_p_gorros_50)$p.value
p_value_60 =shapiro.test(posibles_p_gorros_60)$p.value
p_value_100 =shapiro.test(posibles_p_gorros_100)$p.value
p_value_200 =shapiro.test(posibles_p_gorros_200)$p.value
p_value_500 =shapiro.test(posibles_p_gorros_500)$p.value

vector_p_value = c(p_value_5, p_value_10, p_value_15, p_value_20, p_value_30, p_value_50, p_value_60,p_value_100,p_value_200,p_value_500 )
tamaño_muestra= c(5,10,15,20,30,50,60,100,200,500)
datos_p_value=data.frame(vector_p_value,tamaño_muestra)
g4 = ggplot(datos_p_value, aes(y=vector_p_value, x=tamaño_muestra)) +  geom_line()+ geom_point()
g5 =ggqqplot(posibles_p_gorros_5)
g6 =ggqqplot(posibles_p_gorros_500)
ggarrange(g4,g5, g6,labels = c("D", "E(N=5)", "F(N=500)"),ncol = 3, nrow = 1, widths = c(2, 1.5,1.5))

plot(density(posibles_p_gorros_5),xlim=c(0, 0.5), ylim=c(0, 35))
lines(density(posibles_p_gorros_10), lwd=3, col='blue')
lines(density(posibles_p_gorros_15), lwd=3, col='red')
lines(density(posibles_p_gorros_20), lwd=3, col='green')
lines(density(posibles_p_gorros_30), lwd=3, col='purple')
lines(density(posibles_p_gorros_50), lwd=3, col='orange')
lines(density(posibles_p_gorros_60), lwd=3, col='yellow')
lines(density(posibles_p_gorros_100), lwd=3, col='gray')
lines(density(posibles_p_gorros_200), lwd=3, col='cyan')
lines(density(posibles_p_gorros_500), lwd=3, col='magenta')
legend(x ="topright",legend = c("N=5","N=10","N=15","N=20","N=30","N=50","N=60","N=100","N=200","N=500"), col = c("black","blue", "red","green", "purple", "orange", "yellow","gray","cyan","magenta"),lwd = 4 )

#Población 90% enfermos, 10% sanos
require(ggplot2)
poblacion_90_10 = c(rep("enfermo",900),rep("sano" , 100))
calc_p_gorro_90_10=function(n){
muestra=sample(poblacion_90_10,size = n)
p_gorro=sum(muestra=="enfermo")/n
return(p_gorro)
}

posibles_p_gorros_90_10 = sapply(rep(100,500),calc_p_gorro_90_10)

resumen = summary(posibles_p_gorros_90_10)
Coeficiente_asimetria_Pearson_CAP = 3*(mean(posibles_p_gorros_90_10)-median(posibles_p_gorros_90_10))/sd(posibles_p_gorros_90_10)
Varianza= var(posibles_p_gorros_90_10)
Desviacion = sd(posibles_p_gorros_90_10)
Coeficiente_variación_CV = sd(posibles_p_gorros_90_10)/mean(posibles_p_gorros_90_10)*100
resultados = data.frame(Coeficiente_asimetria_Pearson_CAP, Varianza, Desviacion, Coeficiente_variación_CV)
resumen
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.8200  0.8800  0.9000  0.8992  0.9200  0.9800
resultados
Coeficiente_asimetria_Pearson_CAP Varianza Desviacion Coeficiente_variación_CV
-0.0834032 0.0007872 0.0280565 3.12009 
hist(posibles_p_gorros_90_10)
abline(v=mean(posibles_p_gorros_90_10), col="red", lwd=4)
abline(v=median(posibles_p_gorros_90_10), col="blue", lwd=4)
legend(x = "topright",legend = c("media", "mediana"), col = c("red", "blue"),lwd = 4 )

#---------------------------- 
require(ggpubr)

posibles_p_gorros_5 = sapply(rep(5,500),calc_p_gorro_90_10)
posibles_p_gorros_10 = sapply(rep(10,500),calc_p_gorro_90_10)
posibles_p_gorros_15 = sapply(rep(15,500),calc_p_gorro_90_10)
posibles_p_gorros_20 = sapply(rep(20,500),calc_p_gorro_90_10)
posibles_p_gorros_30 = sapply(rep(30,500),calc_p_gorro_90_10)
posibles_p_gorros_50 = sapply(rep(50,500),calc_p_gorro_90_10)
posibles_p_gorros_60 = sapply(rep(60,500),calc_p_gorro_90_10)
posibles_p_gorros_100 = sapply(rep(100,500),calc_p_gorro_90_10)
posibles_p_gorros_200 = sapply(rep(200,500),calc_p_gorro_90_10)
posibles_p_gorros_500 = sapply(rep(500,500),calc_p_gorro_90_10)

p_value_5 =shapiro.test(posibles_p_gorros_5)$p.value
p_value_10 =shapiro.test(posibles_p_gorros_10)$p.value
p_value_15 =shapiro.test(posibles_p_gorros_15)$p.value
p_value_20 =shapiro.test(posibles_p_gorros_20)$p.value
p_value_30 =shapiro.test(posibles_p_gorros_30)$p.value
p_value_50 =shapiro.test(posibles_p_gorros_50)$p.value
p_value_60 =shapiro.test(posibles_p_gorros_60)$p.value
p_value_100 =shapiro.test(posibles_p_gorros_100)$p.value
p_value_200 =shapiro.test(posibles_p_gorros_200)$p.value
p_value_500 =shapiro.test(posibles_p_gorros_500)$p.value

vector_p_value = c(p_value_5, p_value_10, p_value_15, p_value_20, p_value_30, p_value_50, p_value_60,p_value_100,p_value_200,p_value_500 )
tamaño_muestra= c(5,10,15,20,30,50,60,100,200,500)
datos_p_value=data.frame(vector_p_value,tamaño_muestra)
g7 = ggplot(datos_p_value, aes(y=vector_p_value, x=tamaño_muestra)) +  geom_line()+ geom_point()
g8 =ggqqplot(posibles_p_gorros_5)
g9 =ggqqplot(posibles_p_gorros_500)
ggarrange(g7,g8, g9,labels = c("G", "H(N=5)", "I(N=500)"),ncol = 3, nrow = 1, widths = c(2, 1.5,1.5))

plot(density(posibles_p_gorros_5),xlim=c(0.5,1), ylim=c(0, 35))
lines(density(posibles_p_gorros_10), lwd=3, col='blue')
lines(density(posibles_p_gorros_15), lwd=3, col='red')
lines(density(posibles_p_gorros_20), lwd=3, col='green')
lines(density(posibles_p_gorros_30), lwd=3, col='purple')
lines(density(posibles_p_gorros_50), lwd=3, col='orange')
lines(density(posibles_p_gorros_60), lwd=3, col='yellow')
lines(density(posibles_p_gorros_100), lwd=3, col='gray')
lines(density(posibles_p_gorros_200), lwd=3, col='cyan')
lines(density(posibles_p_gorros_500), lwd=3, col='magenta')
legend(x = "topleft",legend = c("N=5","N=10","N=15","N=20","N=30","N=50","N=60","N=100","N=200","N=500"), col = c("black","blue", "red","green", "purple", "orange", "yellow","gray","cyan","magenta"),lwd = 4 )

Análisis General

ggarrange(g1,g4, g7,labels = c("50-50", "10-90", "90-10"),ncol = 3, nrow = 1, widths = c(1.5, 1.5,1.5))

#--------------------------------------ANÁLISIS---------------------------------------------

#En cuanto a normalidad se puede observar que el p_value incrementa a medida que incrementa el tamaño de muestra como se observa en los gráficos  anteriores para cada proporción. Sin embargo, la proporción 50-50 es la única que evidencia normalidad a partir de tamaños de muestra N>200, mientras que el p_value en las proporciones 10-90 y 90-10 aún con un tamaño de muestra n=500 indica que se rechaza la hipótesis nula, por tanto existe poca probabilidad de que sigan una distribución normal. 

#-------------------------------------------------------------------------------------------
  1. 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.
  1. 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).
Lote1=c(rep("enfermas", 100),rep("sanas", 900))
Lote2=c(rep("enfermas", 150),rep("sanas", 1350))
  1. 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.
diferencia=function(n1){
n2=n1
muestra1= sample(Lote1, size = n1)
P1=sum(muestra1=="enfermas")/n1
muestra2= sample(Lote2, size = n2)
P2=sum(muestra2=="enfermas")/n2
dif=P1-P2
return(dif)
}
  1. 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?
dif_p=sapply(rep(200,500),diferencia)

resumen = summary(dif_p)
Coeficiente_asimetria_Pearson_CAP = 3*(mean(dif_p)-median(dif_p))/sd(dif_p)
Varianza= var(dif_p)
Desviacion = sd(dif_p)
resultados = data.frame(Coeficiente_asimetria_Pearson_CAP, Varianza, Desviacion)
resumen
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.08500 -0.02000  0.00000  0.00006  0.01625  0.08500
resultados
Coeficiente_asimetria_Pearson_CAP Varianza Desviacion
0.0065461 0.0007561 0.0274974 
#--------------------------------------ANÁLISIS---------------------------------------------

# Para una muestra de n=200 y 500 repeticiones, se observa que la media y la mediana son muy similares, esto es un síntoma de simetría que se puede observar en el histograma, en donde la media se representa por la linea de color rojo y la mediana con color azul; el coeficiente de asimetría de Pearson CAP indica que la distribución es simétrica ya que CAP ≈ 0. En el histograma se observa que no existe sesgo (o es muy leve dado que el CAP es muy cercano a cero), la media de la muestra es muy similar a la media real , los datos se distribuyen con una alta simetría. 

#Las diferencias no son siempre cero, de echo existe una desviación (de 0.03 ) aproximadamente y el rango entre el primer y tercer cuartil varia entre -0.07 y -0.04 aproximadamente.

#-------------------------------------------------------------------------------------------

hist(dif_p)
abline(v=median(dif_p), col="blue", lwd=2)
abline(v=mean(dif_p), col="red", lwd=2)
legend(x ="topright",legend = c("media", "mediana"), col = c("red", "blue"),lwd = 4 )

  1. 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?
require(ggplot2)

dif_p_n5 =sapply(rep(5,500),diferencia)
dif_p_n10 =sapply(rep(10,500),diferencia)
dif_p_n15 =sapply(rep(15,500),diferencia)
dif_p_n20 =sapply(rep(20,500),diferencia)
dif_p_n30 =sapply(rep(30,500),diferencia)
dif_p_n50 =sapply(rep(50,500),diferencia)
dif_p_n60 =sapply(rep(60,500),diferencia)
dif_p_n100 =sapply(rep(100,500),diferencia)
dif_p_n200 =sapply(rep(200,500),diferencia)
dif_p_n500 =sapply(rep(500,500),diferencia)

media_dif_n5 = mean(dif_p_n5)
media_dif_n10 = mean(dif_p_n10)
media_dif_n15 = mean(dif_p_n15)
media_dif_n20 = mean(dif_p_n20)
media_dif_n30 = mean(dif_p_n30)
media_dif_n50 = mean(dif_p_n50)
media_dif_n60 = mean(dif_p_n60)
media_dif_n100 = mean(dif_p_n100)
media_dif_n200 = mean(dif_p_n200)
media_dif_n500 = mean(dif_p_n500)

vector_media_dif = c(media_dif_n5, media_dif_n10, media_dif_n15, media_dif_n20, media_dif_n30, media_dif_n50, media_dif_n60,media_dif_n100,media_dif_n200,media_dif_n500 )
tamaño_muestra= c(5,10,15,20,30,50,60,100,200,500)
datos_p_value=data.frame(vector_media_dif,tamaño_muestra)
g10 = ggplot(datos_p_value, aes(x=tamaño_muestra, y=vector_media_dif)) +  geom_line()+ geom_point()

p_value_n5 =shapiro.test(dif_p_n5)$p.value
p_value_n10 =shapiro.test(dif_p_n10)$p.value
p_value_n15 =shapiro.test(dif_p_n15)$p.value
p_value_n20 =shapiro.test(dif_p_n20)$p.value
p_value_n30 =shapiro.test(dif_p_n30)$p.value
p_value_n50 =shapiro.test(dif_p_n50)$p.value
p_value_n60 =shapiro.test(dif_p_n60)$p.value
p_value_n100 =shapiro.test(dif_p_n100)$p.value
p_value_n200 =shapiro.test(dif_p_n200)$p.value
p_value_n500 =shapiro.test(dif_p_n500)$p.value

vector_p_value = c(p_value_n5, p_value_n10, p_value_n15, p_value_n20, p_value_n30, p_value_n50, p_value_n60,p_value_n100,p_value_n200,p_value_n500 )
tamaño_muestra= c(5,10,15,20,30,50,60,100,200,500)
datos_p_value=data.frame(vector_p_value,tamaño_muestra)
g11 = ggplot(datos_p_value, aes(y=vector_p_value, x=tamaño_muestra)) +  geom_line()+ geom_point()
require(ggpubr)
ggarrange(g10, g11,labels = c("media según N", "p-value según N"),ncol = 2, nrow = 1, widths = c(1.5, 1.5))

#--------------------------------------ANÁLISIS---------------------------------------------
#A medida que incrementa el tamaño de muestra N, la media de las diferencias tiende a cero como se observa en el gráfico "media según N", esto es evidencia de que el tamaño de muestra afecta el resultado del analisis y que entre mas grande sea N, existe mas confiabilidad en los resultados. 
#En cuanto a normalidad se puede observar que el p_value incrementa a medida que incrementa el tamaño de muestra como se observa en el gráfico "p-value según N", mientras que entre mas pequeño es el p-value, más se aleja de la hipótesis nula (P=0). Ademas en el gráfico de densidad se observa cómo a medida que aumenta N, se observa forma de campana de Gauss.

#En general se puede concluir con una muestra de N=500 que el tratamiento realizado en el Lote1 presentó un desempeño similar con  respecto al Lote2, ya que al aumentar el tamaño de la muestra la diferencia tiende a cero.
#-------------------------------------------------------------------------------------------
plot(density(dif_p_n5),xlim=c(-0.5, 0.5), ylim=c(0, 30))
lines(density(dif_p_n10), lwd=3, col='blue')
lines(density(dif_p_n15), lwd=3, col='red')
lines(density(dif_p_n20), lwd=3, col='green')
lines(density(dif_p_n30), lwd=3, col='purple')
lines(density(dif_p_n50), lwd=3, col='orange')
lines(density(dif_p_n60), lwd=3, col='yellow')
lines(density(dif_p_n100), lwd=3, col='gray')
lines(density(dif_p_n200), lwd=3, col='cyan')
lines(density(dif_p_n500), lwd=3, col='magenta')
legend(x = "topright",legend = c("N=5","N=10","N=15","N=20","N=30","N=50","N=60","N=100","N=200","N=500"), col = c("black","blue", "red","green", "purple", "orange", "yellow","gray","cyan","magenta"),lwd = 4 )

  1. 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: diferencia de 5%)?
Lote1=c(rep("enfermas", 100),rep("sanas", 900))
Lote2=c(rep("enfermas", 225),rep("sanas", 1275))

require(ggplot2)

dif_p_n5 =sapply(rep(5,500),diferencia)
dif_p_n10 =sapply(rep(10,500),diferencia)
dif_p_n15 =sapply(rep(15,500),diferencia)
dif_p_n20 =sapply(rep(20,500),diferencia)
dif_p_n30 =sapply(rep(30,500),diferencia)
dif_p_n50 =sapply(rep(50,500),diferencia)
dif_p_n60 =sapply(rep(60,500),diferencia)
dif_p_n100 =sapply(rep(100,500),diferencia)
dif_p_n200 =sapply(rep(200,500),diferencia)
dif_p_n500 =sapply(rep(500,500),diferencia)

media_dif_n5 = mean(dif_p_n5)
media_dif_n10 = mean(dif_p_n10)
media_dif_n15 = mean(dif_p_n15)
media_dif_n20 = mean(dif_p_n20)
media_dif_n30 = mean(dif_p_n30)
media_dif_n50 = mean(dif_p_n50)
media_dif_n60 = mean(dif_p_n60)
media_dif_n100 = mean(dif_p_n100)
media_dif_n200 = mean(dif_p_n200)
media_dif_n500 = mean(dif_p_n500)

vector_media_dif = c(media_dif_n5, media_dif_n10, media_dif_n15, media_dif_n20, media_dif_n30, media_dif_n50, media_dif_n60,media_dif_n100,media_dif_n200,media_dif_n500 )
tamaño_muestra= c(5,10,15,20,30,50,60,100,200,500)
datos_p_value=data.frame(vector_media_dif,tamaño_muestra)
g10 = ggplot(datos_p_value, aes(x=tamaño_muestra, y=vector_media_dif)) +  geom_line()+ geom_point()

p_value_n5 =shapiro.test(dif_p_n5)$p.value
p_value_n10 =shapiro.test(dif_p_n10)$p.value
p_value_n15 =shapiro.test(dif_p_n15)$p.value
p_value_n20 =shapiro.test(dif_p_n20)$p.value
p_value_n30 =shapiro.test(dif_p_n30)$p.value
p_value_n50 =shapiro.test(dif_p_n50)$p.value
p_value_n60 =shapiro.test(dif_p_n60)$p.value
p_value_n100 =shapiro.test(dif_p_n100)$p.value
p_value_n200 =shapiro.test(dif_p_n200)$p.value
p_value_n500 =shapiro.test(dif_p_n500)$p.value

vector_p_value = c(p_value_n5, p_value_n10, p_value_n15, p_value_n20, p_value_n30, p_value_n50, p_value_n60,p_value_n100,p_value_n200,p_value_n500 )
tamaño_muestra= c(5,10,15,20,30,50,60,100,200,500)
datos_p_value=data.frame(vector_p_value,tamaño_muestra)
g11 = ggplot(datos_p_value, aes(y=vector_p_value, x=tamaño_muestra)) +  geom_line()+ geom_point()

plot(density(dif_p_n5),xlim=c(-0.5, 0.5), ylim=c(0, 30))
lines(density(dif_p_n10), lwd=3, col='blue')
lines(density(dif_p_n15), lwd=3, col='red')
lines(density(dif_p_n20), lwd=3, col='green')
lines(density(dif_p_n30), lwd=3, col='purple')
lines(density(dif_p_n50), lwd=3, col='orange')
lines(density(dif_p_n60), lwd=3, col='yellow')
lines(density(dif_p_n100), lwd=3, col='gray')
lines(density(dif_p_n200), lwd=3, col='cyan')
lines(density(dif_p_n500), lwd=3, col='magenta')
legend(x = "topright",legend = c("N=5","N=10","N=15","N=20","N=30","N=50","N=60","N=100","N=200","N=500"), col = c("black","blue", "red","green", "purple", "orange", "yellow","gray","cyan","magenta"),lwd = 4 )

require(ggpubr)
ggarrange(g10, g11,labels = c("media según N", "p-value según N"),ncol = 2, nrow = 1, widths = c(1.5, 1.5))

#--------------------------------------ANÁLISIS---------------------------------------------
#A medida que incrementa el tamaño de muestra N, la media de las diferencias tiende a -0.05 como se observa en el gráfico "media según N", esto es evidencia de que el tamaño de muestra afecta el resultado del análisis y que entre mas grande sea N, existe mas confiabilidad en los resultados. 
#En cuanto a normalidad se puede observar que el p_value incrementa a medida que incrementa el tamaño de muestra como se observa en el gráfico "p-value según N", mientras que entre mas pequeño es el p-value, más se aleja de la hipótesis nula (P=-0.05). Ademas en el gráfico de densidad se observa cómo a medida que aumenta N, se observa forma de campana de Gauss.

#En general se puede concluir con una muestra de N=500 que el tratamiento realizado en el Lote1 presentó un mejor desempeño reduciendo alrededor del 5% el porcentaje de enfermos con respecto al Lote2 ya que al aumentar el tamaño de la muestra el valor de la diferencia tiene a -0.05.
#-------------------------------------------------------------------------------------------
summary(dif_p_n500)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.09800 -0.06200 -0.05000 -0.05046 -0.04000  0.00000
  1. 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.

En sus casi nueve décadas de su existencia, los valores p siempre han sufrido criticas relacionadas con lo errores que pueden surgir al utilizar este conocimiento como si fuese una ley. Stephen Ziliak, economista de la Universidad Roosevelt en Chicago, Illinois, afirma que: “Los valores P no están haciendo su trabajo, porque no pueden”. Cuando Ronald Fisher dio a conocer el valor p en la década de 1920 no pretendía que este sea una prueba definitiva, simplemente quería brindar una forma informal que permita saber si la evidencia era significativa para tomar decisiones basadas en datos. Fisher sugirió que cuanto más pequeño era el valor p, mayor era la probabilidad de que la hipótesis nula fuera falsa. Posteriormente, el valor p de 0.05 se consagró como “estadísticamente significativo”.

El termino P-hacking hace referencia a probar cosas hasta obtener el resultado deseado aunque este sea inconsistente. En un análisis se encontró que muchos artículos de psicología publicados reportan valores P que se agrupan de una forma sospechosa alrededor de 0.05, en otras palabras, los investigadores usan p-hacking con el fin de obtener valores p significativos a como de lugar.

Dada la popularidad de usar valor p en investigaciones se ha encontrado que muchos de los hallazgos publicados no son ciertos. Los estadísticos señalan una serie de medidas que pueden ayudar, como conocer el tamaño del efecto y sus intervalos de confianza, ya que esta información transmite lo que no hace el valor p. 

Goodman afirma que los investigadores deben darse cuenta de los límites de las estadísticas convencionales, un método no va a resolver todas las incógnitas que se puedan presentar en un análisis de datos, “Los números son donde la discusión científica debe comenzar, no terminar”.

Como reflexión personal, considero que un análisis estadístico confiable debe hacer uso de todas las herramientas que brinda la estadística, no se puede solo realizar pruebas de hipótesis y concluir o tomar decisiones únicamente con el valor p; es necesario contemplar intervalos de confianza, considerar métodos como la regla de bayes, entre otros, y sobre todo nunca incurrir en p-hacking para arreglar los resultados.