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%.

create_pob=function(N,p){
  
  pob=c(rep("enfermo",N*p),rep("sano",N*(1-p)))
  return(pob)
}

poblacion=create_pob(1000,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.

calc_p_gorro=function(n){
  
muestra=sample(poblacion,n)
p_gorro=sum(muestra=='enfermo')/n
return(p_gorro)

}

calc_p_gorro(60)
## [1] 0.3833333

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?

posibles_p_gorro=sapply(rep(60,500),calc_p_gorro)

descrip_p_0.5 = descriptivas(posibles_p_gorro)

kbl(descrip_p_0.5,caption = "<center><strong>Descriptivas 500 Estimadores P=0.5 t=60
  </strong></center>")%>%   kable_paper(bootstrap_options = "striped")
Descriptivas 500 Estimadores P=0.5 t=60
MEDIDA VALOR
Observaciones 500.00000000
Mínimo 0.25000000
1er Q 0.45000000
Media 0.50030000
Mediana 0.50000000
Desv Est 0.06511572
3er Q 0.55000000
Máximo 0.68333333 
n_muestra=60 

hist(posibles_p_gorro,col="gray",xlab= "proporcion",
     main=paste(c("Histograma Estimación Proporción (P=0.5)","t=",n_muestra,"Rep=",500),collapse = " "), freq = F)

abline(v=0.5,col="red",lwd=4)
abline(v=mean(posibles_p_gorro),col="blue",lwd=4)

Se puede observar gracias a la anterior tabla con información descriptiva y al histograma generado que si usamos un tamaño de muestra equivalente a 60 por ejemplo, encontraremos nuestro \(P\) teórico cercano al \(\hat p\) estimado. Debo aclarar que lo anteriormente encontrado puede verse tambien asociado a coincidencia, pues se esta trabajando un n>30, pero todavia no lo suficientemente grande para afirmar que por teorema de límite central, el \(\hat p\) estimado converge al \(P\) real.

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).

N=500
p_general<-function(t){
  
  p_tabla_gorro=matrix(NA,N)
  for (i in 1:N){
    p_tabla_gorro[i]=sapply(rep(t,N),calc_p_gorro)
    
  }
  return(p_tabla_gorro)  
  
}

P_5=p_general(5)
P_10=p_general(10)
P_15=p_general(15)
P_20=p_general(20)
P_30=p_general(30)
P_50=p_general(50)
P_60=p_general(60)
P_100=p_general(100)
P_200=p_general(200)
P_500=p_general(500)

matrices<-list(P_5,P_10,P_15,P_20,P_30,P_50,P_60,P_100,P_200,P_500)

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

Proporcion (P teorico (rojo) vs P empírico (azul) P=0.5 

par(mfrow=c(2,5))
  for(i in 1:length(t)){
    hist(matrices[[i]][,1], col="gray",xlab= "Proporcion",
         main=paste(c("Hist","t=",t[i],"Rep=",500,"P=",0.5),collapse = " "), freq = F)
    abline(v=0.5,col="red",lwd=4)
    abline(v=mean(matrices[[i]][,1]),col="blue",lwd=4)
  }

Se puede observar en el anterior gráfico que realmente es bastante influyente el ejercicio de aumentar el tamaño de muestra, y mas cuando se desea comparar entre el \(P\) teorico en este caso 0.5, versus el estimado.

Esta diferencia se enmarca mucho mas cuando tenemos por ejemplo tamños de n=5 y n=10. Además que podemos observar que cuando n aumenta, se observa un comportamiento de campana, lo que podría sugerirnos que si se podría comportar de tamaño normal para los gráficos superiores, pero esto evidentemente debe ser sustentado por mas pruebas por lo que acontinuación se muestran un par de pruebas, que nos permitiran concluir lo observado.

Prueba de Normalidad Shapiro-Wilk y Gráficos QQ-Plot

pvalor<-matrix(0,length(t))
  for(i in 1:length(t)){
      pvalor[i]=round(shapiro.test(matrices[[i]][,1])$p.value,3)
  }


prueba=data.frame(pvalor)
prueba$tamaño=c(5,10,15,20,30,50,60,100,200,500)
prueba$proporcion_téorica=c(rep(0.5,10))

colnames(prueba) <- c('p_val_Shapiro','Tamaño n','Proporción téorica')

prueba <- prueba[, c('Proporción téorica', 'Tamaño n', "p_val_Shapiro")]


kbl(prueba,    caption = "<center><strong>Prueba Shapiro </strong></center>") %>%
  kable_paper(bootstrap_options = "striped", full_width = FALSE,position = "left")
Prueba Shapiro
Proporción téorica Tamaño n p_val_Shapiro
0.5 5 0.000
0.5 10 0.000
0.5 15 0.000
0.5 20 0.000
0.5 30 0.000
0.5 50 0.008
0.5 60 0.012
0.5 100 0.033
0.5 200 0.085
0.5 500 0.145 
p_teorico=0.5

par(mfrow=c(2,5))
for(j in 1:length(t)){
  qqnorm(matrices[[j]][,1], col ='blue', main=paste('QQ P=',p_teorico ,'n=', t[j] ), ylab = 'Cuantiles Muestrales',
         xlab = 'Cuantiles Teoricos')
  qqline(matrices[[j]][,1])

}

Como se planteó al inicio de este ejercicio práctico, el propósito es poder demostrar el teorema del límite central, en donde nos sugiere que a medida que nuestro tamaño de muestra aumenta los datos pueden presentar una distribución normal. Se puede observar en los gráficos generados por qqplot,(gráfico en el cual se evalua los valores de la proporción encontrados por simulación vs los valores estandarizados téoricos) que a medida que el tamaño de muestra va incrementando se va observando un mejor comportamiento rectilineo de los puntos, lo que en ese caso para lo observado por estos gráficos, si influye o cobra relvancia trabajar con un tamaño de muestra mas grande y que además se podría sugerir que si presenta una distibución apróximadamente normal.

Para poder acompañar lo que visualmente se puede observar a través del qqplot es importante realizar la prueba de shapiro wilk, que plantea

\(H_0: Sigue \, una\, distribución \,normal\)

\(H_a: No \,sigue \, una\, distribución \,normal\)

Se puede observar que cuando \(n=200\) y \(n=500\) los valores p encontrados estan por encima de un nivel de significancia \(\alpha=0.05\), lo que nos da evidencia estadística para no rechazar \(H_0\), por ende podría pensarse que si presentan una distribución normal.

Aclaración: En lo anteriormente comentado supongo un \(\alpha=0.05\), pero si evaluaramos con un nivel de confianza del 90%, es decir \(\alpha=0.1\) la interpretación podría variar.

Intervalos de confianza para Proporción 0.5 

medias_proporcion<-matrix(NA,10)
for (i in 1:length(t)){
  medias_proporcion[i]=mean(matrices[[i]][,1])
  
}


Ic_confianza<-matrix(NA,10,2)
for (i in 1:length(t)){
  if (t[i]<=30){
  Ic_confianza[i,1]=round(medias_proporcion[i]-round(qt(0.025,t[i]-1 , lower.tail = F),2)*(sqrt((medias_proporcion[i]*(1-medias_proporcion[i])/t[i]))),2)
  Ic_confianza[i,2]=round(medias_proporcion[i]+round(qt(0.025,t[i]-1 , lower.tail = F),2)*(sqrt(((medias_proporcion[i]*(1-medias_proporcion[i]))/t[i]))),2)}
  else {
    Ic_confianza[i,1]=round(medias_proporcion[i]-round(qnorm(0.975,mean=0,sd=1),2)*(sqrt(((medias_proporcion[i]*(1-medias_proporcion[i]))/t[i]))),2)
    Ic_confianza[i,2]=round(medias_proporcion[i]+round(qnorm(0.975,mean=0,sd=1),2)*(sqrt(((medias_proporcion[i]*(1-medias_proporcion[i]))/t[i]))),2)
    
  }
  
}



Ic_confianza=data.frame(Ic_confianza)
Ic_confianza$Tamaño=t
Ic_confianza$P=rep(0.5,10)
colnames(Ic_confianza)<-c("IC_inf",'IC_sup','Tamaño','P')
Ic_confianza$Amplitud_IC=Ic_confianza$IC_sup-Ic_confianza$IC_inf

kbl(Ic_confianza,    caption = "<center><strong>Intervalo de Confianza P=0.5 </strong></center>") %>%
  kable_paper(bootstrap_options = "striped", full_width = FALSE,position = "left")
Intervalo de Confianza P=0.5
IC_inf IC_sup Tamaño P Amplitud_IC
-0.12 1.13 5 0.5 1.25
0.14 0.85 10 0.5 0.71
0.22 0.77 15 0.5 0.55
0.27 0.74 20 0.5 0.47
0.32 0.70 30 0.5 0.38
0.36 0.64 50 0.5 0.28
0.37 0.62 60 0.5 0.25
0.40 0.60 100 0.5 0.20
0.43 0.57 200 0.5 0.14
0.46 0.54 500 0.5 0.08 
plot(x<-c(5  ,10 , 15 , 20 , 30 , 50 , 60, 100, 200 ,500),Ic_confianza$Amplitud_IC,
     xlab='Tamaños Muestrales',ylab='Amplitud IC',main="Amplitud Intervalo vs Tamaño de Muestra P=0.5", 
     cex.main=1,type='b' ,col = "blue")

Una buena forma de poder concluir lo encontrado además del uso de pruebas estadísticas es observar por ejemplo la amplitud del intervalo de confianza. Para este caso decidi calcular el promedio de \(\hat p_{1:500}\) para cada uno de los tamaños evaluados. Recordemos un poco como se define el intervalo de confianza para la proporción.

\(IC_{\hat p} = <\hat p \, -+ \, t_{(\frac{\alpha}{2},gl)}*\sqrt{\frac{\hat p(1-\hat p)}{n}}>\)

Donde \(gl\) los grados de libertad se expresan como \(n-1\).

\(IC_{\hat p} = <\hat p \, -+ \, z_{\frac{\alpha}{2}}*\sqrt{\frac{\hat p(1-\hat p)}{n}}>\)

Como se puede observar tanto en la tabla como en el gráfico anteriormente generado, podríamos decir que para todo los p evaluados, si por ejemplo tuvieramos un tamaño de muestra igual 100, que: el 95% de las veces en que se vuelva a repetir este mismo ejercicio, se espera que la proporción real de la población se encuentre entre 0.40 y 0.60. Cada vez que se incrementa el tamaño de muestra este intervalo suele disminuir, lo que de transfondo explica, por que al aumentar el n, tendremos valores de p estimado mas precisos y mas cercano al p real.

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

Proporcion P=0.1 

N_pob=1000
p=0.1


create_pob=function(N_pob,p){
  
  pob=c(rep("enfermo",N_pob*p),rep("sano",N_pob*(1-p)))
  return(pob)
}

poblacion=create_pob(N_pob,p)


calc_p_gorro=function(n){
  
  muestra=sample(poblacion,n)
  p_gorro=sum(muestra=='enfermo')/n
  return(p_gorro)
  
}

calc_p_gorro(60)
## [1] 0.1166667
n_muestra=60 

posibles_p_gorro=sapply(rep(n_muestra,500),calc_p_gorro)

hist(posibles_p_gorro,col="gray",xlab= "proporcion",
     main=paste(c("Histograma Estimación Proporción (P=0.1)","t=",n_muestra,"Rep=",500),collapse = " "), freq = F)

abline(v=0.1,col="red",lwd=4)
abline(v=mean(posibles_p_gorro),col="blue",lwd=4)

N_veces=500
p_general<-function(t){
  
  p_tabla_gorro=matrix(NA,N_veces)
  for (i in 1:N){
    p_tabla_gorro[i]=sapply(rep(t,N_veces),calc_p_gorro)
    
  }
  return(p_tabla_gorro)  
  
}

P_5=p_general(5)
P_10=p_general(10)
P_15=p_general(15)
P_20=p_general(20)
P_30=p_general(30)
P_50=p_general(50)
P_60=p_general(60)
P_100=p_general(100)
P_200=p_general(200)
P_500=p_general(500)

matrices<-list(P_5,P_10,P_15,P_20,P_30,P_50,P_60,P_100,P_200,P_500)

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



par(mfrow=c(2,5))
#mtext("Proporcion (P teorico (azul) vs P empírico (rojo) P=",p)
for(i in 1:length(t)){
  hist(matrices[[i]][,1], col="gray",xlab= "Proporcion",
       main=paste(c("Hist","t=",t[i],"Rep=",500,"P=",p),collapse = " "), freq = F)
  abline(v=p,col="red",lwd=4)
  abline(v=mean(matrices[[i]][,1]),col="blue",lwd=4)
}

pvalor<-matrix(0,length(t))
for(i in 1:length(t)){
  pvalor[i]=round(shapiro.test(matrices[[i]][,1])$p.value,3)
}

prueba=data.frame(pvalor)
prueba$tamaño=c(5,10,15,20,30,50,60,100,200,500)
prueba$proporcion_téorica=c(rep(p,10))

colnames(prueba) <- c('p_val_Shapiro','Tamaño n','Proporción téorica')

prueba <- prueba[, c('Proporción téorica', 'Tamaño n', "p_val_Shapiro")]

kbl(prueba,    caption = "<center><strong>Prueba Shapiro </strong></center>") %>%
  kable_paper(bootstrap_options = "striped", full_width = FALSE,position = "left")
Prueba Shapiro
Proporción téorica Tamaño n p_val_Shapiro
0.1 5 0.000
0.1 10 0.000
0.1 15 0.000
0.1 20 0.000
0.1 30 0.000
0.1 50 0.000
0.1 60 0.000
0.1 100 0.000
0.1 200 0.000
0.1 500 0.001 
par(mfrow=c(2,5))
for(j in 1:length(t)){
  qqnorm(matrices[[j]][,1], col ='blue', main=paste('QQP=',p ,' n=', t[j], sep = ), ylab = 'Cuantiles Muestrales',
         xlab = 'Cuantiles Teoricos')
  qqline(matrices[[j]][,1])}

Intervalos de confianza para Proporción 0.1 

medias_proporcion<-matrix(NA,10)
for (i in 1:length(t)){
  medias_proporcion[i]=mean(matrices[[i]][,1])
  
}

P=0.1
Q=0.9

Ic_confianza<-matrix(NA,10,2)
for (i in 1:length(t)){
  if (t[i]<=30){
  Ic_confianza[i,1]=round(medias_proporcion[i]-round(qt(0.025,t[i]-1 , lower.tail = F),2)*(sqrt(((medias_proporcion[i]*(1-medias_proporcion[i]))/t[i]))),2)
  Ic_confianza[i,2]=round(medias_proporcion[i]+round(qt(0.025,t[i]-1 , lower.tail = F),2)*(sqrt(((medias_proporcion[i]*(1-medias_proporcion[i]))/t[i]))),2)}
  else {
    Ic_confianza[i,1]=round(medias_proporcion[i]-round(qnorm(0.975,mean=0,sd=1),2)*(sqrt(((medias_proporcion[i]*(1-medias_proporcion[i]))/t[i]))),2)
    Ic_confianza[i,2]=round(medias_proporcion[i]+round(qnorm(0.975,mean=0,sd=1),2)*(sqrt(((medias_proporcion[i]*(1-medias_proporcion[i]))/t[i]))),2)
    
  }
  
}



Ic_confianza=data.frame(Ic_confianza)
Ic_confianza$Tamaño=t
Ic_confianza$P=rep(0.1,10)
colnames(Ic_confianza)<-c("IC_inf",'IC_sup','Tamaño','P')
Ic_confianza$Amplitud_IC=Ic_confianza$IC_sup-Ic_confianza$IC_inf

kbl(Ic_confianza,    caption = "<center><strong>Intervalo de Confianza P=0.1 </strong></center>") %>%
  kable_paper(bootstrap_options = "striped", full_width = FALSE,position = "left")
Intervalo de Confianza P=0.1
IC_inf IC_sup Tamaño P Amplitud_IC
-0.27 0.47 5 0.1 0.74
-0.11 0.31 10 0.1 0.42
-0.07 0.26 15 0.1 0.33
-0.04 0.24 20 0.1 0.28
-0.01 0.21 30 0.1 0.22
0.02 0.19 50 0.1 0.17
0.02 0.17 60 0.1 0.15
0.04 0.16 100 0.1 0.12
0.06 0.14 200 0.1 0.08
0.07 0.13 500 0.1 0.06 
plot(x<-c(5  ,10 , 15 , 20 , 30 , 50 , 60, 100, 200 ,500),Ic_confianza$Amplitud_IC,
     xlab='Tamaños Muestrales',ylab='Amplitud IC',main="Amplitud Intervalo vs Tamaño de Muestra P=0.1", 
     cex.main=1,type='b' ,col = "blue")

Proporcion P=0.9 

N_pob=1000
p=0.9


create_pob=function(N_pob,p){
  
  pob=c(rep("enfermo",N_pob*p),rep("sano",N_pob*(1-p)))
  return(pob)
}

poblacion=create_pob(N_pob,p)


calc_p_gorro=function(n){
  
  muestra=sample(poblacion,n)
  p_gorro=sum(muestra=='enfermo')/n
  return(p_gorro)
  
}

calc_p_gorro(60)
## [1] 0.8666667
n_muestra=60 

posibles_p_gorro=sapply(rep(n_muestra,500),calc_p_gorro)

hist(posibles_p_gorro,col="gray",xlab= "proporcion",
     main=paste(c("Histograma Estimación Proporción (P=0.9)","t=",n_muestra,"Rep=",500),collapse = " "), freq = F)

abline(v=0.9,col="red",lwd=4)
abline(v=mean(posibles_p_gorro),col="blue",lwd=4)

N_veces=500
p_general<-function(t){
  
  p_tabla_gorro=matrix(NA,N_veces)
  for (i in 1:N){
    p_tabla_gorro[i]=sapply(rep(t,N_veces),calc_p_gorro)
    
  }
  return(p_tabla_gorro)  
  
}

P_5=p_general(5)
P_10=p_general(10)
P_15=p_general(15)
P_20=p_general(20)
P_30=p_general(30)
P_50=p_general(50)
P_60=p_general(60)
P_100=p_general(100)
P_200=p_general(200)
P_500=p_general(500)

matrices<-list(P_5,P_10,P_15,P_20,P_30,P_50,P_60,P_100,P_200,P_500)

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



par(mfrow=c(2,5))
#mtext("Proporcion (P teorico (azul) vs P empírico (rojo) P=",p)
for(i in 1:length(t)){
  hist(matrices[[i]][,1], col="gray",xlab= "Proporcion",
       main=paste(c("Hist","t=",t[i],"Rep=",500,"P=",p),collapse = " "), freq = F)
  abline(v=p,col="red",lwd=4)
  abline(v=mean(matrices[[i]][,1]),col="blue",lwd=4)
}

pvalor<-matrix(0,length(t))
for(i in 1:length(t)){
  pvalor[i]=round(shapiro.test(matrices[[i]][,1])$p.value,3)
}

prueba=data.frame(pvalor)
prueba$tamaño=c(5,10,15,20,30,50,60,100,200,500)
prueba$proporcion_téorica=c(rep(p,10))

colnames(prueba) <- c('p_val_Shapiro','Tamaño n','Proporción téorica')

prueba <- prueba[, c('Proporción téorica', 'Tamaño n', "p_val_Shapiro")]

kbl(prueba,    caption = "<center><strong>Prueba Shapiro </strong></center>") %>%
  kable_paper(bootstrap_options = "striped", full_width = FALSE,position = "left")
Prueba Shapiro
Proporción téorica Tamaño n p_val_Shapiro
0.9 5 0.000
0.9 10 0.000
0.9 15 0.000
0.9 20 0.000
0.9 30 0.000
0.9 50 0.000
0.9 60 0.000
0.9 100 0.000
0.9 200 0.001
0.9 500 0.039 
par(mfrow=c(2,5))
for(j in 1:length(t)){
  qqnorm(matrices[[j]][,1], col ='blue', main=paste('QQ P=',p ,' n=', t[j], sep = ), ylab = 'Cuantiles Muestrales',
         xlab = 'Cuantiles Teoricos')
  qqline(matrices[[j]][,1])}

Intervalos de confianza para Proporción 0.9

Después de observar los resultados en conjunto en el cual se evalua tres escenarios, entre los cuales se varia la proporción 0.10, 0.5 y 0.9 Se puede concluir que utilizar un tamaño de muestra mas grande si respalda la idea de que la distribución de proporción se aproxime a una distribución normal.

Evidentemente para poder evaluar que sea un resultado con evidencia estadística, se llevaron a cabo pruebas estadísticas como shapiro wilk, en el cual, destaco que cuando aumentaba la proporción era más dificil encontrar normalidad. Por ejemplo cuando n=500 y P=0.9, nos sugería un valor P menor a cualquier nivel de significancia tradicionalmente usado (0.01, 0.05 ó 0.1).

Para cuando se observa lo encontrado cuando la proporción es igual a 0.1, cuando n=500, nos sugiere un valor p de aproximadamente 0.2, que con este resultado no rechazariamos H0, por ende pdría distribuirse de forma aproximadamente normal.

Es importante aclarar que aunque la proporción como sabemos se mueve entre 0 y 1, cuando se trabaja con tamaños de muestra muy pequeños, se puede correr el riesgo de que, el limite inferior del intervalo este por debajo de 0, pues lo que esta sucediendo es que la proporción estimada este siendo diferente contra la proporcion teorica, es decir el nivel de confianza que estemos usando impacta en los resultados obtenidos. Esto sucede, por ejemplo para P=0.5 cuando se trabaja con un n=5.

Punto 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.

Escenario 1: P1=P2

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 

p_dif=0.1

create_pob=function(N,p){
  
  pob=c(rep("enfermo",N*p),rep("sano",N*(1-p)))
  return(pob)
}

poblacion_1=create_pob(1000,p_dif)
poblacion_2=create_pob(1500,p_dif)

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 

calc_dif_p_gorro=function(n_muestra){
  
  muestra_1=sample(poblacion_1,n_muestra)
  muestra_2=sample(poblacion_2,n_muestra)
  
  p_gorro_1=sum(muestra_1=='enfermo')/n_muestra
  p_gorro_2=sum(muestra_2=='enfermo')/n_muestra
  
  
  dif_p_gorro=p_gorro_1-p_gorro_2
  return(dif_p_gorro)
  
}

calc_dif_p_gorro(60)
## [1] -0.01666667

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? 

posibles_p_gorro_dif=sapply(rep(60,500),calc_dif_p_gorro)

descrip_p_0.1 = descriptivas(posibles_p_gorro_dif)

kbl(descrip_p_0.1,caption = "<center><strong>Descriptivas 500 Estimadores P1=P2=0.1 t=60
  </strong></center>")%>%   kable_paper(bootstrap_options = "striped")
Descriptivas 500 Estimadores P1=P2=0.1 t=60
MEDIDA VALOR
Observaciones 500.0000000000
Mínimo -0.1500000000
1er Q -0.0333333333
Media 0.0003333333
Mediana 0.0000000000
Desv Est 0.0530929368
3er Q 0.0333333333
Máximo 0.2000000000

Aunque el promedio de las diferencia de las proporciones es practicamente 0, se puede observar en la anterior tabla que puede alcanzar valores como minimo de -0.15 y como máximo valores hasta 0.15.(esto puede variar al no tener una muestra fija).

hist(posibles_p_gorro_dif,col="gray",xlab= "diferencia_proporcion",
     main=paste(c("Histograma Diferencia de Proporciones(P1:0.10, P2:0.10)","n=",n_muestra,"Rep=",500),collapse = " "), freq = F)
abline(v=0,col="red",lwd=4)
abline(v=mean(posibles_p_gorro_dif),col="blue",lwd=4)

Se puede observar gracias al histograma y a las descriptivas calculadas que cuando se evalua la diferencia de proporciones en cuanto a las plantas enfermas encontradas en el lote 1 versus el lote 2, que la diferencia es cercana a 0(valor téorico) resultante de \(P1-P2=0\) pero aun se observa un diferencia con respecto a la diferencia de proporción media encontrada con un n=60 y con 500 repeticiones.

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?

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

p_general<-function(t){
  
  p_tabla_dif_gorro=matrix(NA,N_veces)
  for (i in 1:N){
    p_tabla_dif_gorro[i]=sapply(rep(t,N_veces),calc_dif_p_gorro)
    
  }
  return(p_tabla_dif_gorro)  
  
}

P_5=p_general(5)
P_10=p_general(10)
P_15=p_general(15)
P_20=p_general(20)
P_30=p_general(30)
P_50=p_general(50)
P_60=p_general(60)
P_100=p_general(100)
P_200=p_general(200)
P_500=p_general(500)

matrices<-list(P_5,P_10,P_15,P_20,P_30,P_50,P_60,P_100,P_200,P_500)


par(mfrow=c(2,5))
for(i in 1:length(t)){
  hist(matrices[[i]][,1], col="gray",xlab= "diferencia_proporcion",
       main=paste(c("Hist Dif_prop Esc1","n=",t[i],"Rep=",500,"p=",p_dif),collapse = " "), freq = F)
  abline(v=0,col="red",lwd=4)
  abline(v=mean(matrices[[i]][,1]),col="blue",lwd=4)
}

pvalor<-matrix(0,length(t))
for(i in 1:length(t)){
  pvalor[i]=round(shapiro.test(matrices[[i]][,1])$p.value,3)
}

prueba=data.frame(pvalor)
prueba$tamaño=c(5,10,15,20,30,50,60,100,200,500)
prueba$proporcion_téorica=c(rep(p_dif,10))

colnames(prueba) <- c('p_val_Shapiro','Tamaño n','Proporción Enfermos P1=P2')

prueba <- prueba[, c("Proporción Enfermos P1=P2", 'Tamaño n', "p_val_Shapiro")]

kbl(prueba,    caption = "<center><strong>Prueba Shapiro Escenario 1 P1=P2)</strong></center>") %>%
  kable_paper(bootstrap_options = "striped", full_width = FALSE,position = "left")
Prueba Shapiro Escenario 1 P1=P2)
Proporción Enfermos P1=P2 Tamaño n p_val_Shapiro
0.1 5 0.000
0.1 10 0.000
0.1 15 0.000
0.1 20 0.000
0.1 30 0.000
0.1 50 0.000
0.1 60 0.000
0.1 100 0.002
0.1 200 0.004
0.1 500 0.489 
par(mfrow=c(2,5))
for(j in 1:length(t)){
  qqnorm(matrices[[j]][,1], col ='blue', main=paste('QQ P1=P2',p_dif ,' n=', t[j]), ylab = 'Cuantiles Muestrales',
         xlab = 'Cuantiles Teoricos')
  qqline(matrices[[j]][,1])
  
}

Cuando se hace el ejercicio de poder comparar para diferentes tamaños de muestra que van desde 5 hasta 500, que a medida que se va incrementando los tamaños muestrales, la diferencia de la proporción estimada, se va a acercando a la diferencia de proporcion real igual a 0. Algo que podría corroborar lo anterior es si vemos los histogramas, cada vez el rango de valores que se toman en el eje x se va haciendo mas reducido, es decir cercano a 0, lo que estaría disminuyendo la posible variabilidad que podemos encontrar cuando se trabajan con tamaños de muestra pequeños, como 5, 10 o 15 por ejemplo. La distribución muestral se va haciendo mucho más simétrica y con forma mas acampanada.

En segunda instancia, al aplicar la prueba de shapiro wilk, se puede observar, que cuando \(n=500\), nos sugiere un \(valor \;p=0.09\) lo que podríamos sugerir que solo para este tamaño, si se compara con un \(\alpha=0.05\) no se rechazaria la hipotésis nula que sugiere normalidad.

En cuanto a los QQ-Plot se puede observar, que de igual manera cuando se evaluó la proporción, en este caso, en diferencia de proporciones, a medida que aumenta el tamaño de muestra, se evidencia un comportamiento rectilineo, lo que nos podría sugerir un comportamiento normal pero solo cuando n es lo suficientemente grande.

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: diferencia de 5%)? 

#a.

p1_dif=0.1
p2_dif=0.15

create_pob=function(N,p){
  
  pob=c(rep("enfermo",N*p),rep("sano",N*(1-p)))
  return(pob)
}

poblacion_1=create_pob(1000,p1_dif)
poblacion_2=create_pob(1500,p2_dif)



# b.

calc_dif_p_gorro=function(n_muestra){
  
  muestra_1=sample(poblacion_1,n_muestra)
  muestra_2=sample(poblacion_2,n_muestra)
  
  p_gorro_1=sum(muestra_1=='enfermo')/n_muestra
  p_gorro_2=sum(muestra_2=='enfermo')/n_muestra
  
  
  dif_p_gorro=p_gorro_1-p_gorro_2
  return(dif_p_gorro)
  
}

calc_dif_p_gorro(60)
## [1] -0.08333333
# c.

n_muestra=60
posibles_p_gorro_dif=sapply(rep(n_muestra,500),calc_dif_p_gorro)


hist(posibles_p_gorro_dif,col="gray",xlab= "diferencia_proporcion",
     main=paste(c("Histograma Diferencia de Proporciones(P1:0.10, P2:0.15)","n=",n_muestra,"Rep=",500),collapse = " "), freq = F)
abline(v=-0.05,col="red",lwd=4)
abline(v=mean(posibles_p_gorro_dif),col="blue",lwd=4)

#d
t=c(5,10,15,20,30,50,60,100,200,500)
N_veces=500

p_general<-function(t){
  
  p_tabla_dif_gorro=matrix(NA,N_veces)
  for (i in 1:N){
    p_tabla_dif_gorro[i]=sapply(rep(t,N_veces),calc_dif_p_gorro)
    
  }
  return(p_tabla_dif_gorro)  
  
}

P_5=p_general(5)
P_10=p_general(10)
P_15=p_general(15)
P_20=p_general(20)
P_30=p_general(30)
P_50=p_general(50)
P_60=p_general(60)
P_100=p_general(100)
P_200=p_general(200)
P_500=p_general(500)

matrices<-list(P_5,P_10,P_15,P_20,P_30,P_50,P_60,P_100,P_200,P_500)


par(mfrow=c(2,5))
for(i in 1:length(t)){
  hist(matrices[[i]][,1], col="gray",xlab= "diferencia_proporcion",
       main=paste(c("Hist Dif_prop Esc2","n=",t[i]),collapse = " "), freq = F)
  abline(v=-0.05,col="red",lwd=4)
  abline(v=mean(matrices[[i]][,1]),col="blue",lwd=4)
}

pvalor<-matrix(0,length(t))
for(i in 1:length(t)){
  pvalor[i]=round(shapiro.test(matrices[[i]][,1])$p.value,3)
}

prueba=data.frame(pvalor)
prueba$tamaño=c(5,10,15,20,30,50,60,100,200,500)
prueba$proporcion_grupo_1=c(rep(0.1,10))
prueba$proporcion_grupo_2=c(rep(0.15,10))

colnames(prueba) <- c('p_val_Shapiro','Tamaño n','proporcion_grupo_1','proporcion_grupo_2')

prueba <- prueba[, c("proporcion_grupo_1","proporcion_grupo_2", 'Tamaño n', "p_val_Shapiro")]

kbl(prueba,    caption = "<center><strong>Prueba Shapiro Escenario 2 P1<P2)</strong></center>") %>%
  kable_paper(bootstrap_options = "striped", full_width = FALSE,position = "left")
Prueba Shapiro Escenario 2 P1
proporcion_grupo_1 proporcion_grupo_2 Tamaño n p_val_Shapiro
0.1 0.15 5 0.000
0.1 0.15 10 0.000
0.1 0.15 15 0.000
0.1 0.15 20 0.000
0.1 0.15 30 0.000
0.1 0.15 50 0.001
0.1 0.15 60 0.001
0.1 0.15 100 0.076
0.1 0.15 200 0.047
0.1 0.15 500 0.145 
par(mfrow=c(2,5))
for(j in 1:length(t)){
  qqnorm(matrices[[j]][,1], col ='blue', main=paste('QQ P1<P2' ,'n=', t[j], sep = ), ylab = 'Cuantiles Muestrales',
         xlab = 'Cuantiles Teoricos')
  qqline(matrices[[j]][,1])
  
}

En el escenario 2, que nos sugiere que las proporciones teoricas del lote 1 y el lote 2, son diferentes, ($P_1=0.1 , , P_2=0.15 $), se puede evidenciar que en el gráfico del histograma para la diferencia de proporciones con n=60, que la diferencia de proporciones teorica dista un poco de la diferencia de proporciones estimada (esto nos lo permite saber la distancia que hay entre la linea roja y la linea azul). Tengamos en cuenta, que en este caso la diferencia de proporciones real es del 5%.

De igual manera cuando se hace el ejercicio de simulación bajo los diferentes escenarios que nacen de variar el tamaño de muestra, cada vez que se incrementa el tamaño de muestra puede observarse, que efectivamente, parece que la distribucion muestral de la diferencia de proporciones con 500 repeticiones, puede aproximarse a una distribución normal.

Para corroborar lo anteriormente observado, se hace uso de pruebas como shapiro wilk, la cual nos permite evaluar normalidad. Bajo los distintos tamaños de muestra encontrados, se observa que a medida que se incrementa el tamaño de muestra, parece que hubiese evidencia para no rechazar la hipotesis nula, y pensar que hay evidencia sifnificativamente estadística para pensar en normalidad. Pero es curioso que aunque con \(n=200\) sugiere que se rechace \(H_0\) a un nivel de significancia de 0.05, esto no sea lo mismo que se sugiera cuando \(n=500\).

Para los gráficos de QQ-plot se puede observar que a medida que hay un incremento de n, existe un comportamiento rectilineo , entre los cuantiles muestrales y teoricos, lo que podria sugerir, que solo cuando n tiende a ser mas grande (\(n \to \infty\)), podríamos hablar con mas certeza de normalidad, además que las pruebas lo respaldarian.

Para poder decir si hubo una mejora o no en el tratamiento aplicado en el lote 1, versus el lote 2, respecto a la cantidad de plantas enfermas, es necesario hacer medidas como el intervalo de confianza para la diferencia de proporciones que resultados, en el que incluyera el 0 o el no, podría demostrar si hubo diferencia o no en ambas proporciones mas alla del diferencia de la proporcion real, que ya existe.

Punto 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.

El papel del Valor P en la toma de decisiones

A partir de las dos lecturas relacionadas con el valor P, tituladas:

-Statisticians issue warning on P values.

-Statitical Errors: P values, the ‘gold standard’ of statistical validity, are not as reliable as many scientists assume.

Se puede observar que ambos papers buscan hacer una crítica frente al uso del valor P. En el primero, se quiere evidenciar que el uso que se hace del valor P, no debe pensarse como una simple ejecución para poder realizar un contraste de hipotésis o evidenciar si existe una correlación entre un par de variables. Aunque hace muchos años se ha identificado lo peligroso que puede ser trabajar el valor p, como un simple método aislado, no se ha hecho suficiente hincapié en la importancia de entender que hay detrás del número que nos arroja. Es por esto que la Asociación Americana de Estadística, aproximadamente en el año 2016, levanta la mano y decide hacer una aclaración sobre como realmente debe ser su uso. Para ellos, no es suficiente con tratar de demostrar una robustez y reproducibilidad de resultados, con simplemente el uso del valor P. Al contrario, sugieren que se debe ser completamente honestos en lo que se reporta en los diferentes proyectos que se use, además de utilizar el valor P, acompañado de otros análisis descriptivos y cuantitativos que respalden lo encontrado.

Hace énfasis en que tanto puede cambiar nuestro afán por usar y mostrar en diferentes publicaciones el valor p, si de verdad entendiéramos que hay detrás de el. Sugiere que un valor P cuando lo tenemos en niveles llamemos normalmente trabajados (0.01, 0.05 ó 0.1) no significa que la hipótesis sea cierta, por el contrario, nos permite no rechazar la hipótesis nula, y las suposiciones que la acompañan son válidas a su vez, existe una probabilidad del (1% ,5%, 10%)de encontrar un resultado casi imposible como el que se encontró.

Si no hubiese una creencia rotunda sobre lo que nos reporta este valor, es decir hubiesen dudas sobre lo encontrado, permitiría entender la necesidad que hay detrás de hacer más pruebas y de no asumir los primeros hallazgos como determinantes.

En la misma vía el segundo artículo, busca demostrar a partir de diferentes ejemplos y opiniones la relevancia que toma el valor P, cuando se desea hablar de reproducibilidad. Es así como el ejemplo que más me impacto, fue el proyecto llevado a cabo por Matt Motyl, estudiante de Doctorado en Psicología, de una prestigiosa universidad de Estados Unidos, en donde había encontrado con lo que allí puede pensarse como un contundente Valor P (0.01), que había evidencia para pensar que los extremistas ven resultados a blanco o negro, pero que no existen puntos medios. Muy prontamente, se iba a llevar una gran decepción cuando al intentar replicar este estudio, se topó con un valor P de 0.54, lo que prontamente iba a desvirtuar lo primeramente encontrado. Lo que él no se imaginaba, era que basar la reproducibilidad del primer resultado encontrado sobre el Valor P, iba a tener un sustento poco sólido y estable.

¿Cómo hacer para controlar este tipo de situaciones, en las cuales los científicos, se ven persuadidos por utilizar el valor P, de forma poco consciente, y más bien apresurada? Sugieren que se muestren resultados, inclusive como tamaño de muestra a usar, descriptivas sobre los datos que se van a trabajar.

Bien se especifica que la idea de Ronald Fisher, no fue que el valor P, fuera la “última palabra”, sino por el contrario que se considerara como un primer hallazgo, el cual fuera el primer paso, para evaluar de otras formas los resultados. Por esos mismos años Fisher, tenía a su alrededor otros estadísticos que no estaban completamente de acuerdo con su teoría, pero en la cual enfocaron sus esfuerzos y energías en desvirtuar otras ideas, pero no la concebida por el valor P.

Algunos estadísticos como, Cumming cree que los investigadores deben informar el tamaño de los efectos y además el nivel de confianza que hay detrás de los que se está asegurando, pues estas medidas pueden darnos una percepción de la magnitud, y de la importancia relativa de un efecto.

Aunque el valor P se encuentra hoy en día aún muy arraigado, generando en investigadores y científicos una falsa fe en su resultado, la única forma de cambiar este proceder, es definitivamente, generar preguntas asociadas, de cómo se está usando y que se quiere responder cuando hacemos uso de él.