Taller unidad 2

Inferencia estadística y simulación

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

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

lote=c(rep("enfermo",500),rep("sana",500))

Este sera un ejercicio en el cual se supone que se conoce la distribución real de la población para posteriormente mediane simulación ver que tanto nos aproximamos a esta realidad.

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(lote,size = n)

# porcentaje de plantas enfermas
p_gorro=sum(muestra=="enfermo")/n

return(p_gorro)

}

En este caso nuestro estimador de interes será el porcentaje de plantas enfermas dentro de la población, tomaremos diferentes muestras para analizar como es el comportamiento del estimador en relación al parametro de la población, para diferentes tamaños de muestras.

Ejemplo n=5

calc_p_gorro(n = 5)

## [1] 0.2

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?

# muestra de de tamaño 100, repitiendose 500 veces
posibles_p_gorro=sapply(rep(500,500), calc_p_gorro)

hist(posibles_p_gorro)

mean(posibles_p_gorro)

## [1] 0.499688

sd(posibles_p_gorro)

## [1] 0.016522

Conclusión

se evidencia que el estimador presenta un comportamiento cuya función de probabilidad se asemeja a una distribución normal, no presentan sesgos pues estan centrados en el parametro que conocemos que es de 50% (0.5), los posibles valores del estimador presentan desviación baja es decir la variabilidad no es alta y se acercan al valor real, nos daria una buena aproximación.

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

H0: el porcentaje de plantas enfermas sigue una distribución normal.

H1: el porcentaje de plantas enfermas no sigue una distribución normal.

Nivel de significancia: 0.05 (Hipotético).

Gráficos Cuantil-Cuantil (Q-Q plots): Un gráfico Cuantil-Cuantil permite observar cuan cerca está la distribución de un conjunto de datos a alguna distribución ideal

para n=5, proporción población plantas enfermas 50%,50%

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

require(car)

## Loading required package: car

## Loading required package: carData

par(mfrow=c(1,3))
hist(posibles_p_gorro, xlab = "Porcentaje de plantas enfermas", ylab = "Frecuencia", las=1, main = "", col = "gray")
plot(density(posibles_p_gorro), xlab = "Porcentaje de plantas enfermas", ylab = "Densidad", las=1, main = "")
qqPlot(posibles_p_gorro, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales", las=1,main="")

## [1]  15 119

Graficamente se puede ver que se rechaza la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia.

Valor medio y desviación estandar

mean(posibles_p_gorro)

## [1] 0.4872

sd(posibles_p_gorro)

## [1] 0.2238221

Si bien se acerca el valor medio a 0,5 que es el valor real la desviación estandar de los datos es muy alta, lo cual refuerza el rechazo de la hipotesis nula.

shapiro.test(posibles_p_gorro)

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.92495, p-value = 4.205e-15

La hipotesis nula es que se asemejan a una distribución normal, con un nivel de certeza de 0.95, es decir el margen de error es de 0.05, se puede ver que el p value es menor a este, por lo tanto se rechaza la hipotesis de que se asemeja a una función normal

para n=10, proporción población plantas enfermas 50%, sanas 50%

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

require(car)
par(mfrow=c(1,3))
hist(posibles_p_gorro, xlab = "Porcentaje de plantas enfermas", ylab = "Frecuencia", las=1, main = "", col = "gray")
plot(density(posibles_p_gorro), xlab = "Porcentaje de plantas enfermas", ylab = "Densidad", las=1, main = "")
qqPlot(posibles_p_gorro, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales", las=1,main="")

## [1]  46 193

Graficamente se puede ver que se rechaza la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia.

Valor medio y desviación estandar

mean(posibles_p_gorro)

## [1] 0.506

sd(posibles_p_gorro)

## [1] 0.1594203

Si bien se acerca el valor medio a 0,5 que es el valor real la desviación estandar de los datos es muy alta, lo cual refuerza el rechazo de la hipotesis nula.

shapiro.test(posibles_p_gorro)

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.96442, p-value = 1.21e-09

Comportamiento similar al escenario anterior, se rechaza la hipotesis nula, con el test shapiro

para n=15, proporción población plantas enfermas 50%, sanas 50%

## [1]  34 314

Graficamente se puede ver que se rechaza la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia.

Valor medio y desviación estandar

## [1] 0.5004

## [1] 0.1297783

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.97435, p-value = 1.108e-07

omportamiento similar al escenario anterior, se rechaza la hipotesis nula, con el test shapiro

para n=20, proporción población plantas enfermas 50%, sanas 50%

## [1] 178 291

Graficamente se puede ver que se rechaza la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia

Valor medio y desviación estandar

## [1] 0.4997

## [1] 0.1103143

Si bien se acerca el valor medio a 0,5 que es el valor real la desviación estandar de los datos es muy alta, lo cual refuerza el rechazo de la hipotesis nula. Cade vez que se aumenta el tamaño de la muestra se disminuye la desviación estandar.

para n=30, proporción población plantas enfermas 50%, sanas 50%

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

require(car)
par(mfrow=c(1,3))
hist(posibles_p_gorro, xlab = "Porcentaje de plantas enfermas", ylab = "Frecuencia", las=1, main = "", col = "gray")
plot(density(posibles_p_gorro), xlab = "Porcentaje de plantas enfermas", ylab = "Densidad", las=1, main = "")
qqPlot(posibles_p_gorro, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales", las=1,main="")

## [1] 127 290

Graficamente se puede ver que se rechaza la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.98658, p-value = 0.0001451

omportamiento similar al escenario anterior, se rechaza la hipotesis nula, con el test shapiro

para n=50, proporción población plantas enfermas 50%, sanas 50%

## [1] 177 208

Graficamente se puede ver que se rechaza la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.98745, p-value = 0.0002673

omportamiento similar al escenario anterior, se rechaza la hipotesis nula, con el test shapiro

para n=100, proporción población plantas enfermas 50%, sanas 50%

## [1]  68 373

Graficamente se puede ver que se rechaza la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.99324, p-value = 0.02423

Comportamiento similar al escenario anterior, se rechaza la hipotesis nula, con el test shapiro. Aunque cada vez que aumentamos el tamaño de la muestra nos acercamos a que sea superior a nuestro margen de error de 0,05, par poder aceptar la hipotesis.

para n=200, proporción población plantas enfermas 50%, sanas 50%

## [1] 383 441

Graficamente se puede ver que NO SE RECHAZA la hipotesis nula, la grafica q-q sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.99353, p-value = 0.03077

En esta ocasión pero el método de Shapiro podemos ver que el p Value es mayor a 0,05, por lo tanto NO PODEMOS rechazar la hipotesis nula que la población de plantas enfermas sigue una distribución normal.

para n=500, proporción población plantas enfermas 50%, sanas 50%

## [1] 428  42

Graficamente se puede ver que NO SE RECHAZA la hipotesis nula, la grafica q-q sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.99522, p-value = 0.1265

Conclusión

En esta ocasión por el método de Shapiro podemos ver que el p Value es mayor a 0,05, por lo tanto NO PODEMOS rechazar la hipotesis nula que la población de plantas enfermas sigue una distribución normal.

Para un tamaño de muestras n=500 y n=200 tanto gráficamente como numéricamente se puede llegar a comprobar que el porcentaje de plantas enfermas sigue el comportamiento de una distribución normal.

Vemos como coincide también con un valor medio cercanos al valor real y una deviación estándar menor.

## [1] 0.50022

## [1] 0.01612401

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

##Lote 10 % Enfermas

lote=c(rep("enfermo",100),rep("sana",900))

calc_p_gorro=function(n){

muestra=sample(lote,size = n)

# porcentaje de plantas enfermas
p_gorro=sum(muestra=="enfermo")/n

return(p_gorro)

}

para n=10, proporción población plantas enfermas 10%, sanas 90%

H0: el porcentaje de plantas enfermas sigue una distribución normal.

H1: el porcentaje de plantas enfermas no sigue una distribución normal.

Nivel de significancia: 0.05 (Hipotético).

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

require(car)
par(mfrow=c(1,3))
hist(posibles_p_gorro, xlab = "Porcentaje de plantas enfermas", ylab = "Frecuencia", las=1, main = "", col = "gray")
plot(density(posibles_p_gorro), xlab = "Porcentaje de plantas enfermas", ylab = "Densidad", las=1, main = "")
qqPlot(posibles_p_gorro, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales", las=1,main="")

## [1] 327 498

Graficamente se puede ver que SE RECHAZA la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia. Incluso para el mismo tamaño de muestra la desciación a una distribución normal es mucho mayor que cuando la proporción de enfermas era del 50%.

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.8279, p-value < 2.2e-16

La hipotesis nula es que se asemejan a una distribución normal, con un nivel de certeza de 0.95, es decir el margen de error es de 0.05, se puede ver que el p value es menor a este, por lo tanto se rechaza la hipotesis de que se asemeja a una función normal

para n=15, proporción población plantas enfermas 10%, sanas 90%

## [1]  82 434

Graficamente se puede ver que SE RECHAZA la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia. Incluso para el mismo tamaño de muestra la desciación a una distribución normal es mucho mayor que cuando la proporción de enfermas era del 50%.

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.89337, p-value < 2.2e-16

dado el valor p, se rechaza la hipotesis nula, en mucho menor a 0.05

para n=20, proporción población plantas enfermas 10%, sanas 90%

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

require(car)
par(mfrow=c(1,3))
hist(posibles_p_gorro, xlab = "Porcentaje de plantas enfermas", ylab = "Frecuencia", las=1, main = "", col = "gray")
plot(density(posibles_p_gorro), xlab = "Porcentaje de plantas enfermas", ylab = "Densidad", las=1, main = "")
qqPlot(posibles_p_gorro, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales", las=1,main="")

## [1] 29 84

Graficamente se puede ver que SE RECHAZA la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia. Incluso para el mismo tamaño de muestra la desciación a una distribución normal es mucho mayor que cuando la proporción de enfermas era del 50%.

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.93502, p-value = 6.087e-14

La hipotesis nula es que se asemejan a una distribución normal, con un nivel de certeza de 0.95, es decir el margen de error es de 0.05, se puede ver que el p value es menor a este, por lo tanto se rechaza la hipotesis de que se asemeja a una función normal

para n=30, proporción población plantas enfermas 10%, sanas 90%

## [1] 458  90

Graficamente se puede ver que SE RECHAZA la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia. Incluso para el mismo tamaño de muestra la desciación a una distribución normal es mucho mayor que cuando la proporción de enfermas era del 50%.

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.94787, p-value = 2.889e-12

dado el valor p, se rechaza la hipotesis nula, en mucho menor a 0.05

para n=50, proporción población plantas enfermas 10%, sanas 90%

## [1] 223  44

Graficamente se puede ver que SE RECHAZA la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia. Incluso para el mismo tamaño de muestra la desciación a una distribución normal es mucho mayor que cuando la proporción de enfermas era del 50%.

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.97432, p-value = 1.088e-07

dado el valor p, se rechaza la hipotesis nula, en mucho menor a 0.05

para n=60, proporción población plantas enfermas 10%, sanas 90%

## [1]  28 204

Graficamente se puede ver que SE RECHAZA la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia. Incluso para el mismo tamaño de muestra la desciación a una distribución normal es mucho mayor que cuando la proporción de enfermas era del 50%.

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.97752, p-value = 5.782e-07

dado el valor p, se rechaza la hipotesis nula, en mucho menor a 0.05

para n=60, proporción población plantas enfermas 10%, sanas 90%

## [1] 366  50

Graficamente se puede ver que SE RECHAZA la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia. Incluso para el mismo tamaño de muestra la desciación a una distribución normal es mucho mayor que cuando la proporción de enfermas era del 50%.

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.97534, p-value = 1.833e-07

dado el valor p, se rechaza la hipotesis nula, en mucho menor a 0.05

para n=100, proporción población plantas enfermas 10%, sanas 90%

## [1] 210 335

Graficamente se puede ver que SE RECHAZA la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia. Incluso para el mismo tamaño de muestra la desciación a una distribución normal es mucho mayor que cuando la proporción de enfermas era del 50%.

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.98306, p-value = 1.435e-05

dado el valor p, se rechaza la hipotesis nula, en mucho menor a 0.05

para n=200, proporción población plantas enfermas 10%, sanas 90%

## [1] 169 435

Graficamente se puede ver que SE RECHAZA la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia. Incluso para el mismo tamaño de muestra la desciación a una distribución normal es mucho mayor que cuando la proporción de enfermas era del 50%.

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.99191, p-value = 0.008047

dado el valor p, se rechaza la hipotesis nula, en mucho menor a 0.05.

Conclusión

Lo interesante es que en el caso en que la proporción enfermos-sanos era 50%-50%, para una muestra n=200 nos dio como resultado que no se rechazaba la hipotesis, pero cuando la proporción de enfermos es menor 10%, para el mismo tamaño de muesra se rechaza la hipotesis. Lo que se puede concluir es que cuando la proporción que estamos buscando caracterizar es muy pequeña en la población requeriremos tamaños mayores de muestras para poder tener una caracterización adecuada.

para n=500, proporción población plantas enfermas 10%, sanas 90%

## [1] 488 366

Graficamente se puede ver que SE RECHAZA la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, sigue la linea de referencia.

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.99434, p-value = 0.06067

Se rechaza la hipotesis nula,esto nos permite concluir que para tener una caracterización adecuada cuando la proporción de plantas enfermas es pequeña, requeriremos muestras cada vez mas grandes.

##Lote 90 % Enfermas

lote=c(rep("enfermo",900),rep("sana",100))

calc_p_gorro=function(n){

muestra=sample(lote,size = n)

# porcentaje de plantas enfermas
p_gorro=sum(muestra=="enfermo")/n

return(p_gorro)

}

para n=10, proporción población plantas enfermas 90%, sanas 10%

H0: el porcentaje de plantas enfermas sigue una distribución normal.

H1: el porcentaje de plantas enfermas no sigue una distribución normal.

Nivel de significancia: 0.05 (Hipotético).

## [1] 172 161

Graficamente se puede ver que SE RECHAZA la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia.

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.83791, p-value < 2.2e-16

La hipotesis nula es que se asemejan a una distribución normal, con un nivel de certeza de 0.95, es decir el margen de error es de 0.05, se puede ver que el p value es menor a este, por lo tanto se rechaza la hipotesis de que se asemeja a una función normal

para n=20, proporción población plantas enfermas 90%, sanas 10%

## [1] 242 345

Graficamente se puede ver que SE RECHAZA la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia.

shapiro.test(posibles_p_gorro)

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.90738, p-value < 2.2e-16

por el valor inferior a 0,05 de p, se rechaza la hipotesis.

para n=30, proporción población plantas enfermas 90%, sanas 10%

## [1]  95 247

Graficamente se puede ver que SE RECHAZA la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia.

shapiro.test(posibles_p_gorro)

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.9525, p-value = 1.354e-11

Se rechaza la hipotesis

para n=50, proporción población plantas enfermas 90%, sanas 10%

## [1] 18 45

Graficamente se puede ver que SE RECHAZA la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia.

shapiro.test(posibles_p_gorro)

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.97594, p-value = 2.492e-07

Se rechaza la hipotesis

para n=100, proporción población plantas enfermas 90%, sanas 10%

## [1] 223 230

Graficamente se puede ver que SE RECHAZA la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia.

shapiro.test(posibles_p_gorro)

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.98657, p-value = 0.0001444

Se rechaza la hipotesis

para n=200, proporción población plantas enfermas 90%, sanas 10%

## [1]  32 206

Graficamente se puede ver que SE RECHAZA la hipotesis nula, la grafica q-q no sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia.

shapiro.test(posibles_p_gorro)

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.98563, p-value = 7.609e-05

Se rechaza la hipotesis

para n=500, proporción población plantas enfermas 90%, sanas 10%

## [1] 133  13

Graficamente se puede ver que NO SE RECHAZA la hipotesis nula, la grafica q-q sigue el comportamiento de una distribución normal ideal, no sigue la linea de referencia.

shapiro.test(posibles_p_gorro)

## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro
## W = 0.99128, p-value = 0.004842

Conslusión

Pero númericamente SE RECHAZA la hipotesis, se puede ver que en los casos extremos donde las proporciones estan muy desvalanceadas caso (10% enfermas y 90% enfermas), se requiere muestras mas grande para poder obtener una caracterización adecuada.

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.

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

#Lote de población N1 y N2
loteN1 = c(rep("Enfermo", 100), rep("Sana",900))

loteN2 = c(rep("Enfermo", 150), rep("Sana",1350))

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.

##Generar muestras de lote1 (n1) y lote 2 (n2)

calc_dif_p=function(n1){

n2=n1

muestra1=sample(loteN1,n1)
p1=sum(muestra1=="Enfermo")/n1

muestra2=sample(loteN2,n2)
p2=sum(muestra2=="Enfermo")/n2

dif_p=p1-p2
return(dif_p)
}

ejemplo:

calc_dif_p(n1=100)

## [1] -0.04

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?

## para una tamaño de muestra n1=n2=300
dif_p=sapply(rep(300,500), calc_dif_p)
table(dif_p==0)

## 
## FALSE  TRUE 
##   467    33

summary(dif_p)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -0.073333 -0.013333  0.003333  0.001687  0.016667  0.066667

hist(dif_p)
abline(v=mean(dif_p), col="red", lwd=3)

Cuando representamos una distribución podemos analizar su nivel de simetría: una distribución es simétrica si en relación a un valor central la distribución se distribuye un 50% a la derecha y otro 50% a la izquierda, presentando una forma similar a ambos lados del valor central.

Podemos decir que no es perfectamente simetrica dado que la media y mediana no tienen el mismo valor. Como se indica en la tabla y en la distribución las diferencias no siempre son cero, aunque la realidad de la población si lo es.

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?

simulación para n1=n2=5, proporción de efectividad tratamiento 10%/ 10%

dif_p=sapply(rep(5,500), calc_dif_p)
table(dif_p==0)

## 
## FALSE  TRUE 
##   256   244

require(car)
par(mfrow=c(1,3))
hist(dif_p, xlab = " diferencia Porcentaje % desempeño ", ylab = "Frecuencia", las=1, main = "", col = "gray")
abline(v=mean(dif_p), col="red", lwd=3)
plot(density(dif_p), xlab = " diferencia Porcentaje % desempeño ", ylab = "Densidad", las=1, main = "")
qqPlot(dif_p, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales", las=1,main="")

## [1] 286  15

evaluar normalidad

H0: la diferencias de los porcentajes de desempeño sigue una distribución normal

H1: la diferencias de los porcentajes de desempeño sigue una distribución normal

Nivel de significancia: 0.05 (Hipotético).

shapiro.test(dif_p)

## 
##  Shapiro-Wilk normality test
## 
## data:  dif_p
## W = 0.89167, p-value < 2.2e-16

simulación para n1=n2=10, proporción de efectividad tratamiento 10%/ 10%

dif_p=sapply(rep(10,500), calc_dif_p)
table(dif_p==0)

## 
## FALSE  TRUE 
##   339   161

require(car)
par(mfrow=c(1,3))
hist(dif_p, xlab = " diferencia Porcentaje % desempeño ", ylab = "Frecuencia", las=1, main = "", col = "gray")
abline(v=mean(dif_p), col="red", lwd=3)
plot(density(dif_p), xlab = " diferencia Porcentaje % desempeño ", ylab = "Densidad", las=1, main = "")
qqPlot(dif_p, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales", las=1,main="")

## [1] 422  71

evaluar normalidad

H0: la diferencias de los porcentajes de desempeño sigue una distribución normal

H1: la diferencias de los porcentajes de desempeño sigue una distribución normal

Nivel de significancia: 0.05 (Hipotético).

shapiro.test(dif_p)

## 
##  Shapiro-Wilk normality test
## 
## data:  dif_p
## W = 0.94946, p-value = 4.847e-12

Se rechaza la hipotesis

simulación para n1=n2=15, proporción de efectividad tratamiento 10%/ 10%

dif_p=sapply(rep(15,500), calc_dif_p)
table(dif_p==0)

## 
## FALSE  TRUE 
##   353   147

require(car)
par(mfrow=c(1,3))
hist(dif_p, xlab = " diferencia Porcentaje % desempeño ", ylab = "Frecuencia", las=1, main = "", col = "gray")
abline(v=mean(dif_p), col="red", lwd=3)
plot(density(dif_p), xlab = " diferencia Porcentaje % desempeño ", ylab = "Densidad", las=1, main = "")
qqPlot(dif_p, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales", las=1,main="")

## [1] 211 125

evaluar normalidad

H0: la diferencias de los porcentajes de desempeño sigue una distribución normal

H1: la diferencias de los porcentajes de desempeño sigue una distribución normal

Nivel de significancia: 0.05 (Hipotético).

shapiro.test(dif_p)

## 
##  Shapiro-Wilk normality test
## 
## data:  dif_p
## W = 0.96313, p-value = 7.113e-10

Se rechaza la hipotesis

simulación para n1=n2=20, proporción de efectividad tratamiento 10%/ 10%

## 
## FALSE  TRUE 
##   400   100

## [1] 49 72

evaluar normalidad

H0: la diferencias de los porcentajes de desempeño sigue una distribución normal

H1: la diferencias de los porcentajes de desempeño sigue una distribución normal

Nivel de significancia: 0.05 (Hipotético).

## 
##  Shapiro-Wilk normality test
## 
## data:  dif_p
## W = 0.97545, p-value = 1.938e-07

Se rechaza la hipotesis

simulación para n1=n2=30, proporción de efectividad tratamiento 10%/ 10%

dif_p=sapply(rep(30,500), calc_dif_p)
table(dif_p==0)

## 
## FALSE  TRUE 
##   402    98

require(car)
par(mfrow=c(1,3))
hist(dif_p, xlab = " diferencia Porcentaje % desempeño ", ylab = "Frecuencia", las=1, main = "", col = "gray")
abline(v=mean(dif_p), col="red", lwd=3)
plot(density(dif_p), xlab = " diferencia Porcentaje % desempeño ", ylab = "Densidad", las=1, main = "")
qqPlot(dif_p, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales", las=1,main="")

## [1] 226 356

evaluar normalidad

H0: la diferencias de los porcentajes de desempeño sigue una distribución normal

H1: la diferencias de los porcentajes de desempeño sigue una distribución normal

Nivel de significancia: 0.05 (Hipotético).

shapiro.test(dif_p)

## 
##  Shapiro-Wilk normality test
## 
## data:  dif_p
## W = 0.97893, p-value = 1.253e-06

Se rechaza la hipotesis

simulación para n1=n2=50, proporción de efectividad tratamiento 10%/ 10%

## 
## FALSE  TRUE 
##   452    48

## [1] 196 145

evaluar normalidad

H0: la diferencias de los porcentajes de desempeño sigue una distribución normal

H1: la diferencias de los porcentajes de desempeño sigue una distribución normal

Nivel de significancia: 0.05 (Hipotético).

## 
##  Shapiro-Wilk normality test
## 
## data:  dif_p
## W = 0.9858, p-value = 8.528e-05

Se rechaza la hipotesis

Pero es claro que a medida que aumentamos el tamaño de la muestra nos acercamos a no rechazar la hipotesis.

simulación para n1=n2=100, proporción de efectividad tratamiento 10%/ 10%

## 
## FALSE  TRUE 
##   457    43

## [1] 206  61

evaluar normalidad

H0: la diferencias de los porcentajes de desempeño sigue una distribución normal

H1: la diferencias de los porcentajes de desempeño sigue una distribución normal

Nivel de significancia: 0.05 (Hipotético).

## 
##  Shapiro-Wilk normality test
## 
## data:  dif_p
## W = 0.99279, p-value = 0.01665

Se rechaza la hipotesis

simulación para n1=n2=200, proporción de efectividad tratamiento 10%/ 10%

## 
## FALSE  TRUE 
##   468    32

## [1] 234 122

Graficamente podriamos afirmar que sigue una distribución normal, y para muestras n1=n2 superiores a 200 se esta logrando simetria en los distribución con un margen de error pequeño.

evaluar normalidad

H0: la diferencias de los porcentajes de desempeño sigue una distribución normal

H1: la diferencias de los porcentajes de desempeño sigue una distribución normal

Nivel de significancia: 0.05 (Hipotético).

## 
##  Shapiro-Wilk normality test
## 
## data:  dif_p
## W = 0.99425, p-value = 0.05617

NO se rechaza la hipotesis

simulación para n1=n2=500, proporción de efectividad tratamiento 10%/ 10%

## 
## FALSE  TRUE 
##   475    25

## [1] 441  56

Graficamente podriamos afirmar que sigue una distribución normal, y para muestras n1=n2 superiores a 200 se esta logrando simetria en los distribución con un margen de error pequeño.

evaluar normalidad

H0: la diferencias de los porcentajes de desempeño sigue una distribución normal

H1: la diferencias de los porcentajes de desempeño sigue una distribución normal

Nivel de significancia: 0.05 (Hipotético).

## 
##  Shapiro-Wilk normality test
## 
## data:  dif_p
## W = 0.99611, p-value = 0.2591

No se rechaza la hipotesis

##Conclusión

Se puede observar que a medida que aumentamos los tamaños de las muestras en las mismas proporción en las dos poblaciones de plantas con tratamientos, se mejora la simetría en los resultados, con una media muy ajustada a 0% que es la realidad, además los datos tienen poca dispersión lo cual nos da señales que para muestras mayores o iguales a 200 plantas podríamos lograr buena inferencia sobre la población.

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%)?

#Lote de población N1 y N2
loteN1 = c(rep("Enfermo", 100), rep("Sana",900))

loteN2 = c(rep("Enfermo", 225), rep("Sana",1275))

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.

##Generar muestras de lote1 (n1) y lote 2 (n2)

calc_dif_p=function(n1){

n2=n1

muestra1=sample(loteN1,n1)
p1=sum(muestra1=="Enfermo")/n1

muestra2=sample(loteN2,n2)
p2=sum(muestra2=="Enfermo")/n2

dif_p=p1-p2
return(dif_p)
}

simulación para n1=n2=15, proporción de efectividad tratamiento 10%/ 15%

dif_p=sapply(rep(15,500), calc_dif_p)
table(dif_p==0)

## 
## FALSE  TRUE 
##   403    97

require(car)
par(mfrow=c(1,3))
hist(dif_p, xlab = " diferencia Porcentaje % desempeño ", ylab = "Frecuencia", las=1, main = "", col = "gray")
abline(v=mean(dif_p), col="red", lwd=3)
plot(density(dif_p), xlab = " diferencia Porcentaje % desempeño ", ylab = "Densidad", las=1, main = "")
qqPlot(dif_p, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales", las=1,main="")

## [1]  45 244

## 
##  Shapiro-Wilk normality test
## 
## data:  dif_p
## W = 0.97011, p-value = 1.446e-08

Se rechaza la hipotesis graficamente y númericamente.

simulación para n1=n2=30, proporción de efectividad tratamiento 10%/ 15%

## 
## FALSE  TRUE 
##   439    61

## [1] 114 206

## 
##  Shapiro-Wilk normality test
## 
## data:  dif_p
## W = 0.98696, p-value = 0.00019

Se rechaza la hipotesis graficamente y númericamente.

simulación para n1=n2=50, proporción de efectividad tratamiento 10%/ 15%

## 
## FALSE  TRUE 
##   451    49

## [1]  94 268

## 
##  Shapiro-Wilk normality test
## 
## data:  dif_p
## W = 0.98677, p-value = 0.0001666

Se rechaza la hipotesis graficamente y númericamente.

simulación para n1=n2=50, proporción de efectividad tratamiento 10%/ 15%

## 
## FALSE  TRUE 
##   478    22

## [1] 412 279

## 
##  Shapiro-Wilk normality test
## 
## data:  dif_p
## W = 0.99241, p-value = 0.01209

Se rechaza la hipotesis graficamente y númericamente.

simulación para n1=n2=200, proporción de efectividad tratamiento 10%/ 15%

## 
## FALSE  TRUE 
##   494     6

## [1] 381 336

## 
##  Shapiro-Wilk normality test
## 
## data:  dif_p
## W = 0.99394, p-value = 0.04333

Se rechaza la hipotesis graficamente y númericamente.

simulación para n1=n2=50, proporción de efectividad tratamiento 10%/ 15%

## 
## FALSE 
##   500

## [1] 106 160

## 
##  Shapiro-Wilk normality test
## 
## data:  dif_p
## W = 0.99639, p-value = 0.3194

Se rechaza la hipotesis graficamente y númericamente.

Conclusión

Cuando se presenta diferencia en las proporciones de efectividad de los tratamientos en cada población, requerimos muestras mas grandes para lograr una inferencia mejor, para el caso de proporciones 10% -15 %, requeriamos muestras mayores a 500 plantas, caso diferente cuando la proporciones eran iguales 10%-10% donde era suficiente desde una muestra mayor a 200 plantas.

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.

Resumen: Matt Motyl un estudiante de doctorado dentro de su investigación planteo que los políticos de izquierda o de derecha tienden a ser extremistas y ver el mundo en tonos de negro o blanco, al contrario, los de centro tienen den a ver o analizar el mundo en tonos grises, su estudio se baso en los resultados de 2000 personas. Su hipotesis tomaba fuerza pues los estudios numéricos del valor p eran 0.01, lo cual se considera muy significativo. Lo curioso es que el estudio se repitió con información extra y el valor P resulto ser de 0,59, lejos del valor de significancia [1]. ¿qué sucedió? ¿cambio la situación analizada? ¿fueron manipulados mal los datos? Lo cierto es que no tiene que ver con los datos o con el análisis matemático. Es la interpretación que se le da al valor P como certero e incuestionable.

La verdad es que la hipotesis de Motyl sobre los políticos de izquierda o derecha, nos aplica a todos, cuando las personas buscamos respuestas a algo, queremos respuestas del tipo (SI, NO), (Se rechaza, no se rechaza), (Pertences , no pertenece) y así muchos ejemplos en diferentes contextos donde buscamos respuestas, buscamos un indicador mágico que nos ayude a dar claridad de los problemas con respuestas binarias, no nos gusta los matices. En esto radica la tergiversación que se le ha dado al valor P como el “estándar de oro” en la valides estadística. En parte como se indica en [1], el valor P a tenido muchas criticas y se le considera como los mosquitos molestos e imposibles de aplastar, considero en parte por la tendencia de las personas a buscar indicadores incuestionables y sin matices.

Cuando el estadístico británico Ronald Fisher introdujo el valor P en la década de 1920, no pretendía que fuera una prueba definitiva. Lo pensó simplemente como una forma informal de juzgar si la evidencia era significativa pero sujeto a revisión que involucra combinación de datos y conocimiento previo, para finalmente si llegar a una conclusión científica. En muchos estudios científicos se ha utilizado únicamente el valor P para validar hipotesis, lo cual esta causando que muchos estudios e inferencias realizadas requieran ser revaluadas. No es sorpresa que la la Asociación Estadounidense de Estadística (ASA, por sus siglas en inglés) aconseja a los investigadores que eviten sacar conclusiones científicas o tomar decisiones políticas basándose únicamente en los valores P [2]. Retomando el ejemplo de Motyl cuando obtuvo el resultado de P igual a 0,01 no significaba que pudiera un 1% de probabilidad de que su resultado fuera una falsa alarma, El valor P no puede decir esto: todo lo que puede hacer es resumir los datos asumiendo una hipótesis nula específica, no puede hacer afirmaciones sobre la realidad.

Considero que es importante tener conciencia del contexto de donde se toman los datos, realizar exploración de datos y evaluar que tan plausible es la hipotesis que voy a validar, el valor P no nos dará la plausibilidad de la hipotesis. Una buena exploración de datos, conocimiento del contexto y una hipotesis plausible nos darán buenos resultados de inferencia cuando utilizamos el valor P para validar nuestra hipotesis. Lo que en realidad obtuvo Motyl con su valor de P de 0,01 es que tenia un chance de tener hasta un 11% de falsos positivos, eso asumiendo que efectivamente su hipotesis tenía un efecto real.

La conciencia sobre utilizar el valor P como única herramienta de validación de hipotesis viene creciendo, aunque no al ritmo deseado, en parte por la cultura arraigada de los científicos de buscar respuestas con un solo indicador, en parte influye como se viene enseñando estadística, seguramente con el interés cada vez mas grande en explorar información oculta en gran cantidad de datos, con las diferentes herramientas que nos pone al servicio la ciencia de datos, se logrará el objetivo de transformar el valor de P como la única herramienta de valoración de valides.

Referencias: [1] “Statistical Errors: P values, the gold standard of statistical validity, are not as reliable as many scientists assume” [2] “Statisticians issue warning on P values: Statement aims to halt missteps in the quest for certainty”

Taller unidad 2

ELKIN CANTOR

23/3/2022

Inferencia estadística y simulación

Punto 1

Conclusión

para n=5, proporción población plantas enfermas 50%,50%

para n=10, proporción población plantas enfermas 50%, sanas 50%

para n=15, proporción población plantas enfermas 50%, sanas 50%

para n=20, proporción población plantas enfermas 50%, sanas 50%

para n=30, proporción población plantas enfermas 50%, sanas 50%

para n=50, proporción población plantas enfermas 50%, sanas 50%

para n=100, proporción población plantas enfermas 50%, sanas 50%

para n=200, proporción población plantas enfermas 50%, sanas 50%

para n=500, proporción población plantas enfermas 50%, sanas 50%

Conclusión

para n=10, proporción población plantas enfermas 10%, sanas 90%

para n=15, proporción población plantas enfermas 10%, sanas 90%

para n=20, proporción población plantas enfermas 10%, sanas 90%

para n=30, proporción población plantas enfermas 10%, sanas 90%

para n=50, proporción población plantas enfermas 10%, sanas 90%

para n=60, proporción población plantas enfermas 10%, sanas 90%

para n=60, proporción población plantas enfermas 10%, sanas 90%

para n=100, proporción población plantas enfermas 10%, sanas 90%

para n=200, proporción población plantas enfermas 10%, sanas 90%

Conclusión

para n=500, proporción población plantas enfermas 10%, sanas 90%

para n=10, proporción población plantas enfermas 90%, sanas 10%

para n=20, proporción población plantas enfermas 90%, sanas 10%

para n=30, proporción población plantas enfermas 90%, sanas 10%

para n=50, proporción población plantas enfermas 90%, sanas 10%

para n=100, proporción población plantas enfermas 90%, sanas 10%

para n=200, proporción población plantas enfermas 90%, sanas 10%

para n=500, proporción población plantas enfermas 90%, sanas 10%

Conslusión

Punto 2

simulación para n1=n2=5, proporción de efectividad tratamiento 10%/ 10%

simulación para n1=n2=10, proporción de efectividad tratamiento 10%/ 10%

simulación para n1=n2=15, proporción de efectividad tratamiento 10%/ 10%

simulación para n1=n2=20, proporción de efectividad tratamiento 10%/ 10%

simulación para n1=n2=30, proporción de efectividad tratamiento 10%/ 10%

simulación para n1=n2=50, proporción de efectividad tratamiento 10%/ 10%

simulación para n1=n2=100, proporción de efectividad tratamiento 10%/ 10%

simulación para n1=n2=200, proporción de efectividad tratamiento 10%/ 10%

simulación para n1=n2=500, proporción de efectividad tratamiento 10%/ 10%

simulación para n1=n2=15, proporción de efectividad tratamiento 10%/ 15%

simulación para n1=n2=30, proporción de efectividad tratamiento 10%/ 15%

simulación para n1=n2=50, proporción de efectividad tratamiento 10%/ 15%

simulación para n1=n2=50, proporción de efectividad tratamiento 10%/ 15%

simulación para n1=n2=200, proporción de efectividad tratamiento 10%/ 15%

simulación para n1=n2=50, proporción de efectividad tratamiento 10%/ 15%

Conclusión

Punto 3