Taller 2: Inferencia Estadística y Simulación
John Jairo Bedoya Saenz y Liz Mary Gutierrez Rendón
2022-09-17
PUNTO 1
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 es del 50%.
lote_pl = c(rep("Plantas_Buenas",500), rep("Plantas_Enfermas",500))
#lote_pl
table(lote_pl)/1000## lote_pl
## Plantas_Buenas Plantas_Enfermas
## 0.5 0.5head(lote_pl,10)## [1] "Plantas_Buenas" "Plantas_Buenas" "Plantas_Buenas" "Plantas_Buenas"
## [5] "Plantas_Buenas" "Plantas_Buenas" "Plantas_Buenas" "Plantas_Buenas"
## [9] "Plantas_Buenas" "Plantas_Buenas"tail(lote_pl,10)## [1] "Plantas_Enfermas" "Plantas_Enfermas" "Plantas_Enfermas" "Plantas_Enfermas"
## [5] "Plantas_Enfermas" "Plantas_Enfermas" "Plantas_Enfermas" "Plantas_Enfermas"
## [9] "Plantas_Enfermas" "Plantas_Enfermas"Porcentaje real de plantas enfermas = P = 50% = 500/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_enfermas=function(n){
muestrap=sample(lote_pl,size = n)
return(sum(muestrap=="Plantas_Enfermas")/n)
}
calc_enfermas(n=100)## [1] 0.49c. 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 que pasa en cuanto a variabilidad?.
estimadores_p= sapply(rep(100,500), calc_enfermas)
estimadores_p## [1] 0.53 0.55 0.56 0.48 0.46 0.54 0.50 0.54 0.50 0.52 0.56 0.59 0.47 0.52 0.57
## [16] 0.50 0.53 0.48 0.53 0.47 0.52 0.42 0.50 0.56 0.52 0.44 0.43 0.50 0.55 0.51
## [31] 0.57 0.51 0.49 0.55 0.51 0.41 0.53 0.55 0.53 0.63 0.40 0.61 0.48 0.56 0.50
## [46] 0.43 0.47 0.47 0.42 0.44 0.51 0.41 0.52 0.53 0.52 0.53 0.41 0.47 0.50 0.43
## [61] 0.47 0.48 0.52 0.52 0.52 0.44 0.59 0.45 0.48 0.50 0.47 0.48 0.50 0.48 0.44
## [76] 0.49 0.57 0.54 0.57 0.56 0.55 0.51 0.50 0.55 0.56 0.52 0.46 0.54 0.44 0.45
## [91] 0.54 0.50 0.54 0.51 0.55 0.47 0.59 0.50 0.54 0.49 0.45 0.54 0.54 0.47 0.53
## [106] 0.44 0.48 0.40 0.48 0.48 0.45 0.58 0.57 0.48 0.45 0.56 0.50 0.56 0.47 0.54
## [121] 0.50 0.49 0.44 0.55 0.46 0.56 0.43 0.49 0.65 0.49 0.50 0.50 0.54 0.41 0.49
## [136] 0.52 0.61 0.56 0.49 0.42 0.50 0.51 0.41 0.48 0.49 0.48 0.54 0.50 0.52 0.45
## [151] 0.50 0.43 0.60 0.51 0.41 0.57 0.52 0.46 0.58 0.45 0.45 0.41 0.44 0.54 0.51
## [166] 0.47 0.53 0.55 0.45 0.47 0.60 0.44 0.48 0.50 0.51 0.47 0.52 0.45 0.51 0.62
## [181] 0.52 0.41 0.56 0.52 0.52 0.51 0.50 0.45 0.53 0.46 0.42 0.56 0.48 0.50 0.46
## [196] 0.46 0.55 0.51 0.53 0.49 0.50 0.52 0.39 0.44 0.49 0.51 0.55 0.45 0.51 0.56
## [211] 0.48 0.44 0.54 0.41 0.48 0.53 0.54 0.49 0.47 0.50 0.45 0.55 0.56 0.49 0.54
## [226] 0.51 0.50 0.49 0.47 0.54 0.53 0.51 0.51 0.51 0.45 0.46 0.54 0.53 0.56 0.45
## [241] 0.52 0.51 0.42 0.48 0.46 0.48 0.65 0.49 0.42 0.43 0.39 0.50 0.57 0.53 0.49
## [256] 0.47 0.53 0.46 0.51 0.46 0.45 0.44 0.43 0.54 0.49 0.51 0.55 0.53 0.49 0.50
## [271] 0.48 0.46 0.50 0.47 0.56 0.42 0.54 0.39 0.51 0.48 0.52 0.52 0.55 0.48 0.50
## [286] 0.57 0.50 0.52 0.48 0.43 0.50 0.46 0.49 0.43 0.46 0.50 0.52 0.51 0.51 0.47
## [301] 0.48 0.51 0.49 0.49 0.56 0.47 0.46 0.47 0.42 0.51 0.49 0.50 0.48 0.50 0.56
## [316] 0.50 0.42 0.49 0.47 0.42 0.55 0.47 0.53 0.51 0.52 0.48 0.58 0.48 0.49 0.49
## [331] 0.56 0.48 0.51 0.52 0.54 0.45 0.49 0.50 0.53 0.49 0.50 0.52 0.55 0.44 0.57
## [346] 0.49 0.55 0.47 0.55 0.52 0.55 0.54 0.48 0.52 0.53 0.58 0.56 0.47 0.55 0.52
## [361] 0.44 0.48 0.50 0.58 0.52 0.57 0.41 0.47 0.47 0.46 0.51 0.54 0.41 0.49 0.45
## [376] 0.43 0.53 0.53 0.44 0.41 0.60 0.51 0.44 0.58 0.55 0.64 0.45 0.52 0.49 0.46
## [391] 0.48 0.53 0.54 0.56 0.47 0.54 0.43 0.56 0.53 0.58 0.50 0.51 0.46 0.54 0.45
## [406] 0.52 0.43 0.41 0.44 0.54 0.47 0.45 0.54 0.46 0.49 0.47 0.46 0.52 0.44 0.50
## [421] 0.40 0.53 0.44 0.45 0.50 0.60 0.47 0.48 0.60 0.46 0.46 0.52 0.45 0.47 0.47
## [436] 0.53 0.50 0.50 0.46 0.55 0.53 0.47 0.58 0.47 0.49 0.51 0.53 0.49 0.47 0.47
## [451] 0.45 0.55 0.51 0.49 0.53 0.49 0.51 0.48 0.51 0.60 0.55 0.50 0.38 0.51 0.47
## [466] 0.52 0.52 0.50 0.53 0.51 0.44 0.48 0.54 0.43 0.49 0.51 0.46 0.53 0.49 0.44
## [481] 0.51 0.54 0.49 0.54 0.47 0.49 0.50 0.49 0.53 0.42 0.36 0.42 0.52 0.57 0.50
## [496] 0.44 0.54 0.56 0.49 0.44hist(estimadores_p, col="darkblue", main = "Histograma Estimadores (500 veces)")
summary(estimadores_p)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.3600 0.4700 0.5000 0.4989 0.5300 0.6500sd(estimadores_p)## [1] 0.04758065Simetría: El histograma de las estimaciones realizadas parece ser simétrico dado que que los valores se concentran al rededor de 0.5 y comparando la media con un 0.5 versus la mediana con un 0.49, se observa que estan muy cerca una de la otra, lo que se infiere que tiene una distribución simétrica.
Sesgo: Se observa que la concentración de los valores se encuentran alrededor del parámetro real 0.5, por tanto hay poco sesgo, lo que quiere decir que las muestras fueron tomadas adecuadamente.
Variablidad: El rango (max-min), muestra que se presenta una dispersión de los datos en una longitud angosta que es donde se hallan todos los datos muestrales. Con respecto a la desviación estandar, se observa que la desviación de plantas malas que se presentan en los lotes con respecto al promedio es muy bajo.
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).
estimadores_p5=sapply(rep(5,500), calc_enfermas)
estimadores_p10=sapply(rep(10,500), calc_enfermas)
estimadores_p15=sapply(rep(15,500), calc_enfermas)
estimadores_p20=sapply(rep(20,500), calc_enfermas)
estimadores_p30=sapply(rep(30,500), calc_enfermas)
estimadores_p50=sapply(rep(50,500), calc_enfermas)
estimadores_p60=sapply(rep(60,500), calc_enfermas)
estimadores_p100=sapply(rep(100,500), calc_enfermas)
estimadores_p200=sapply(rep(200,500), calc_enfermas)
estimadores_p500=sapply(rep(500,500), calc_enfermas)
res_p=data.frame(estimadores_p5,estimadores_p10,estimadores_p15,estimadores_p20,estimadores_p30,estimadores_p50,estimadores_p60,estimadores_p100,estimadores_p200,estimadores_p500)
#summary(res_p)
# Ajuste para la tabla:
indicadores <- cbind(apply(res_p, 2, mean),
apply(res_p, 2, median),
apply(res_p, 2, sd),
apply(res_p, 2, min),
apply(res_p, 2, max))
colnames(indicadores) <- c("media","mediana","sd","min","max")
round(indicadores,3)## media mediana sd min max
## estimadores_p5 0.506 0.600 0.220 0.000 1.000
## estimadores_p10 0.501 0.500 0.157 0.000 0.900
## estimadores_p15 0.496 0.467 0.125 0.133 0.867
## estimadores_p20 0.491 0.500 0.112 0.200 0.800
## estimadores_p30 0.499 0.500 0.089 0.233 0.767
## estimadores_p50 0.497 0.500 0.071 0.280 0.720
## estimadores_p60 0.497 0.500 0.066 0.317 0.700
## estimadores_p100 0.504 0.500 0.047 0.380 0.650
## estimadores_p200 0.500 0.500 0.031 0.405 0.590
## estimadores_p500 0.501 0.500 0.016 0.458 0.546boxplot(res_p)
abline(h=0.5,col="red",lwd=1)
En los calculos de media, mediana, sd, maximo, minimo y el diagrama de cajas, se observa que entre mayor es el tamaño de la muestra, menor es la dispersión de los datos y simetría de los mismos, que la media y mediana tienden a ser iguales al mismo tiempo que disminuye el valor de la desviavión estandar al igual que la longitud del rango, se evidencia tambien que es mas precisa la información con las muesras mas grandes, para el caso de n=500 observamos que el 50% de los datos estan mas cercanos al 50% e inclusive los bigotes del grafico se encuentran cercanos al estimador del 50%. Por otro lado, para las muestras de tamaño 5 hasta el tamaño de muestras de 60 datos, se presenta muchas dispersión de los datos, asimetría de los mismos, datos atípicos y sesgos alejados del parametro con una gran desviación estandar.
Pruebas de bondad y ajuste (shapiro wilks)
Hipótesis nula: las estimaciones provienen de una distribución Normal
Hipótesis alternativa: Las estimaciones NO provienen de una distribución Normal
Criterio de rechazo: Dado una significancia (Asumiendo alpha = 0.05)
NO Se rechaza H_0 si el p-valor de la prueba es > a 0.05
Se rechaza H_0 si el p-valor de la prueba es <= a 0.05
shapiro.test(estimadores_p5) ##
## Shapiro-Wilk normality test
##
## data: estimadores_p5
## W = 0.92649, p-value = 6.221e-15shapiro.test(estimadores_p10) ##
## Shapiro-Wilk normality test
##
## data: estimadores_p10
## W = 0.96339, p-value = 7.93e-10shapiro.test(estimadores_p15) ##
## Shapiro-Wilk normality test
##
## data: estimadores_p15
## W = 0.97346, p-value = 7.091e-08shapiro.test(estimadores_p20) ##
## Shapiro-Wilk normality test
##
## data: estimadores_p20
## W = 0.97973, p-value = 1.968e-06shapiro.test(estimadores_p30) ##
## Shapiro-Wilk normality test
##
## data: estimadores_p30
## W = 0.98652, p-value = 0.0001394shapiro.test(estimadores_p50) ##
## Shapiro-Wilk normality test
##
## data: estimadores_p50
## W = 0.99172, p-value = 0.006892shapiro.test(estimadores_p60) ##
## Shapiro-Wilk normality test
##
## data: estimadores_p60
## W = 0.99226, p-value = 0.01072shapiro.test(estimadores_p100)##
## Shapiro-Wilk normality test
##
## data: estimadores_p100
## W = 0.99242, p-value = 0.01223shapiro.test(estimadores_p200)##
## Shapiro-Wilk normality test
##
## data: estimadores_p200
## W = 0.99473, p-value = 0.08455shapiro.test(estimadores_p500)##
## Shapiro-Wilk normality test
##
## data: estimadores_p500
## W = 0.99546, p-value = 0.1549De acuerdo con las pruebas de bondad y ajuste (shapiro wilks), tiende a ser rechazada la hipotesis nula en los estimadores de muestra de 5 a 60 dado que el p-valor se encuentrapor debajo de la significancia asumida del 0.05, evidenciandose que dichas estimaciones NO provienen de una distribución Normal. Por el contrario, la hipotesis nula tiende a no ser rechazada para los estimadores de muestras 100, 200 y 500 con un p-valor por encima del alpha = 0.05, evidenciandose que dichas estimaciones provienen de una distribución Normal.
Gráficos Cuantil-Cuantil (Q-Q plots) - grafico qq de normalidad
par(mfrow=c(2,2))
qqnorm(estimadores_p5, main = "QQ-plot para n=5")
qqline(estimadores_p5, col="red")
qqnorm(estimadores_p10, main = "QQ-plot para n=10")
qqline(estimadores_p10, col="red")
qqnorm(estimadores_p15, main = "QQ-plot para n=15")
qqline(estimadores_p15, col="red")
qqnorm(estimadores_p20, main = "QQ-plot para n=20")
qqline(estimadores_p20, col="red")
qqnorm(estimadores_p30, main = "QQ-plot para n=30")
qqline(estimadores_p30, col="red")
qqnorm(estimadores_p50, main = "QQ-plot para n=50")
qqline(estimadores_p50, col="red")
qqnorm(estimadores_p60, main = "QQ-plot para n=60")
qqline(estimadores_p60, col="red")
qqnorm(estimadores_p100, main = "QQ-plot para n=100")
qqline(estimadores_p100, col="red")
qqnorm(estimadores_p200, main = "QQ-plot para n=200")
qqline(estimadores_p200, col="red")
qqnorm(estimadores_p500, main = "QQ-plot para n=500")
qqline(estimadores_p500, col="red")
Al momento de graficar un vector de estimaciones vs los valores teóricos de una distribución normal, se observa entonces que a medida que las estimaciones van incrementando el tamaño de muestra, el vector se distribuye Normal y el gráfico qq se asemeja cada vez más a la de la forma X=Y
e. Repita toda la simulación (puntos a – d) pero ahora con lotes con 10% y 90% de plantas enfermas. Concluya todo el ejercicio.
Plantas enfermas = 10%
#a.
lote_pl10 = c(rep("Plantas_Buenas",900), rep("Plantas_Enfermas",100))
#lote_pl10
table(lote_pl10)/1000## lote_pl10
## Plantas_Buenas Plantas_Enfermas
## 0.9 0.1#b.
calc_enfermas10=function(n){
muestrap10=sample(lote_pl10,size = n)
return(sum(muestrap10=="Plantas_Enfermas")/n)
}
calc_enfermas10(n=100)## [1] 0.11#c.
estimadores_p_p10= sapply(rep(100,500), calc_enfermas10)
estimadores_p_p10## [1] 0.13 0.11 0.07 0.10 0.08 0.16 0.13 0.09 0.14 0.11 0.06 0.11 0.10 0.13 0.10
## [16] 0.09 0.13 0.09 0.10 0.12 0.13 0.09 0.10 0.02 0.07 0.12 0.11 0.11 0.09 0.12
## [31] 0.06 0.08 0.06 0.11 0.10 0.07 0.10 0.09 0.08 0.03 0.10 0.08 0.06 0.12 0.09
## [46] 0.09 0.10 0.07 0.10 0.05 0.13 0.12 0.08 0.14 0.13 0.09 0.18 0.11 0.08 0.06
## [61] 0.08 0.10 0.11 0.13 0.08 0.11 0.10 0.10 0.11 0.10 0.07 0.09 0.07 0.08 0.12
## [76] 0.10 0.13 0.08 0.09 0.08 0.08 0.09 0.10 0.03 0.08 0.08 0.08 0.03 0.11 0.10
## [91] 0.11 0.07 0.12 0.08 0.11 0.10 0.09 0.13 0.17 0.13 0.06 0.12 0.13 0.10 0.11
## [106] 0.09 0.10 0.09 0.10 0.10 0.08 0.10 0.10 0.11 0.12 0.09 0.05 0.11 0.05 0.13
## [121] 0.06 0.17 0.12 0.14 0.09 0.07 0.10 0.16 0.15 0.12 0.07 0.10 0.08 0.05 0.12
## [136] 0.14 0.15 0.10 0.12 0.09 0.11 0.09 0.06 0.20 0.09 0.08 0.09 0.08 0.13 0.19
## [151] 0.14 0.08 0.13 0.15 0.10 0.09 0.13 0.10 0.09 0.06 0.10 0.09 0.12 0.11 0.07
## [166] 0.06 0.11 0.09 0.16 0.13 0.11 0.13 0.05 0.04 0.09 0.10 0.11 0.11 0.05 0.11
## [181] 0.08 0.06 0.10 0.17 0.08 0.07 0.13 0.09 0.11 0.12 0.14 0.09 0.12 0.05 0.12
## [196] 0.07 0.11 0.10 0.09 0.13 0.08 0.11 0.07 0.08 0.11 0.09 0.11 0.11 0.11 0.13
## [211] 0.07 0.09 0.08 0.04 0.09 0.08 0.12 0.13 0.09 0.14 0.07 0.11 0.17 0.16 0.14
## [226] 0.07 0.05 0.12 0.08 0.16 0.16 0.06 0.12 0.06 0.11 0.08 0.12 0.08 0.09 0.09
## [241] 0.08 0.13 0.08 0.11 0.09 0.09 0.06 0.07 0.11 0.10 0.09 0.11 0.14 0.05 0.06
## [256] 0.12 0.14 0.13 0.10 0.06 0.06 0.08 0.09 0.06 0.09 0.10 0.09 0.13 0.14 0.11
## [271] 0.13 0.08 0.11 0.16 0.10 0.05 0.06 0.08 0.06 0.12 0.19 0.09 0.12 0.11 0.07
## [286] 0.12 0.08 0.11 0.09 0.13 0.08 0.05 0.05 0.12 0.09 0.07 0.09 0.09 0.10 0.09
## [301] 0.12 0.09 0.09 0.08 0.12 0.12 0.08 0.12 0.07 0.12 0.09 0.11 0.09 0.15 0.13
## [316] 0.08 0.05 0.08 0.11 0.08 0.20 0.07 0.07 0.04 0.10 0.09 0.07 0.10 0.11 0.08
## [331] 0.08 0.09 0.10 0.13 0.07 0.09 0.10 0.07 0.12 0.10 0.06 0.08 0.13 0.09 0.10
## [346] 0.11 0.08 0.05 0.14 0.09 0.10 0.10 0.08 0.17 0.09 0.13 0.17 0.09 0.12 0.12
## [361] 0.09 0.10 0.07 0.12 0.12 0.05 0.10 0.06 0.08 0.10 0.14 0.10 0.10 0.08 0.09
## [376] 0.05 0.04 0.11 0.08 0.15 0.13 0.13 0.12 0.12 0.10 0.12 0.14 0.02 0.09 0.17
## [391] 0.11 0.12 0.12 0.07 0.09 0.10 0.04 0.09 0.06 0.14 0.10 0.07 0.14 0.06 0.07
## [406] 0.06 0.14 0.06 0.13 0.07 0.11 0.12 0.12 0.09 0.11 0.05 0.12 0.07 0.03 0.10
## [421] 0.09 0.10 0.11 0.13 0.11 0.12 0.06 0.06 0.12 0.13 0.12 0.06 0.12 0.08 0.10
## [436] 0.15 0.12 0.10 0.08 0.07 0.07 0.11 0.07 0.09 0.16 0.07 0.13 0.09 0.11 0.11
## [451] 0.06 0.07 0.09 0.11 0.09 0.16 0.11 0.15 0.14 0.11 0.10 0.15 0.09 0.07 0.08
## [466] 0.12 0.09 0.07 0.07 0.09 0.08 0.11 0.15 0.14 0.11 0.07 0.06 0.14 0.13 0.06
## [481] 0.08 0.05 0.08 0.08 0.10 0.11 0.10 0.05 0.12 0.14 0.09 0.09 0.09 0.16 0.09
## [496] 0.08 0.11 0.12 0.10 0.13hist(estimadores_p_p10, main = "Histograma Estimadores (500 veces)")
summary(estimadores_p_p10)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0200 0.0800 0.1000 0.0991 0.1200 0.2000sd(estimadores_p_p10)## [1] 0.02990617#d.
estimadores_p5_10=sapply(rep(5,500), calc_enfermas10)
estimadores_p10_10=sapply(rep(10,500), calc_enfermas10)
estimadores_p15_10=sapply(rep(15,500), calc_enfermas10)
estimadores_p20_10=sapply(rep(20,500), calc_enfermas10)
estimadores_p30_10=sapply(rep(30,500), calc_enfermas10)
estimadores_p50_10=sapply(rep(50,500), calc_enfermas10)
estimadores_p60_10=sapply(rep(60,500), calc_enfermas10)
estimadores_p100_10=sapply(rep(100,500), calc_enfermas10)
estimadores_p200_10=sapply(rep(200,500), calc_enfermas10)
estimadores_p500_10=sapply(rep(500,500), calc_enfermas10)
res_p_10=data.frame(estimadores_p5_10,estimadores_p10_10,estimadores_p15_10,estimadores_p20_10,estimadores_p30_10,estimadores_p50_10,estimadores_p60_10,estimadores_p100_10,estimadores_p200_10,estimadores_p500_10)
# Ajuste para tabla:
indicadores10 <- cbind(apply(res_p_10, 2, mean),
apply(res_p_10, 2, median),
apply(res_p_10, 2, sd),
apply(res_p_10, 2, min),
apply(res_p_10, 2, max))
colnames(indicadores10) <- c("media","mediana","sd","min","max")
round(indicadores10,3)## media mediana sd min max
## estimadores_p5_10 0.095 0.000 0.131 0.000 0.600
## estimadores_p10_10 0.098 0.100 0.095 0.000 0.500
## estimadores_p15_10 0.096 0.067 0.077 0.000 0.400
## estimadores_p20_10 0.099 0.100 0.069 0.000 0.350
## estimadores_p30_10 0.100 0.100 0.055 0.000 0.333
## estimadores_p50_10 0.101 0.100 0.041 0.000 0.240
## estimadores_p60_10 0.099 0.100 0.035 0.000 0.217
## estimadores_p100_10 0.100 0.100 0.027 0.030 0.180
## estimadores_p200_10 0.099 0.100 0.019 0.045 0.155
## estimadores_p500_10 0.100 0.100 0.010 0.068 0.138boxplot(res_p_10)
abline(h=0.1,col="red",lwd=1)
# Pruebas de bondad y ajuste (shapiro wilks)
# Hipótesis nula: las estimaciones provienen de una distribución Normal
# Hipótesis alternativa: Las estimaciones NO provienen de una distribución Normal
# Criterio de rechazo: Dado una significancia (Asumiendo alpha = 0.05)
# NO Se rechaza H_0 si el p-valor de la prueba es > a 0.05
# Se rechaza H_0 si el p-valor de la prueba es <= a 0.05
shapiro.test(estimadores_p5_10) ##
## Shapiro-Wilk normality test
##
## data: estimadores_p5_10
## W = 0.69638, p-value < 2.2e-16shapiro.test(estimadores_p10_10) ##
## Shapiro-Wilk normality test
##
## data: estimadores_p10_10
## W = 0.83198, p-value < 2.2e-16shapiro.test(estimadores_p15_10) ##
## Shapiro-Wilk normality test
##
## data: estimadores_p15_10
## W = 0.89143, p-value < 2.2e-16shapiro.test(estimadores_p20_10) ##
## Shapiro-Wilk normality test
##
## data: estimadores_p20_10
## W = 0.91591, p-value = 4.689e-16shapiro.test(estimadores_p30_10) ##
## Shapiro-Wilk normality test
##
## data: estimadores_p30_10
## W = 0.94149, p-value = 3.954e-13shapiro.test(estimadores_p50_10) ##
## Shapiro-Wilk normality test
##
## data: estimadores_p50_10
## W = 0.97489, p-value = 1.451e-07shapiro.test(estimadores_p60_10) ##
## Shapiro-Wilk normality test
##
## data: estimadores_p60_10
## W = 0.97943, p-value = 1.659e-06shapiro.test(estimadores_p100_10)##
## Shapiro-Wilk normality test
##
## data: estimadores_p100_10
## W = 0.98554, p-value = 7.147e-05shapiro.test(estimadores_p200_10)##
## Shapiro-Wilk normality test
##
## data: estimadores_p200_10
## W = 0.98791, p-value = 0.0003739shapiro.test(estimadores_p500_10)##
## Shapiro-Wilk normality test
##
## data: estimadores_p500_10
## W = 0.99239, p-value = 0.01192par(mfrow=c(2,2))
qqnorm(estimadores_p5_10, main = "QQ-plot para n=5")
qqline(estimadores_p5_10, col="red")
qqnorm(estimadores_p10_10, main = "QQ-plot para n=10")
qqline(estimadores_p10_10, col="red")
qqnorm(estimadores_p15_10, main = "QQ-plot para n=15")
qqline(estimadores_p15_10, col="red")
qqnorm(estimadores_p20_10, main = "QQ-plot para n=20")
qqline(estimadores_p20_10, col="red")
qqnorm(estimadores_p30_10, main = "QQ-plot para n=30")
qqline(estimadores_p30_10, col="red")
qqnorm(estimadores_p50_10, main = "QQ-plot para n=50")
qqline(estimadores_p50_10, col="red")
qqnorm(estimadores_p60_10, main = "QQ-plot para n=60")
qqline(estimadores_p60_10, col="red")
qqnorm(estimadores_p100_10, main = "QQ-plot para n=100")
qqline(estimadores_p100_10, col="red")
qqnorm(estimadores_p200_10, main = "QQ-plot para n=200")
qqline(estimadores_p200_10, col="red")
qqnorm(estimadores_p500_10, main = "QQ-plot para n=500")
qqline(estimadores_p500_10, col="red")
Plantas enfermas = 90%
#a.
lote_pl90 = c(rep("Plantas_Buenas",100), rep("Plantas_Enfermas",900))
#lote_pl90
table(lote_pl90)/1000## lote_pl90
## Plantas_Buenas Plantas_Enfermas
## 0.1 0.9#b.
calc_enfermas90=function(n){
muestrap90=sample(lote_pl90,size = n)
return(sum(muestrap90=="Plantas_Enfermas")/n)
}
calc_enfermas90(n=100)## [1] 0.91#c.
estimadores_p_p90= sapply(rep(100,500), calc_enfermas90)
estimadores_p_p90## [1] 0.92 0.86 0.96 0.88 0.86 0.93 0.88 0.92 0.88 0.97 0.92 0.91 0.89 0.91 0.98
## [16] 0.86 0.87 0.94 0.91 0.90 0.89 0.87 0.90 0.92 0.87 0.92 0.91 0.88 0.87 0.87
## [31] 0.88 0.92 0.91 0.86 0.81 0.93 0.87 0.94 0.90 0.91 0.97 0.90 0.88 0.87 0.84
## [46] 0.91 0.93 0.94 0.96 0.95 0.90 0.87 0.87 0.93 0.85 0.95 0.94 0.90 0.88 0.91
## [61] 0.92 0.87 0.94 0.92 0.91 0.92 0.86 0.89 0.94 0.89 0.89 0.91 0.89 0.93 0.86
## [76] 0.93 0.92 0.88 0.91 0.90 0.89 0.93 0.91 0.91 0.92 0.90 0.88 0.82 0.90 0.90
## [91] 0.91 0.89 0.91 0.86 0.89 0.88 0.93 0.93 0.91 0.92 0.91 0.88 0.83 0.96 0.92
## [106] 0.95 0.86 0.92 0.93 0.86 0.88 0.93 0.91 0.87 0.91 0.90 0.90 0.92 0.85 0.95
## [121] 0.93 0.87 0.91 0.93 0.84 0.88 0.90 0.89 0.90 0.92 0.90 0.90 0.89 0.92 0.89
## [136] 0.91 0.88 0.90 0.92 0.90 0.94 0.92 0.87 0.84 0.92 0.90 0.89 0.93 0.96 0.89
## [151] 0.90 0.89 0.93 0.91 0.90 0.91 0.93 0.96 0.92 0.92 0.96 0.91 0.89 0.91 0.95
## [166] 0.90 0.86 0.91 0.91 0.91 0.86 0.87 0.87 0.92 0.91 0.89 0.89 0.93 0.91 0.96
## [181] 0.83 0.88 0.89 0.89 0.91 0.90 0.88 0.89 0.90 0.90 0.90 0.84 0.87 0.92 0.84
## [196] 0.89 0.87 0.87 0.92 0.92 0.92 0.90 0.88 0.90 0.96 0.88 0.89 0.87 0.96 0.88
## [211] 0.94 0.88 0.87 0.92 0.94 0.88 0.94 0.95 0.87 0.86 0.89 0.89 0.92 0.94 0.90
## [226] 0.88 0.89 0.90 0.88 0.93 0.92 0.91 0.91 0.91 0.92 0.85 0.88 0.85 0.89 0.88
## [241] 0.91 0.94 0.91 0.88 0.88 0.91 0.93 0.82 0.86 0.84 0.93 0.88 0.92 0.89 0.93
## [256] 0.91 0.85 0.90 0.87 0.91 0.93 0.89 0.87 0.92 0.92 0.93 0.88 0.88 0.90 0.94
## [271] 0.88 0.87 0.89 0.91 0.88 0.95 0.85 0.89 0.89 0.91 0.91 0.86 0.89 0.94 0.89
## [286] 0.91 0.91 0.90 0.91 0.90 0.87 0.88 0.92 0.91 0.89 0.86 0.90 0.92 0.94 0.90
## [301] 0.90 0.93 0.88 0.89 0.91 0.88 0.89 0.91 0.92 0.91 0.93 0.90 0.87 0.94 0.95
## [316] 0.92 0.87 0.90 0.89 0.92 0.92 0.95 0.91 0.89 0.85 0.91 0.90 0.90 0.86 0.86
## [331] 0.87 0.93 0.95 0.84 0.90 0.92 0.90 0.93 0.92 0.92 0.89 0.91 0.87 0.89 0.90
## [346] 0.88 0.90 0.92 0.92 0.90 0.86 0.88 0.86 0.91 0.95 0.96 0.89 0.95 0.91 0.87
## [361] 0.92 0.91 0.90 0.94 0.85 0.95 0.92 0.87 0.93 0.90 0.88 0.95 0.92 0.93 0.89
## [376] 0.95 0.97 0.88 0.93 0.95 0.88 0.96 0.91 0.93 0.87 0.96 0.96 0.88 0.89 0.88
## [391] 0.88 0.87 0.88 0.91 0.88 0.92 0.88 0.90 0.90 0.92 0.89 0.87 0.94 0.90 0.83
## [406] 0.90 0.93 0.85 0.91 0.91 0.89 0.94 0.91 0.91 0.91 0.93 0.91 0.93 0.91 0.87
## [421] 0.95 0.88 0.90 0.88 0.90 0.93 0.91 0.90 0.90 0.84 0.90 0.86 0.91 0.91 0.86
## [436] 0.91 0.93 0.92 0.93 0.82 0.88 0.90 0.89 0.89 0.93 0.90 0.86 0.88 0.91 0.90
## [451] 0.92 0.90 0.92 0.87 0.91 0.85 0.85 0.85 0.93 0.89 0.91 0.93 0.95 0.92 0.87
## [466] 0.89 0.89 0.92 0.88 0.90 0.90 0.88 0.94 0.87 0.96 0.94 0.89 0.91 0.89 0.90
## [481] 0.85 0.92 0.97 0.94 0.87 0.88 0.87 0.91 0.93 0.90 0.88 0.88 0.93 0.92 0.92
## [496] 0.92 0.90 0.96 0.85 0.89hist(estimadores_p_p90, main = "Histograma Estimadores (500 veces)")
summary(estimadores_p_p90)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.8100 0.8800 0.9000 0.9015 0.9200 0.9800sd(estimadores_p_p90)## [1] 0.0294272#d.
estimadores_p5_90=sapply(rep(5,500), calc_enfermas90)
estimadores_p10_90=sapply(rep(10,500), calc_enfermas90)
estimadores_p15_90=sapply(rep(15,500), calc_enfermas90)
estimadores_p20_90=sapply(rep(20,500), calc_enfermas90)
estimadores_p30_90=sapply(rep(30,500), calc_enfermas90)
estimadores_p50_90=sapply(rep(50,500), calc_enfermas90)
estimadores_p60_90=sapply(rep(60,500), calc_enfermas90)
estimadores_p100_90=sapply(rep(100,500), calc_enfermas90)
estimadores_p200_90=sapply(rep(200,500), calc_enfermas90)
estimadores_p500_90=sapply(rep(500,500), calc_enfermas90)
res_p_90=data.frame(estimadores_p5_90,estimadores_p10_90,estimadores_p15_90,estimadores_p20_90,estimadores_p30_90,estimadores_p50_90,estimadores_p60_90,estimadores_p100_90,estimadores_p200_90,estimadores_p500_90)
# Ajuste para tabla:
indicadores90 <- cbind(apply(res_p_90, 2, mean),
apply(res_p_90, 2, median),
apply(res_p_90, 2, sd),
apply(res_p_90, 2, min),
apply(res_p_90, 2, max))
colnames(indicadores90) <- c("media","mediana","sd","min","max")
round(indicadores90,3)## media mediana sd min max
## estimadores_p5_90 0.905 1.000 0.134 0.400 1.000
## estimadores_p10_90 0.895 0.900 0.101 0.600 1.000
## estimadores_p15_90 0.897 0.933 0.076 0.600 1.000
## estimadores_p20_90 0.897 0.900 0.068 0.650 1.000
## estimadores_p30_90 0.898 0.900 0.055 0.700 1.000
## estimadores_p50_90 0.896 0.900 0.041 0.780 0.980
## estimadores_p60_90 0.901 0.900 0.038 0.783 0.983
## estimadores_p100_90 0.900 0.900 0.029 0.800 0.960
## estimadores_p200_90 0.900 0.900 0.019 0.845 0.950
## estimadores_p500_90 0.900 0.900 0.009 0.872 0.926boxplot(res_p_90)
abline(h=0.9,col="red",lwd=1)
# Pruebas de bondad y ajuste (shapiro wilks)
# Hipótesis nula: las estimaciones provienen de una distribución Normal
# Hipótesis alternativa: Las estimaciones NO provienen de una distribución Normal
# Criterio de rechazo: Dado una significancia (Asumiendo alpha = 0.05)
# NO Se rechaza H_0 si el p-valor de la prueba es > a 0.05
# Se rechaza H_0 si el p-valor de la prueba es <= a 0.05
shapiro.test(estimadores_p5_90)##
## Shapiro-Wilk normality test
##
## data: estimadores_p5_90
## W = 0.69393, p-value < 2.2e-16shapiro.test(estimadores_p10_90)##
## Shapiro-Wilk normality test
##
## data: estimadores_p10_90
## W = 0.84118, p-value < 2.2e-16shapiro.test(estimadores_p15_90)##
## Shapiro-Wilk normality test
##
## data: estimadores_p15_90
## W = 0.90013, p-value < 2.2e-16shapiro.test(estimadores_p20_90)##
## Shapiro-Wilk normality test
##
## data: estimadores_p20_90
## W = 0.92573, p-value = 5.122e-15shapiro.test(estimadores_p30_90)##
## Shapiro-Wilk normality test
##
## data: estimadores_p30_90
## W = 0.95516, p-value = 3.442e-11shapiro.test(estimadores_p50_90)##
## Shapiro-Wilk normality test
##
## data: estimadores_p50_90
## W = 0.97049, p-value = 1.721e-08shapiro.test(estimadores_p60_90)##
## Shapiro-Wilk normality test
##
## data: estimadores_p60_90
## W = 0.9757, p-value = 2.201e-07shapiro.test(estimadores_p100_90)##
## Shapiro-Wilk normality test
##
## data: estimadores_p100_90
## W = 0.9798, p-value = 2.051e-06shapiro.test(estimadores_p200_90)##
## Shapiro-Wilk normality test
##
## data: estimadores_p200_90
## W = 0.99053, p-value = 0.002682shapiro.test(estimadores_p500_90)##
## Shapiro-Wilk normality test
##
## data: estimadores_p500_90
## W = 0.99363, p-value = 0.03342require(ggplot2)## Loading required package: ggplot2
require(ggpubr)## Loading required package: ggpubrge5_90=ggplot(res_p, aes(sample = estimadores_p5_90)) + stat_qq(distribution = qnorm, dparams = 5) + ggtitle("QQ-plot para n=5")+theme_bw()
ge10_90=ggplot(res_p, aes(sample = estimadores_p10_90)) + stat_qq(distribution = qnorm, dparams = 5) + ggtitle("QQ-plot para n=10")+theme_bw()
ge15_90=ggplot(res_p, aes(sample = estimadores_p15_90)) + stat_qq(distribution = qnorm, dparams = 5) + ggtitle("QQ-plot para n=15")+theme_bw()
ge20_90=ggplot(res_p, aes(sample = estimadores_p20_90)) + stat_qq(distribution = qnorm, dparams = 5) + ggtitle("QQ-plot para n=20")+theme_bw()
ge30_90=ggplot(res_p, aes(sample = estimadores_p30_90)) + stat_qq(distribution = qnorm, dparams = 5) + ggtitle("QQ-plot para n=30")+theme_bw()
ge50_90=ggplot(res_p, aes(sample = estimadores_p50_90)) + stat_qq(distribution = qnorm, dparams = 5) + ggtitle("QQ-plot para n=50")+theme_bw()
ge60_90=ggplot(res_p, aes(sample = estimadores_p60_90)) + stat_qq(distribution = qnorm, dparams = 5) + ggtitle("QQ-plot para n=60")+theme_bw()
ge100_90=ggplot(res_p, aes(sample = estimadores_p100_90)) + stat_qq(distribution = qnorm, dparams = 5) + ggtitle("QQ-plot para n=100")+theme_bw()
ge200_90=ggplot(res_p, aes(sample = estimadores_p200_90)) + stat_qq(distribution = qnorm, dparams = 5) + ggtitle("QQ-plot para n=200")+theme_bw()
ge500_90=ggplot(res_p, aes(sample = estimadores_p500_90)) + stat_qq(distribution = qnorm, dparams = 5) + ggtitle("QQ-plot para n=500")+theme_bw()
ggarrange(ge5_90, ge10_90, ge15_90, ge20_90, ge30_90, ge50_90, ge60_90, ge100_90, ge200_90, ge500_90, ncol = 2, nrow = 3)## $`1`
##
## $`2`
##
## attr(,"class")
## [1] "list" "ggarrange"Al realizar las pruebas de bondad y ajuste (shapiro wilks) para comprobar la normalidad de los lotes creados con 10% y 90% de plantas enfermas y con los diferentes tamaños de muestra, se obtiene como resultado para todos los casos que se rechaza la hipotesis nula: las estimaciones provienen de una distribución Normal, dado que en todas las pruebas se obtuvo un p-valor menor que la significancia cuyo valor es del 5%, esto nos indica entonces que las muestras NO tienen una distribución normal.
PUNTO 2
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 utilizara 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).
N1 = c(rep("Plantas_Buenas",900), rep("Plantas_Enfermas",100))
N2 = c(rep("Plantas_Buenas",1350), rep("Plantas_Enfermas",150))
table(N1)/1000## N1
## Plantas_Buenas Plantas_Enfermas
## 0.9 0.1table(N2)/1500## N2
## Plantas_Buenas Plantas_Enfermas
## 0.9 0.1b. 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=function(n){
muestraN1=sample(N1,size = n)
muestraN2=sample(N2,size = n)
p1=sum(muestraN1=="Plantas_Enfermas")/n
p2=sum(muestraN2=="Plantas_Enfermas")/n
y = p1-p2
return(y)
}
calc_dif(n=100)## [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?.
par(mfrow=c(1,2))
est_dif= sapply(rep(100,500), calc_dif)
hist(est_dif, main = "Hist. 500 estimadores de p1-p2")
boxplot(est_dif, main = "G. de cajas 500 est. de p1-p2")
abline(h=0,col="red",lwd=1)
summary(est_dif)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.1200 -0.0200 0.0000 0.0011 0.0200 0.1600sd(est_dif)## [1] 0.03877625De acuerdo con los graficos de histograma y de cajas, se evidencia que con respecto a las diferencias de los 500 estimadores de p1-p2 se observa que la concentración de los datos se encuentran alrededor del cero, no siempre la diferencia entre el numero de plantas enfermas de ambas muestras es cero. Presenta una desviación estandar relativamente pequeña, igualmente la media y la mediana tienden al valor de cero, esto nos permite inferir que tiene una distribución simétrica.
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 evalué si. ¿Considera que es mas probable concluir que existen diferencias entre los tratamientos con muestras grandes que pequeñas, es decir cual considera usted que es el efecto del tamaño de muestra en el caso de la comparación de proporciones?.
est_dif5= sapply(rep(5,500), calc_dif)
est_dif10= sapply(rep(10,500), calc_dif)
est_dif15= sapply(rep(15,500), calc_dif)
est_dif20= sapply(rep(20,500), calc_dif)
est_dif30= sapply(rep(30,500), calc_dif)
est_dif50= sapply(rep(50,500), calc_dif)
est_dif60= sapply(rep(60,500), calc_dif)
est_dif100= sapply(rep(100,500), calc_dif)
est_dif200= sapply(rep(200,500), calc_dif)
est_dif500= sapply(rep(500,500), calc_dif)
res_dif=data.frame(est_dif5, est_dif10, est_dif15, est_dif20, est_dif30, est_dif50, est_dif60, est_dif100, est_dif200, est_dif500)
boxplot(res_dif)
abline(h=0, col="red")
De acuerdo con el diagrama de cajas, se evidencia que entre mayor es el tamaño de la muestra, menor es la dispersión de los datos, se percibe que los datos son simetricos, presentan una mediana de 0 y la agrupación del 50% de los datos tienden a 0.
Pruebas de normalidad de shapiro wilk:
Hipótesis nula: Las estimaciones provienen de una distribución Normal
Hipótesis alternativa: Las estimaciones NO provienen de una distribución Normal
Criterio de rechazo: Dado una significancia (Asumiendo alpha = 0.05)
NO Se rechaza H_0 si el p-valor de la prueba es > a 0.05
Se rechaza H_0 si el p-valor de la prueba es <= a 0.05
shapiro.test(est_dif5)##
## Shapiro-Wilk normality test
##
## data: est_dif5
## W = 0.90602, p-value < 2.2e-16shapiro.test(est_dif10)##
## Shapiro-Wilk normality test
##
## data: est_dif10
## W = 0.95092, p-value = 7.905e-12shapiro.test(est_dif15)##
## Shapiro-Wilk normality test
##
## data: est_dif15
## W = 0.96864, p-value = 7.421e-09shapiro.test(est_dif20)##
## Shapiro-Wilk normality test
##
## data: est_dif20
## W = 0.97678, p-value = 3.883e-07shapiro.test(est_dif30)##
## Shapiro-Wilk normality test
##
## data: est_dif30
## W = 0.97792, p-value = 7.188e-07shapiro.test(est_dif50)##
## Shapiro-Wilk normality test
##
## data: est_dif50
## W = 0.98785, p-value = 0.0003576shapiro.test(est_dif60)##
## Shapiro-Wilk normality test
##
## data: est_dif60
## W = 0.99076, p-value = 0.003211shapiro.test(est_dif100)##
## Shapiro-Wilk normality test
##
## data: est_dif100
## W = 0.99145, p-value = 0.005561shapiro.test(est_dif200)##
## Shapiro-Wilk normality test
##
## data: est_dif200
## W = 0.99173, p-value = 0.006947shapiro.test(est_dif500)##
## Shapiro-Wilk normality test
##
## data: est_dif500
## W = 0.99536, p-value = 0.1424par(mfrow=c(2,2))
qqnorm(est_dif5)
qqline(est_dif5, col="red")
qqnorm(est_dif10)
qqline(est_dif10, col="red")
qqnorm(est_dif15)
qqline(est_dif15, col="red")
qqnorm(est_dif20)
qqline(est_dif20, col="red")
qqnorm(est_dif30)
qqline(est_dif30, col="red")
qqnorm(est_dif50)
qqline(est_dif50, col="red")
qqnorm(est_dif60)
qqline(est_dif60, col="red")
qqnorm(est_dif100)
qqline(est_dif100, col="red")
qqnorm(est_dif200)
qqline(est_dif200, col="red")
qqnorm(est_dif500)
qqline(est_dif500, col="red")
La prueba de normalidad de Shapiro Wilk tiene mejor aplicabilidad cuando se analizan tamaños de muestras pequeñas. De acuerdo con los ejercicios realizados, parece ser que una muestra ideal para este ejercicio de de tamaño 200.
En cuanto a la normalidad, se rechaza la hipotesis nula en las muestras pequeñas, se aprueba la hipotesis nula para las muestras mas grandes, en los graficos qq se observa que va normalizandose a medida que crece la muestra, por el contrario para las muestras pequeñas se evidencia sesgo en los datos.
En terminos de proporción para ambas muestras de los lotes de plantas, entre mayor es el tamaño de las muestras, se concentran mas los datos hacia una diferencia de 0, es decir que el porcentaje de plantas enfermas de la muestra del lote 1 tiende a ser igual al del lote 2, por el contrario, si se toma la decisión sobre una muestra pequeña, puede caer en el error de inferir que el desempeño de un lote es meor que el del otro, cuando realmente el desmpeño es el mismo.
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 presento 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.
Na = c(rep("Plantas_Buenas",900), rep("Plantas_Enfermas",100))
Nb = c(rep("Plantas_Buenas",1275), rep("Plantas_Enfermas",225))
table(Na)/1000## Na
## Plantas_Buenas Plantas_Enfermas
## 0.9 0.1table(Nb)/1500## Nb
## Plantas_Buenas Plantas_Enfermas
## 0.85 0.15#b.
calc_dif_1=function(n){
muestra_Na=sample(Na,size = n)
muestra_Nb=sample(Nb,size = n)
p_1=sum(muestra_Na=="Plantas_Enfermas")/n
p_2=sum(muestra_Nb=="Plantas_Enfermas")/n
y = p_1-p_2
return(y)
}
calc_dif_1(n=100)## [1] -0.13# c.
par(mfrow=c(1,2))
est_dif_a= sapply(rep(100,500), calc_dif_1)
hist(est_dif_a, main = "Hist. 500 estimadores de p1-p2")
boxplot(est_dif_a, main = "G. de cajas 500 est. de p1-p2")
abline(h=-0.05,col="red",lwd=1)
summary(est_dif_a)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.19000 -0.08000 -0.05000 -0.04952 -0.02000 0.08000sd(est_dif_a)## [1] 0.04620432# d.
est_dif5_a= sapply(rep(5,500), calc_dif_1)
est_dif10_a= sapply(rep(10,500), calc_dif_1)
est_dif15_a= sapply(rep(15,500), calc_dif_1)
est_dif20_a= sapply(rep(20,500), calc_dif_1)
est_dif30_a= sapply(rep(30,500), calc_dif_1)
est_dif50_a= sapply(rep(50,500), calc_dif_1)
est_dif60_a= sapply(rep(60,500), calc_dif_1)
est_dif100_a= sapply(rep(100,500), calc_dif_1)
est_dif200_a= sapply(rep(200,500), calc_dif_1)
est_dif500_a= sapply(rep(500,500), calc_dif_1)
res_dif_a=data.frame(est_dif5_a, est_dif10_a, est_dif15_a, est_dif20_a, est_dif30_a, est_dif50_a, est_dif60_a, est_dif100_a, est_dif200_a, est_dif500_a)
boxplot(res_dif_a)
abline(h=-0.05, col="red")
shapiro.test(est_dif5_a)##
## Shapiro-Wilk normality test
##
## data: est_dif5_a
## W = 0.91006, p-value < 2.2e-16shapiro.test(est_dif10_a)##
## Shapiro-Wilk normality test
##
## data: est_dif10_a
## W = 0.96202, p-value = 4.56e-10shapiro.test(est_dif15_a)##
## Shapiro-Wilk normality test
##
## data: est_dif15_a
## W = 0.97259, p-value = 4.651e-08shapiro.test(est_dif20_a)##
## Shapiro-Wilk normality test
##
## data: est_dif20_a
## W = 0.97555, p-value = 2.037e-07shapiro.test(est_dif30_a)##
## Shapiro-Wilk normality test
##
## data: est_dif30_a
## W = 0.98388, p-value = 2.41e-05shapiro.test(est_dif50_a)##
## Shapiro-Wilk normality test
##
## data: est_dif50_a
## W = 0.98839, p-value = 0.0005268shapiro.test(est_dif60_a)##
## Shapiro-Wilk normality test
##
## data: est_dif60_a
## W = 0.99133, p-value = 0.005028shapiro.test(est_dif100_a)##
## Shapiro-Wilk normality test
##
## data: est_dif100_a
## W = 0.99221, p-value = 0.01025shapiro.test(est_dif200_a)##
## Shapiro-Wilk normality test
##
## data: est_dif200_a
## W = 0.99679, p-value = 0.4267shapiro.test(est_dif500_a)##
## Shapiro-Wilk normality test
##
## data: est_dif500_a
## W = 0.99486, p-value = 0.09371par(mfrow=c(2,2))
qqnorm(est_dif5_a)
qqline(est_dif5_a, col="red")
qqnorm(est_dif10_a)
qqline(est_dif10_a, col="red")
qqnorm(est_dif15_a)
qqline(est_dif15_a, col="red")
qqnorm(est_dif20_a)
qqline(est_dif20_a, col="red")
qqnorm(est_dif30_a)
qqline(est_dif30_a, col="red")
qqnorm(est_dif50_a)
qqline(est_dif50_a, col="red")
qqnorm(est_dif60_a)
qqline(est_dif60_a, col="red")
qqnorm(est_dif100_a)
qqline(est_dif100_a, col="red")
qqnorm(est_dif200_a)
qqline(est_dif200_a, col="red")
qqnorm(est_dif500_a)
qqline(est_dif500_a, col="red")
Al realizar el ejercicio con proporciones diferentes de plantas
enfermas, se observa que al obtener la diferncia entre las muestras de
cada población p1-p2, tiende a ser negativo el valor, pues vemos un
promedio y una mediana cercanos al -0.05 con una baja desviación
estandar, permitiendo inferir para el escenario 2 que se escoge el lote
1 de 1000 plantas con el 10% de plantas enfermas, dado que presenta un
mejor desempeño en el control de la plaga, mientras que en el escenario
1 ambos lotes presentan el mismo desempeño frente al control de la
plaga.
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” escriba un resumen (máximo 2 paginas) sobre el artículos e incluya en este sus opiniones en cuanto al uso del valor p.
Ensayo John Jairo Bedoya:
La evolución de la estadística y la prueba P Valor
El P valor fue creado por el estadístico Ronald Fisher como una forma informal de juzgar si la evidencia era significativa, no pretendía que se convirtiera en una prueba definitiva. Aunque los valores p han tenido críticas, siempre necesitaremos de ellos, permite que la toma de decisiones basadas en evidencia sea más rigurosa y objetiva posible.
El valor p evolucionó y los expertos en el tema no lo previeron solo hasta que se dieron cuenta que los resultados obtenidos en ciertos estudios no eran confiables. En el mundo, la tecnología ha hecho que todos los sectores de la economía avancen rápidamente, pero los estadísticos, aunque cada vez utilizan herramientas más sofisticadas para sus análisis, algunos criterios estadísticos como el p valor no ha tenido mayor avance en nueve décadas y a su vez ha sido sobrevalorado.
Se ha sobrevalorado esta herramienta estadística para el cálculo de probabilidades de ocurrencia pues este no puede trabajar hacia atrás y hacer afirmaciones sobre la realidad subyacente, lo que puede hacer es resumir datos partiendo de una hipótesis nula específica. Para calcular probabilidades se requiere contar con otra información u herramientas que permitan hacer un análisis más amplio de la situación, de lo contrario podríamos inferir por ejemplo que una persona tiene un tumor cerebral partiendo solo de un síntoma de dolor de cabeza. Los cambios en algunas decisiones de análisis de datos han demostrado que pueden aumentar la tasa de falsos positivos en un solo estudio al 60%.
Algunos académicos piensan que se necesita cambiar la forma como se enseña la estadística, cómo se realiza el análisis de datos y cómo se informan y se interpretan los resultados. Uno de los motivos a los que llevó esto es que el p valor nunca se tenía que utilizar de la forma en que hoy se usa, pues se tiende a caer en la trampa de pensar en los resultados como significativos y no significativos. Según el estadístico Richard Royall un científico debería responder estas preguntas al hacer un estudio, ¿cuál es la evidencia?, ¿qué debo creer? y ¿qué debo hacer?, para esto un solo método no es suficiente para dar respuesta, los números son donde la discusión científica debe comenzar, no terminar.
A pesar que el p valor es una prueba de mucha validez estadística, es inexacta en algunos casos, es el momento de avanzar en la aplicación de nuevos métodos estadísticos que cambien el paradigma y complementen estos resultados.
Ensayo Liz Gutierrez:
La prueba P Valor no es tan fiable como se cree
La prueba p valor lleva nueve décadas desde que el científico Ronald Fisher la creó utilizándola de forma informal para sus pruebas, nunca pensó que se convirtiera en una prueba definitiva. Hoy parece tener la última palabra, pero no es así, se necesita usar otros criterios como el intervalo de confianza y formulaciones de pruebas de hipótesis que proporcionen más profundidad en el análisis de los resultados obtenidos. El “estándar de oro” de la validez estadística, no es tan confiable como se supone.
Las pruebas de hipótesis de estudios que han sido hitos históricos, al evaluarlas hoy presentan inconsistencias teóricas en la técnica y de uso incorrecto de los procedimientos. En el caso de la biomedicina, refleja cierto subjetivismo en el análisis de datos relacionado con la plausibilidad biológica, a lo que se le suma un manejo deficiente de la información que arrojan los intervalos de confianza.
Estamos en la era del paso de las cosas a las no-cosas, estás últimas se refieren a informaciones, nuestra concentración está puesta en la información y en los datos, cada vez hay volúmenes de datos mucho más grandes y complejos, pero hemos visto como avanza la tecnología y las herramientas que usan los académicos pero el p valor no hizo esta misma evolución. Por esto se hace cada vez más necesario avanzar en los métodos estadísticos como el p valor o alguno en su reemplazo, para realizar análisis más precisos que permitan la toma de decisiones acertadas.
Un error muy común en el uso de esta prueba es dicotomizar la decisión, si se consulta la teoría, la recomendación es utilizar un valor de 0,05 como valor de umbral y enfocar el análisis de hipótesis si p< 0,05 se rechaza la hipótesis nula y en caso contrario se acepta, lo cual es incorrecto ya que si un estudio tiene un p igual a 0,04999 rechaza la hipótesis nula pero si el resultado es p igual a 0,05001 se acepta la hipótesis nula porque p es mayor que 0,05 lo cual muestra una clara inconsistencia pues la diferencia entre los dos valores es de 0,00002 lo cual es muy pequeña, sin embargo llevan a dos conclusiones diferentes. Lo que aportaría en este caso es mostrar el valor p en la cuantía que se obtuvo y analizarlo.
Por casos como el anterior algunos académicos manifiestan que debería cambiarse la forma en la que se analizan los datos, se recogen los resultados y la manera de enseñar la estadística misma. Al adelantar una investigación deberíamos comenzarla con los números más no debería ser al terminar como se suele hacer. Se deben aplicar varios métodos que en su conjunto lleven a dar respuesta a tres preguntas: ¿cuál es la evidencia?, ¿qué debo creer? y ¿qué debo hacer?. De esta manera se podrá poco a poco cambiar el paradigma de los estadísticos actuales y evitar el sesgo de los futuros.