Unidad 2 - Taller de simulación en R
Universidad Pontificia Javeriana Cali
Asignatura métodos y simulación estadística
Catalina Gómez Vallejo
1. Teorema del Límite Central
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%.
#### A. Parametro = P = 500/1000 = 0.50 = 50%
set.seed(0)
Lote<-c(rep("Enfermas",500),rep("Sanas",500))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.
Se empleo un tamaño de muestra n=100.
set.seed(0)
plantas_enfermas=function(n){
muestra=sample(Lote,size = n)
return(sum(muestra=="Enfermas")/n)
}
plantas_enfermas(n = 100)## [1] 0.46c). Repita el escenario anterior (b) 500 veces y analice los resultados en cuanto al comportamiento de los 500 estimadores.
set.seed(0)
plantas_enfermas=function(n){
muestra=sample(Lote,size = n)
return(sum(muestra=="Enfermas")/n)
}
Estimadores=sapply(rep(100,500), plantas_enfermas)¿Qué tan simétricos son los datos?
library("moments")
skewness(Estimadores)## [1] 0.01749204kurtosis(Estimadores)## [1] 3.140607par(mfrow=c(1,2))
hist(Estimadores,main="Histograma de estimadores P=0.50",ylab="Frecuencia",col="pink")
boxplot(Estimadores,main="Boxplot de estimadores")
Dado que el indicador de asimetría es de 0.01749204 y la curtosis 3.140607, se puede considerar que la distribución de los datos es lepticúrtica y presenta una leve asimetría positiva. En los graficos se puede observar que los datos se concentran alrededor del valor del parámetro (0.50).
¿Son sesgados y qué pasa en cuanto a variabilidad?
Sesgo=mean(Estimadores)-0.50
Sesgo## [1] 0.00114sd(Estimadores)## [1] 0.04614281El sesgo del estimador es de 0.00114, y la desviación estandar es baja por lo cual se puede considerar que los datos se concentran alrededor de la media.
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).
# Literal B.
set.seed(0)
plantas_enfermas=function(n){
muestra=sample(Lote,size = n)
return(sum(muestra=="Enfermas")/n)
}
S<-sapply(c(5,10,15,20,30,50,60,100,200,500), plantas_enfermas)
data.frame(Tamaño_muestra_n=c(5,10,15,20,30,50,60,100,200,500),Estimador= S)## Tamaño_muestra_n Estimador
## 1 5 0.4000000
## 2 10 0.6000000
## 3 15 0.3333333
## 4 20 0.3500000
## 5 30 0.6333333
## 6 50 0.4400000
## 7 60 0.5333333
## 8 100 0.4600000
## 9 200 0.5150000
## 10 500 0.4960000Se puede observar que a mayor tamaño de muestra (n) el estimador tiende a aproximarse con mayor exactitud al valor del parámetro.
# Literal C.
set.seed(0)
plantas_enfermas=function(n){
muestra=sample(Lote,size = n)
return(sum(muestra=="Enfermas")/n)
}
ntamaños=c(5,10,15,20,30,50,60,100,200,500)
estimadores<-data.frame(matrix(NA,nrow=500,ncol=length(ntamaños)))
for (i in 1:10) {
estimadores[,i]=sapply(rep(ntamaños[i],500), plantas_enfermas)
}
# 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)
St<-c()
for (i in 1:10) {
St[i]<-shapiro.test(estimadores[,i])$p.value
}
data.frame(Tamaño_muestra_n=c(5,10,15,20,30,50,60,100,200,500),Valor_P=St)## Tamaño_muestra_n Valor_P
## 1 5 4.804487e-15
## 2 10 1.196627e-09
## 3 15 3.075849e-08
## 4 20 7.170556e-06
## 5 30 1.802550e-04
## 6 50 5.310096e-04
## 7 60 4.993082e-04
## 8 100 5.015948e-02
## 9 200 1.409295e-01
## 10 500 9.194540e-02A partir del tamaño de muestra n=100, se puede apreciar que el Valor P de la prueba de normalidad shapiro wilks es mayor al nivel de significancia 0,05 por lo cual, se puede considerar que los estimadores obtenidos en la simulación tienen una distribución normal cuando en n es mayor o igual a 100.
par(mfrow=c(2, 5))
v<-c(5,10,15,20,30,50,60,100,200,500)
for (l in v) {
for (i in 1:10) {
l=v[i]
qqnorm(estimadores[,i], main = bquote(~ Tamaño_muestra_n== .(l)))
qqline(estimadores[,i],col="red")
}
}
En el grafico anterior se puede apreciar que a menor tamaño de muestra los estimadores no se ajustan a la linea roja, es decir, parecen no estar distribuidos normalmente, mientras que a mayor tamaño de muestra los puntos de los estimadores se ajustan a la recta (Son normales).
par(mfrow=c(2, 5))
v<-c(5,10,15,20,30,50,60,100,200,500)
for (l in v) {
for (i in 1:10) {
l=v[i]
hist(estimadores[,i],main = bquote(~ Tamaño_muestra_n== .(l)), xlab = "Estimadores")
}
}
En el grafico anterior se puede apreciar que a mayor tamaño de muestra los estimadores parecen comportarse como una distribución normal, y se concentran en mayor cantidad alrededor del parámetro (0.50).
e). LOTE 10% PLANTAS ENFERMAS
# A.
set.seed(0)
Lote2<-c(rep("Enfermas",100),rep("Sanas",900))
# B. Se selecciona una muestra aleatoria de n=100 y se calcula el estimador de la proporción muestral
set.seed(0)
plantas_enfermas2=function(n){
muestra=sample(Lote2,size = n)
return(sum(muestra=="Enfermas")/n)
}
plantas_enfermas2(n = 100)## [1] 0.07# C. Repita el escenario anterior 500 veces y analice los resultados en cuanto al comportamiento de los 500 estimadores.
set.seed(0)
plantas_enfermas2=function(n){
muestra=sample(Lote2,size = n)
return(sum(muestra=="Enfermas")/n)
}
Estimadores2=sapply(rep(100,500), plantas_enfermas2)¿Qué tan simétricos son los datos?
skewness(Estimadores2)## [1] 0.2316654kurtosis(Estimadores2)## [1] 3.363282par(mfrow=c(1,2))
hist(Estimadores2, main="Histograma de estimadores P=0.10",ylab="Frecuencia",col="Blue")
boxplot(Estimadores2, main="Boxplot de estimadores")
Dado que el indicador de asimetría es de 0.2316654 y la curtosis 3.363282, se puede considerar que la distribución de los datos es lepticúrtica y presenta asimetría positiva. En los graficos se puede observar que los datos se concentran alrededor del valor del parámetro (0.10).
¿Son sesgados y qué pasa en cuanto a variabilidad?
Sesgo=mean(Estimadores2)-0.10
Sesgo## [1] -0.00056sd(Estimadores2)## [1] 0.02829997El sesgo del estimador es de -0.00056, y la desviación estandar es baja por lo cual se puede considerar que los datos se concentran alrededor de la media.
# D.
# Literal B.
set.seed(0)
plantas_enfermas2=function(n){
muestra=sample(Lote2,size = n)
return(sum(muestra=="Enfermas")/n)
}
S2<-sapply(c(5,10,15,20,30,50,60,100,200,500), plantas_enfermas2)
data.frame(Tamaño_muestra_n=c(5,10,15,20,30,50,60,100,200,500),Estimador= S2)## Tamaño_muestra_n Estimador
## 1 5 0.00000000
## 2 10 0.00000000
## 3 15 0.06666667
## 4 20 0.00000000
## 5 30 0.13333333
## 6 50 0.06000000
## 7 60 0.11666667
## 8 100 0.09000000
## 9 200 0.10000000
## 10 500 0.10200000Se puede observar que a mayor tamaño de muestra (n) el estimador tiende a aproximarse con mayor exactitud al valor del parámetro (0.10).
# Literal C.
set.seed(0)
plantas_enfermas2=function(n){
muestra=sample(Lote2,size = n)
return(sum(muestra=="Enfermas")/n)
}
ntamaños=c(5,10,15,20,30,50,60,100,200,500)
estimadores2<-data.frame(matrix(NA,nrow=500,ncol=length(ntamaños)))
for (i in 1:10) {
estimadores2[,i]=sapply(rep(ntamaños[i],500), plantas_enfermas2)
}
# 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)
St2<-c()
for (i in 1:10) {
St2[i]<-shapiro.test(estimadores2[,i])$p.value
}
data.frame(Tamaño_muestra_n=c(5,10,15,20,30,50,60,100,200,500),Valor_P=St2)## Tamaño_muestra_n Valor_P
## 1 5 5.747995e-29
## 2 10 3.888739e-22
## 3 15 1.774752e-18
## 4 20 2.894111e-15
## 5 30 1.027430e-14
## 6 50 7.744160e-10
## 7 60 8.875991e-08
## 8 100 8.862828e-06
## 9 200 6.088116e-03
## 10 500 1.446940e-02A partir del tamaño de muestra n=500, se puede apreciar que el Valor P de la prueba de normalidad shapiro wilks es mayor al nivel de significancia 0,05 por lo cual, se puede considerar que los estimadores obtenidos en la simulación tienen una distribución normal cuando en n es mayor o igual a 500.
par(mfrow=c(2, 5))
v<-c(5,10,15,20,30,50,60,100,200,500)
for (l in v) {
for (i in 1:10) {
l=v[i]
qqnorm(estimadores2[,i], main = bquote(~ Tamaño_muestra_n== .(l)))
qqline(estimadores2[,i],col="red")
}
}
En el grafico anterior se puede apreciar que a menor tamaño de muestra
los estimadores no se ajustan a la linea roja, es decir, parecen no
estar distribuidos normalmente, mientras que a mayor tamaño de muestra
los puntos de los estimadores se ajustan a la recta (Son normales).
par(mfrow=c(2, 5))
v<-c(5,10,15,20,30,50,60,100,200,500)
for (l in v) {
for (i in 1:10) {
l=v[i]
hist(estimadores2[,i],main = bquote(~ Tamaño_muestra_n== .(l)), xlab = "Estimadores")
}
}
En el grafico anterior se puede apreciar que a mayor tamaño de muestra los estimadores parecen comportarse como una distribución normal, y se concentran en mayor cantidad alrededor del parámetro (0.10).
e). LOTE 90% PLANTAS ENFERMAS
# A.
set.seed(0)
Lote3<-c(rep("Enfermas",900),rep("Sanas",100))
# B. Se selecciona una muestra aleatoria de n=100 y se calcula el estimador de la proporción muestral
set.seed(0)
plantas_enfermas3=function(n){
muestra=sample(Lote3,size = n)
return(sum(muestra=="Enfermas")/n)
}
plantas_enfermas3(n = 100)## [1] 0.9# C. Repita el escenario anterior 500 veces y analice los resultados en cuanto al comportamiento de los 500 estimadores.
set.seed(0)
plantas_enfermas3=function(n){
muestra=sample(Lote3,size = n)
return(sum(muestra=="Enfermas")/n)
}
Estimadores3=sapply(rep(100,500), plantas_enfermas3)¿Qué tan simétricos son los datos?
skewness(Estimadores3)## [1] -0.1934003kurtosis(Estimadores3)## [1] 3.792par(mfrow=c(1,2))
hist(Estimadores3, main="Histograma de estimadores P=0.90",ylab="Frecuencia",col="Green")
boxplot(Estimadores3, main="Boxplot de estimadores")
Dado que el indicador de asimetría es de -0.1934003 y la curtosis 3.792, se puede considerar que la distribución de los datos es lepticúrtica y presenta una asimetría negativa. En los graficos se puede observar que los datos se concentran alrededor del valor del parámetro (0.90).
¿Son sesgados y qué pasa en cuanto a variabilidad?
Sesgo=mean(Estimadores3)-0.90
Sesgo## [1] 0.00232sd(Estimadores3)## [1] 0.02622199El sesgo del estimador es de 0.00232, y la desviación estandar es baja por lo cual se puede considerar que los datos se concentran alrededor de la media.
# D.
# Literal B.
set.seed(0)
plantas_enfermas3=function(n){
muestra=sample(Lote3,size = n)
return(sum(muestra=="Enfermas")/n)
}
S3<-sapply(c(5,10,15,20,30,50,60,100,200,500), plantas_enfermas3)
data.frame(Tamaño_muestra_n=c(5,10,15,20,30,50,60,100,200,500),Estimador= S3)## Tamaño_muestra_n Estimador
## 1 5 0.8000000
## 2 10 0.8000000
## 3 15 0.9333333
## 4 20 0.8000000
## 5 30 1.0000000
## 6 50 0.9200000
## 7 60 0.9166667
## 8 100 0.9100000
## 9 200 0.8850000
## 10 500 0.8940000Se puede observar que a mayor tamaño de muestra (n) el estimador tiende a aproximarse con mayor exactitud al valor del parámetro (0.90).
# Literal C.
set.seed(0)
plantas_enfermas3=function(n){
muestra=sample(Lote3,size = n)
return(sum(muestra=="Enfermas")/n)
}
ntamaños=c(5,10,15,20,30,50,60,100,200,500)
estimadores3<-data.frame(matrix(NA,nrow=500,ncol=length(ntamaños)))
for (i in 1:10) {
estimadores3[,i]=sapply(rep(ntamaños[i],500), plantas_enfermas3)
}
# 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)
St3<-c()
for (i in 1:10) {
St3[i]<-shapiro.test(estimadores3[,i])$p.value
}
data.frame(Tamaño_muestra_n=c(5,10,15,20,30,50,60,100,200,500),Valor_P=St3)## Tamaño_muestra_n Valor_P
## 1 5 8.704355e-30
## 2 10 1.826071e-22
## 3 15 1.855704e-18
## 4 20 2.187717e-15
## 5 30 1.474995e-12
## 6 50 9.417340e-09
## 7 60 6.179770e-08
## 8 100 7.557976e-06
## 9 200 4.311948e-03
## 10 500 8.063048e-03A partir de los resultados anteriores (Valor P de la prueba de normalidad shapiro wilks) con base en el Valor p se puede concluir que para los diferentes tamaños de muestra planteador los estimadores no siguen una distribución normal.
par(mfrow=c(2, 5))
v<-c(5,10,15,20,30,50,60,100,200,500)
for (l in v) {
for (i in 1:10) {
l=v[i]
qqnorm(estimadores3[,i], main = bquote(~ Tamaño_muestra_n== .(l)))
qqline(estimadores3[,i],col="red")
}
}
par(mfrow=c(2, 5))
v<-c(5,10,15,20,30,50,60,100,200,500)
for (l in v) {
for (i in 1:10) {
l=v[i]
hist(estimadores3[,i],main = bquote(~ Tamaño_muestra_n== .(l)), xlab = "Estimadores")
}
}
En los graficos anteriores se puede apreciar que a menor tamaño de
muestra la distribución registra cola indicando asimetría negativa, a
medida que aumenta el tamaño de muestra los estimadores graficamente se
ajustan con mayor exactitud a una distribución normal, y se concentran
en mayor cantidad alrededor del parámetro (0.90).
Para los diferentes lotes simulados (Proporción plantas enfermas: 0.5, 0.10 y 0.90) se puede observar que a mayor tamaño de muestra los estimadores tienden a su respectivo parámetro poblacional. Los estimadores obtenidos con los lotes simulados de 0,10 y 0,90 proporción de plantas enfermas presentaron distribución asimetrica la cual se puede apreciar en los graficos. Se pudo evidenciar que apesar de que los estimadores de lo lotes simulados con 0,10 y 0,90 proporción de enfermos para los diferentes tamaños de muestra no siguieren una distribución normal de acuerdo a la prueba shapiro wilks, a medida que se incrementaron los tamaños de muestra propuestos, graficamente se converge a la distribución normal (Teorema del límite central).
2. La comparación de tratamientos es una práctica fundamental en las ciencias agropecuarias y para esto a nivel estadístico se cuenta con algunas herramientas para apoyar el proceso de toma de decisiones y lograr concluir con algún grado de confianza que los resultados observados en una muestra son representativos y se pueden asociar a los tratamientos y no se deben únicamente al azar. Por medio una simulación validemos algunos de estos resultados.
a). Suponga un escenario en el cual usted aplicó tratamientos diferentes a dos lotes y desea analizar si alguno de los dos presenta un mejor desempeño en el control de una plaga presente en ambos al momento inicial. Para ello utilizará como criterio de desempeño el tratamiento que menor % de plantas enfermas presente después de un tiempo de aplicación (es decir, si se presentan o no diferencias en las proporciones de enfermos P1 y P2). Realice una simulación en la cual genere dos poblaciones de N1=1000 (Lote1) y N2=1500 (Lote2), además asuma que el porcentaje de individuos (plantas) enfermas en ambos lotes sea la misma 10% (es decir, sin diferencias entre los tratamientos).
#### A. Plantas enfermas=10%
set.seed(0)
LoteN1=c(rep("P_Enfermas",100),rep("P_Sanas",900))
LoteN2=c(rep("P_Enfermas",150),rep("P_Sanas",1350))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.
set.seed(0)
muestra_LoteN1=function(n){
muestraL1=sample(LoteN1,size = n)
return(sum(muestraL1=="P_Enfermas")/n)
}
set.seed(0)
muestra_LoteN2=function(n){
muestraL2=sample(LoteN2,size = n)
return(sum(muestraL2=="P_Enfermas")/n)
}
ML1<-muestra_LoteN1(n = 100)
ML2<-muestra_LoteN2(n = 100)
ML1-ML2## [1] -0.03c). 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?
set.seed(0)
muestra_LoteN1=function(n){
muestraL1=sample(LoteN1,size = n)
return(sum(muestraL1=="P_Enfermas"))
}
set.seed(0)
muestra_LoteN2=function(n){
muestraL2=sample(LoteN2,size = n)
return(sum(muestraL2=="P_Enfermas"))
}
EstimadoresL1=sapply(rep(100,500), muestra_LoteN1)
EstimadoresL2=sapply(rep(100,500), muestra_LoteN2)
Dif<-EstimadoresL1-EstimadoresL2¿Qué tan simétricos son los datos?
par(mfrow=c(1,3))
boxplot(EstimadoresL1, main="Estimadores Lote1",col="lightpink")
boxplot(EstimadoresL2, main="Estimadores Lote2",col="cyan")
boxplot(Dif, main="Diferencia estimadores Lote1-Lote2",col="purple")
En los graficos anteriores se puede apreciar que los estimadores de los Lotes 1 y 2 se concentran alredor de 0.10, mientras que la diferencia de los estimadores se concentra al rededor de 0.
skewness(EstimadoresL1); skewness(EstimadoresL2); skewness(Dif)## [1] 0.2316654## [1] 0.1931634## [1] -0.05401725Los indicadores de asimetría de los estimadores para los Lotes 1 y 2 sugieren una distribución asimetríca positiva, miestras que la diferencia de los estimadores tiene una distribución levemente asimetríca negativa.
¿Son siempre cero las diferencias? Hipotesis nula: P1-P2=0
prop.test(x=c(mean(EstimadoresL1),mean(EstimadoresL2)),n=c(100,100),conf.level = 0.95)##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(mean(EstimadoresL1), mean(EstimadoresL2)) out of c(100, 100)
## X-squared = 1.1951e-30, df = 1, p-value = 1
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## -0.08757818 0.08333818
## sample estimates:
## prop 1 prop 2
## 0.09944 0.10156La proporción de plantas enfermas en el Lote 1 es en promedio del 9,944% y en el lote 2 es en promedio del 10,156%, considerando que el intervalo de confianza para la diferencia de proporciones contiene el valor 0, y que el Valor-P>0,05 se puede considerar que la diferencia entre los estimadores de los lotes puede llegar a ser cero, sin embargo, no se puede afirmar con certeza que siempre esta diferencia sera cero.
d). Realice los puntos b y c para tamaños de muestra n1=n2=5, 10, 15, 20, 30, 50, 60, 100, 200, 500.
### Literal B.
set.seed(0)
muestra_LoteN1=function(n){
muestraL1=sample(LoteN1,size = n)
return(sum(muestraL1=="P_Enfermas")/n)
}
set.seed(0)
muestra_LoteN2=function(n){
muestraL2=sample(LoteN2,size = n)
return(sum(muestraL2=="P_Enfermas")/n)
}
MLS1<-sapply(c(5,10,15,20,30,50,60,100,200,500), muestra_LoteN1)
MLS2<-sapply(c(5,10,15,20,30,50,60,100,200,500), muestra_LoteN2)
data.frame(Tamaño_muestra_n=c(5,10,15,20,30,50,60,100,200,500),Dif_estimadores= MLS1-MLS2)## Tamaño_muestra_n Dif_estimadores
## 1 5 -0.20000000
## 2 10 0.00000000
## 3 15 0.06666667
## 4 20 -0.15000000
## 5 30 0.10000000
## 6 50 -0.04000000
## 7 60 -0.03333333
## 8 100 0.02000000
## 9 200 0.02500000
## 10 500 -0.01600000### Literal C.
set.seed(0)
muestra_LoteN1=function(n){
muestraL1=sample(LoteN1,size = n)
return(sum(muestraL1=="P_Enfermas"))
}
n=c(5,10,15,20,30,50,60,100,200,500)
estimadoresLN1<-data.frame(matrix(NA,nrow=500,ncol=length(n)))
for (i in 1:10) {
estimadoresLN1[,i]=sapply(rep(n[i],500), muestra_LoteN1)
}
set.seed(0)
muestra_LoteN2=function(n){
muestraL2=sample(LoteN2,size = n)
return(sum(muestraL2=="P_Enfermas"))
}
n=c(5,10,15,20,30,50,60,100,200,500)
estimadoresLN2<-data.frame(matrix(NA,nrow=500,ncol=length(n)))
for (i in 1:10) {
estimadoresLN2[,i]=sapply(rep(n[i],500), muestra_LoteN2)
}
Dif12<-estimadoresLN1-estimadoresLN2Compare los resultados de los estimadores (p1-p2) en cuanto a la normalidad
StLote12<-c()
for (i in 1:10) {
StLote12[i]<-shapiro.test(Dif12[,i])$p.value
}
data.frame(Tamaño_muestra_n=c(5,10,15,20,30,50,60,100,200,500),Valor_P=StLote12)## Tamaño_muestra_n Valor_P
## 1 5 3.666629e-16
## 2 10 3.214408e-12
## 3 15 3.459638e-09
## 4 20 9.870193e-09
## 5 30 1.245186e-06
## 6 50 2.956282e-04
## 7 60 4.294763e-04
## 8 100 1.840590e-02
## 9 200 6.336315e-02
## 10 500 1.218896e-01A partir del tamaño de muestra n=200, se puede apreciar que el Valor P de la prueba de normalidad shapiro wilks es mayor al nivel de significancia 0,05 por lo cual, se puede considerar que los estimadores obtenidos en la simulación tienen una distribución normal cuando en n es mayor o igual a 200.
A mayor tamaño de muestra los estimadores convergen al parametro poblacional, por lo cual la diferencia entre los estimadores deberia aproximarse a la diferencia real de los parámetros poblacionales.
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. Plantas enfermas L3=0.10 y L4=0.15
set.seed(0)
LoteN3=c(rep("P_Enfermas",100),rep("P_Sanas",900))
LoteN4=c(rep("P_Enfermas",225),rep("P_Sanas",1275))
#### B.
set.seed(0)
muestra_LoteN3=function(n){
muestraL3=sample(LoteN3,size = n)
return(sum(muestraL3=="P_Enfermas")/n)
}
set.seed(0)
muestra_LoteN4=function(n){
muestraL4=sample(LoteN4,size = n)
return(sum(muestraL4=="P_Enfermas")/n)
}
ML3<-muestra_LoteN3(n = 100)
ML4<-muestra_LoteN4(n = 100)
ML3-ML4## [1] -0.04#### C.
set.seed(0)
muestra_LoteN3=function(n){
muestraL3=sample(LoteN3,size = n)
return(sum(muestraL3=="P_Enfermas"))
}
set.seed(0)
muestra_LoteN4=function(n){
muestraL4=sample(LoteN4,size = n)
return(sum(muestraL4=="P_Enfermas"))
}
EstimadoresL3=sapply(rep(100,500), muestra_LoteN3)
EstimadoresL4=sapply(rep(100,500), muestra_LoteN4)
Dif2<-EstimadoresL3-EstimadoresL4¿Qué tan simétricos son los datos?
par(mfrow=c(1,3))
boxplot(EstimadoresL3, main="Estimadores Lote3",col="lightpink")
boxplot(EstimadoresL4, main="Estimadores Lote4",col="cyan")
boxplot(Dif2, main="Diferencia estimadores Lote3-Lote4",col="purple")
En los graficos anteriores se puede apreciar que los estimadores del Lote 3 se concentran alredor de 0.10 (parámetro), los estimadores del Lote 4 se concentran alredor de 0.15 (parámetro), mientras que la diferencia de los estimadores se concentra al rededor de -5.
skewness(EstimadoresL3); skewness(EstimadoresL4); skewness(Dif2)## [1] 0.2316654## [1] 0.2497016## [1] -0.09049871Los indicadores de asimetría de los estimadores para los Lotes 3 y 4 sugieren una distribución asimetríca positiva, miestras que la diferencia de los estimadores tiene una distribución levemente asimetríca negativa.
¿Son siempre cero las diferencias? Hipotesis nula: P3-P4=0
prop.test(x=c(mean(EstimadoresL3),mean(EstimadoresL4)),n=c(100,100),conf.level = 0.95)##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(mean(EstimadoresL3), mean(EstimadoresL4)) out of c(100, 100)
## X-squared = 0.83768, df = 1, p-value = 0.3601
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## -0.15460096 0.04872096
## sample estimates:
## prop 1 prop 2
## 0.09944 0.15238La proporción de plantas enfermas en el Lote 3 es en promedio del 9,944% y en el lote 4 es en promedio del 15,238%, considerando que el intervalo de confianza para la diferencia de proporciones contiene el valor 0, y que el Valor-P>0,05 se puede considerar que la diferencia entre los estimadores de los lotes puede llegar a ser cero, sin embargo, no se puede afirmar con certeza que siempre esta diferencia sera cero.
#### D.
### Literal B.
set.seed(0)
muestra_LoteN3=function(n){
muestraL3=sample(LoteN3,size = n)
return(sum(muestraL3=="P_Enfermas")/n)
}
set.seed(0)
muestra_LoteN4=function(n){
muestraL4=sample(LoteN4,size = n)
return(sum(muestraL4=="P_Enfermas")/n)
}
MLS3<-sapply(c(5,10,15,20,30,50,60,100,200,500), muestra_LoteN3)
MLS4<-sapply(c(5,10,15,20,30,50,60,100,200,500), muestra_LoteN4)
data.frame(Tamaño_muestra_n=c(5,10,15,20,30,50,60,100,200,500),Dif_estimadores= MLS3-MLS4)## Tamaño_muestra_n Dif_estimadores
## 1 5 -0.20000000
## 2 10 0.00000000
## 3 15 0.00000000
## 4 20 -0.20000000
## 5 30 0.03333333
## 6 50 -0.04000000
## 7 60 -0.10000000
## 8 100 -0.04000000
## 9 200 -0.03500000
## 10 500 -0.06800000#### D.
### Literal C.
set.seed(0)
muestra_LoteN3=function(n){
muestraL3=sample(LoteN3,size = n)
return(sum(muestraL3=="P_Enfermas"))
}
n=c(5,10,15,20,30,50,60,100,200,500)
estimadoresLN3<-data.frame(matrix(NA,nrow=500,ncol=length(n)))
for (i in 1:10) {
estimadoresLN3[,i]=sapply(rep(n[i],500), muestra_LoteN3)
}
set.seed(0)
muestra_LoteN4=function(n){
muestraL4=sample(LoteN4,size = n)
return(sum(muestraL4=="P_Enfermas"))
}
n=c(5,10,15,20,30,50,60,100,200,500)
estimadoresLN4<-data.frame(matrix(NA,nrow=500,ncol=length(n)))
for (i in 1:10) {
estimadoresLN4[,i]=sapply(rep(n[i],500), muestra_LoteN4)
}
Dif34<-estimadoresLN3-estimadoresLN4Compare los resultados de los estimadores (p1-p2) en cuanto a la normalidad
StLote34<-c()
for (i in 1:10) {
StLote34[i]<-shapiro.test(Dif34[,i])$p.value
}
data.frame(Tamaño_muestra_n=c(5,10,15,20,30,50,60,100,200,500),Valor_P=StLote34)## Tamaño_muestra_n Valor_P
## 1 5 5.268020e-15
## 2 10 8.582363e-11
## 3 15 3.600233e-08
## 4 20 3.589930e-07
## 5 30 2.505092e-05
## 6 50 6.078612e-04
## 7 60 1.157319e-03
## 8 100 9.170776e-03
## 9 200 6.656074e-02
## 10 500 2.744798e-01A partir del tamaño de muestra n=200, se puede apreciar que el Valor P de la prueba de normalidad shapiro wilks es mayor al nivel de significancia 0,05 por lo cual, se puede considerar que los estimadores obtenidos en la simulación tienen una distribución normal cuando en n es mayor o igual a 200.
El actual ejercicio respecto al ejercicio con los dos lotes de propocion de plantas enfermas de 0.10 radica en que al evaluar las diferencias de los estimadores se puede eviderciar que para el primer ejecicio a mayor tamaño de muestra la diferencia entre los estimadores converge a 0, mientras que en el ultimo ejercicio con lotes de proporciones de plantas enfermas diferentes a mayor tamaño de muestra esa diferencia converge a -5.
3. Resumen artículos “Statistical Errors: P values, the gold standard of statistical validity, are not as reliable as many scientists assume” y “Statisticians issue warning on P values: Statement aims to halt missteps in the quest for certainty”
En la década de 1920 el estadístico británico Ronald Fisher introdujo el Valor P como una manera informal para juzgar si la evidencia era significativa; primero se debía establecer una “hipótesis nula” que se quería refutar, se realizaba el experimento y se evaluaba si los resultados eran consistentes o se podrían deber al azar; partiendo de la suposición que la hipótesis fuera cierta, se calculaba la posibilidad de obtener resultados tan extremos como los que se observaban (Valor P), de acuerdo a Fisher cuanto más pequeño era el Valor P mayor era la probabilidad de que la hipótesis nula fuera falsa. Fisher pretendía que el Valor P fuera un método informal, sin embargo, pronto se convirtió en un movimiento para la toma de decisiones basado en evidencia lo más rigurosa y objetiva posible, fue entonces cuando el Valor P de 0,05 se consagro como “Estadísticamente significativo”.
Desde su existencia el Valor P ha sido fuertemente cuestionado por no ser tan fiable ni tan objetivo como suponen la mayoría de científicos; una de las mayores controversias gira en torno a la reproducibilidad de los estudios, en el año 2005 el epidemiólogo John Loannidis sugirió que la mayoría de los hallazgos publicados eran falsos, obligando así a los científicos a repensar como evaluar los resultados con métodos diferentes al Valor P debido a una serie de problemas de replicación de alto perfil, posteriormente la Asociación Estadounidense de Estadística (ASA) emitió recomendaciones explicitas preocupados por el uso incorrecto de este, sugiriendo no sacar conclusiones basadas únicamente en Valor p; los críticos también resaltan que el Valor P puede fomentar un pensamiento confuso, pues se podría llegar a pensar que su valor representa la probabilidad de que un resultado fuera una falsa alarma, cuando realmente solo se pueden concluir los datos entorno a una hipótesis nula específica, cabe mencionar que también se debe tener en cuenta el tamaño real de un efecto y evaluar si este realmente existe. Actualmente el Valor P se ha vuelto un método tan riguroso y concluyente que quizás se podría llegar a manipular y monitorear los datos hasta llegar al resultado deseado, estas prácticas tienen el efecto de convertir los descubrimientos exploratorios en afirmaciones sólidas, las cuales se desvaneces tras la replicación.
Un Valor P de 0,05 no significa que haya un 95% de posibilidad de que una hipótesis dada sea correcta, en cambio, significa que la hipótesis nula es verdadera y todas las demás suposiciones son válidas, y que hay un 5% de posibilidades de obtener un resultado al menos tan extremo como el observado.
A pesar de las diversas críticas por cambiar el concepto y restar importancia al Valor P, es un proceso largo y demorado, pues el marco básico de la estadística no ha cambiado desde que fue introducido por Fisher, Neyman y Pearson. Entre las soluciones pensadas para esta problemática, Cumming consideran importante que los investigadores informen el tamaño de los efectos y los intervalos de confianza, otros estadísticos abogan por reemplazar el Valor P con métodos basados en la regla de Bayes, otros abogan a probar diferentes metodología para probar los datos, y Simonsohn argumenta que es necesario que los investigadores especifiquen las inclusiones y exclusiones de las muestras, sus tamaños y cualquier modificación que haya sido relevante para llegar a los resultados deseados.
El uso del Valor P como criterio de decisión en inferencia estadística, debe contemplar las diferentes problemáticas expuestas anteriormente, es fundamental evaluar la selección de las muestras y ser explícitos en el marco muestral de los datos, se deben tener en cuenta el contexto de la información y las diversas metodologías a parte del Valor P para evaluar la significancia de los resultados, y no se deben hacer conclusiones a priori de la población a partir de los resultados del análisis exploratorio muestral.