Inferencia Estadistica y Simulacion
Yeffer Edilberto Diaz Cuesta
2022-09-12
Taller Inferencia Estadística y Simulación
1. El Teorema del Limite Central es uno de los mas importantes en la inferencia estadística y habla sobre la convergencia de los estimadores como la proporción muestral a la distribución normal. Algunos autores afirman que esta aproximación es bastante buena a partir del umbral n>30.
1. Desarrollo primer punto
SIMULACIÓN CON POBLACION ENFERMA DEL 50%
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%.
lotePlantas=c(rep("sanas",500),rep("enfermas",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.
plan_aleat=function(n){
muestra=sample(lotePlantas, size=n)
return(sum(muestra=="enfermas")/n)
}
n=10
plan_aleat(n)## [1] 0.4estimadores=sapply(rep(n,500), plan_aleat)
hist(estimadores,col = "lightblue")
mean(estimadores)## [1] 0.494sd(estimadores)## [1] 0.1585379c. 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=sapply(rep(500,500), plan_aleat)
hist(estimadores,col = "lightgoldenrod")
mean(estimadores)## [1] 0.500348sd(estimadores)## [1] 0.01609438Respuesta. En términos de simetría los histogramas muestran que en los dos casos tiene un comportamiento normal, sin embargo, cuando las muestras realizadas son mayores por ejemplo el de 500 estimadores, se visualiza que el conjunto de datos ese acumula sobre el valor medio real del porcentaje de plantas enfermas. Esto también es evidente en términos de la variabilidad, en donde se pasa de tener aproximadamente un 16% cuando se realiza solamente 10 muestras a un 1.6% cuando se realizaron las 50 muestras.
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).
est_5=sapply(rep(5,500), plan_aleat)
est_10=sapply(rep(10,500), plan_aleat)
est_15=sapply(rep(15,500), plan_aleat)
est_20=sapply(rep(20,500), plan_aleat)
est_30=sapply(rep(30,500), plan_aleat)
est_50=sapply(rep(50,500), plan_aleat)
est_60=sapply(rep(60,500), plan_aleat)
est_100=sapply(rep(100,500), plan_aleat)
est_200=sapply(rep(200,500), plan_aleat)
est_500=sapply(rep(500,500), plan_aleat)
res=data.frame(est_5,est_10,est_15,est_20,est_30,est_50,est_60,est_100,est_200,est_500)
sapply(res, mean)## est_5 est_10 est_15 est_20 est_30 est_50 est_60 est_100
## 0.4992000 0.4954000 0.5040000 0.5094000 0.4997333 0.4983600 0.4992333 0.5006200
## est_200 est_500
## 0.4997400 0.5007400sapply(res, sd)## est_5 est_10 est_15 est_20 est_30 est_50 est_60
## 0.22490114 0.14615595 0.12438940 0.10707273 0.09071781 0.06621540 0.06775145
## est_100 est_200 est_500
## 0.04788320 0.02973214 0.01627729sapply(res,shapiro.test)## est_5 est_10
## statistic 0.9290194 0.9580274
## p.value 1.199336e-14 9.811173e-11
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## est_15 est_20
## statistic 0.9700955 0.9772095
## p.value 1.433704e-08 4.882931e-07
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## est_30 est_50
## statistic 0.9864278 0.9888612
## p.value 0.0001309478 0.00074821
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## est_60 est_100
## statistic 0.9901146 0.9943122
## p.value 0.001936034 0.05941312
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## est_200 est_500
## statistic 0.9950524 0.9965906
## p.value 0.1103715 0.3708957
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"require(ggplot2)## Loading required package: ggplot2require (ggpubr)## Loading required package: ggpubrg1=ggplot(res,aes(x=est_5))+geom_histogram(fill = "lightblue",col = 1)+geom_density(col = 2,lwd=0.8)+xlab("Muestra n=5")+theme_bw()
g2=ggplot(res,aes(x=est_5))+geom_boxplot(fill = "lightblue")+xlab("Muestra n=5")+theme_bw()
g3=ggplot(res,aes(x=est_10))+geom_histogram(fill = "lightblue",col = 1)+geom_density(col = 2,lwd=0.8)+xlab("Muestra n=10")+theme_bw()
g4=ggplot(res,aes(x=est_10))+geom_boxplot(fill = "lightblue")+xlab("Muestra n=10")+theme_bw()
g5=ggplot(res,aes(x=est_15))+geom_histogram(fill = "lightgreen",col = 1)+geom_density(col = 2,lwd=0.8)+xlab("Muestra n=15")+theme_bw()
g6=ggplot(res,aes(x=est_15))+geom_boxplot(fill = "lightgreen")+xlab("Muestra n=15")+theme_bw()
g7=ggplot(res,aes(x=est_20))+geom_histogram(fill = "lightgreen",col = 1)+geom_density(col = 2,lwd=0.8)+xlab("Muestra n=20")+theme_bw()
g8=ggplot(res,aes(x=est_20))+geom_boxplot(fill = "lightgreen")+xlab("Muestra n=20")+theme_bw()
g9=ggplot(res,aes(x=est_30))+geom_histogram(fill = "wheat",col = 1)+geom_density(col = 2,lwd=0.8)+xlab("Muestra n=30")+theme_bw()
g10=ggplot(res,aes(x=est_30))+geom_boxplot(fill = "wheat")+xlab("Muestra n=30")+theme_bw()
g11=ggplot(res,aes(x=est_50))+geom_histogram(fill = "wheat",col = 1)+geom_density(col = 2,lwd=0.8)+xlab("Muestra n=50")+theme_bw()
g12=ggplot(res,aes(x=est_50))+geom_boxplot(fill = "wheat")+xlab("Muestra n=50")+theme_bw()
g13=ggplot(res,aes(x=est_60))+geom_histogram(fill = "ivory",col = 1)+geom_density(col = 2,lwd=0.8)+xlab("Muestra n=60")+theme_bw()
g14=ggplot(res,aes(x=est_60))+geom_boxplot(fill = "ivory")+xlab("Muestra n=60")+theme_bw()
g15=ggplot(res,aes(x=est_100))+geom_histogram(fill = "ivory",col = 1)+geom_density(col = 2,lwd=0.8)+xlab("Muestra n=100")+theme_bw()
g16=ggplot(res,aes(x=est_100))+geom_boxplot(fill = "ivory")+xlab("Muestra n=100")+theme_bw()
g17=ggplot(res,aes(x=est_200))+geom_histogram(fill = "olivedrab1",col = 1)+geom_density(col = 2,lwd=0.8)+xlab("Muestra n=200")+theme_bw()
g18=ggplot(res,aes(x=est_200))+geom_boxplot(fill = "olivedrab1")+xlab("Muestra n=200")+theme_bw()
g19=ggplot(res,aes(x=est_500))+geom_histogram(fill = "olivedrab1",col = 1)+geom_density(col = 2,lwd=0.8)+xlab("Muestra n=500")+theme_bw()
g20=ggplot(res,aes(x=est_500))+geom_boxplot(fill = "olivedrab1")+xlab("Muestra n=500")+theme_bw()
GRAFICOS DE NORMALIDAD QQ
qqnorm(est_5,main='Grafico QQ de Normalidad estimador 5',col="lightcoral")
qqline(est_5)
qqnorm(est_50,main='Grafico QQ de Normalidad estimador 50',col="lightcoral")
qqline(est_50)
qqnorm(est_100,main='Grafico QQ de Normalidad estimador 100',col="lightcoral")
qqline(est_100)
qqnorm(est_500,main='Grafico QQ de Normalidad estimador 500',col="lightcoral")
qqline(est_500)
Respuesta. En términos de porcentaje se observa que cuando es mayor el número de muestras realizadas el resultado se acerca cada vez más al valor real (50%) de la población simulada en este caso pasando del 51.5% del estimador 5 al 49.9% del estimador 500, en términos de variabilidad también se evidencia que entre mas muestras los valores obtenidos son menos dispersos con menor coeficiente de desviación estándar. Por último, los valores obtenidos con el test de normalidad muestran que los valores con una confianza superior en su valor p al 5% son el de los estimadores 200 y 500, los demás valores serian rechazados.
e. Repita toda la simulación (puntos a – d) pero ahora con lotes con 10% y 90% de plantas enfermas. Concluya todo el ejercicio. :::
SIMULACIÓN CON POBLACION ENFERMA DEL 10%
lotePlantas=c(rep("sanas",900),rep("enfermas",100))
h=10
plan_aleat(h)## [1] 0.1estimadores=sapply(rep(h,500), plan_aleat)
hist(estimadores,col = "lightblue")
mean(estimadores)## [1] 0.1024sd(estimadores)## [1] 0.09344358estimadores=sapply(rep(500,500), plan_aleat)
hist(estimadores,col = "lightgoldenrod")
mean(estimadores)## [1] 0.100204sd(estimadores)## [1] 0.009189956Tamaños de muestra n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500. Y comparacion de 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).
est_5=sapply(rep(5,500), plan_aleat)
est_10=sapply(rep(10,500), plan_aleat)
est_15=sapply(rep(15,500), plan_aleat)
est_20=sapply(rep(20,500), plan_aleat)
est_30=sapply(rep(30,500), plan_aleat)
est_50=sapply(rep(50,500), plan_aleat)
est_60=sapply(rep(60,500), plan_aleat)
est_100=sapply(rep(100,500), plan_aleat)
est_200=sapply(rep(200,500), plan_aleat)
est_500=sapply(rep(500,500), plan_aleat)
res=data.frame(est_5,est_10,est_15,est_20,est_30,est_50,est_60,est_100,est_200,est_500)
sapply(res, mean)## est_5 est_10 est_15 est_20 est_30 est_50 est_60 est_100
## 0.0956000 0.1032000 0.1025333 0.1106000 0.1006000 0.1014800 0.1000000 0.1008600
## est_200 est_500
## 0.0996400 0.1003960sapply(res, sd)## est_5 est_10 est_15 est_20 est_30 est_50
## 0.134822083 0.095959243 0.080095579 0.066982186 0.055489923 0.042716026
## est_60 est_100 est_200 est_500
## 0.036884940 0.027264392 0.018358771 0.009725449sapply(res,shapiro.test)## est_5 est_10
## statistic 0.6894479 0.8446211
## p.value 2.697672e-29 9.484268e-22
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## est_15 est_20
## statistic 0.8962517 0.9337445
## p.value 6.732848e-18 4.273184e-14
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## est_30 est_50
## statistic 0.9507338 0.9675273
## p.value 7.42223e-12 4.531388e-09
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## est_60 est_100
## statistic 0.9808698 0.9864886
## p.value 3.812614e-06 0.0001365544
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## est_200 est_500
## statistic 0.989879 0.9949834
## p.value 0.001614596 0.1042235
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"GRAFICOS DE NORMALIDAD QQ
qqnorm(est_5,main='Grafico QQ de Normalidad estimador 5',col="lightcoral")
qqline(est_5)
qqnorm(est_50,main='Grafico QQ de Normalidad estimador 50',col="lightcoral")
qqline(est_50)
qqnorm(est_100,main='Grafico QQ de Normalidad estimador 100',col="lightcoral")
qqline(est_100)
qqnorm(est_500,main='Grafico QQ de Normalidad estimador 500',col="lightcoral")
qqline(est_500)
SIMULACIÓN CON POBLACION ENFERMA DEL 90%
lotePlantas=c(rep("sanas",100),rep("enfermas",900))
h=10
plan_aleat(h)## [1] 1estimadores=sapply(rep(h,500), plan_aleat)
hist(estimadores,col = "lightblue")
mean(estimadores)## [1] 0.896sd(estimadores)## [1] 0.1014129estimadores=sapply(rep(500,500), plan_aleat)
hist(estimadores,col = "lightgoldenrod")
mean(estimadores)## [1] 0.899944sd(estimadores)## [1] 0.009854089Tamaños de muestra n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500. Y comparacion de 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).
est_5=sapply(rep(5,500), plan_aleat)
est_10=sapply(rep(10,500), plan_aleat)
est_15=sapply(rep(15,500), plan_aleat)
est_20=sapply(rep(20,500), plan_aleat)
est_30=sapply(rep(30,500), plan_aleat)
est_50=sapply(rep(50,500), plan_aleat)
est_60=sapply(rep(60,500), plan_aleat)
est_100=sapply(rep(100,500), plan_aleat)
est_200=sapply(rep(200,500), plan_aleat)
est_500=sapply(rep(500,500), plan_aleat)
res=data.frame(est_5,est_10,est_15,est_20,est_30,est_50,est_60,est_100,est_200,est_500)
sapply(res, mean)## est_5 est_10 est_15 est_20 est_30 est_50 est_60 est_100
## 0.9012000 0.8996000 0.9052000 0.9015000 0.8984667 0.9002400 0.9005000 0.9030200
## est_200 est_500
## 0.8996800 0.8996440sapply(res, sd)## est_5 est_10 est_15 est_20 est_30 est_50
## 0.133093907 0.096011856 0.073948233 0.068132474 0.053676709 0.042031842
## est_60 est_100 est_200 est_500
## 0.038686330 0.027455185 0.018332205 0.008982351sapply(res,shapiro.test)## est_5 est_10
## statistic 0.705719 0.841678
## p.value 1.148656e-28 6.141099e-22
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## est_15 est_20
## statistic 0.8916457 0.9277094
## p.value 2.712873e-18 8.521129e-15
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## est_30 est_50
## statistic 0.9563286 0.9729975
## p.value 5.247464e-11 5.665959e-08
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## est_60 est_100
## statistic 0.9695755 0.9877146
## p.value 1.130218e-08 0.0003237556
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## est_200 est_500
## statistic 0.9897397 0.9935501
## p.value 0.001451216 0.03130384
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"GRAFICOS DE NORMALIDAD QQ
qqnorm(est_5,main='Grafico QQ de Normalidad estimador 5',col="lightcoral")
qqline(est_5)
qqnorm(est_50,main='Grafico QQ de Normalidad estimador 50',col="lightcoral")
qqline(est_50)
qqnorm(est_100,main='Grafico QQ de Normalidad estimador 100',col="lightcoral")
qqline(est_100)
qqnorm(est_500,main='Grafico QQ de Normalidad estimador 500',col="lightcoral")
qqline(est_500)
Concluciones. Al igual que el anterior, se observa que cuando es mayor el número de muestras realizadas el resultado se acerca cada vez más al valor real de la población simulada y a su vez entre más aumenta el estimador la variabilidad en los resultados va a ser menor. Sin embargo, en las pruebas de normalidad realizadas a las poblaciones con muestras del 90% y 10% ninguno de los estimadores seria aprobado, esto se puede deber principalmente a que el porcentaje de plantas enfermas es del 10% y 90% respectivamente, es decir el valor es muy pequeño o muy grande lo que da un margen de combinación entre las muestras que no permite saber con severidad si se trata de un error o del valor real. Los gráficos de QQ normalidad en el estimador mas bajo muestra en principio una gran dispersión de datos cerca del valor 0 de la grafica.
2. La comparación de tratamientos es una practica 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.
2. Desarrollo segundo punto
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).
lotePlantasN1=c(rep("sanas",900),rep("enfermas",100))
lotePlantasN2=c(rep("sanas",1350),rep("enfermas",150))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.
plan_aleat_N1=function(i){
muestra=sample(lotePlantasN1, size=i)
return(sum(muestra=="enfermas")/i)
}
plan_aleat_N2=function(j){
muestra=sample(lotePlantasN2, size=j)
return(sum(muestra=="enfermas")/j)
}
i=j=72
plan_aleat_N1(i)## [1] 0.09722222plan_aleat_N2(j)## [1] 0.125estimadoresN1=sapply(rep(i,200), plan_aleat_N1)
estimadoresN2=sapply(rep(j,200), plan_aleat_N2)
resu_segundoP=data.frame(estimadoresN1,estimadoresN2)
sapply(resu_segundoP, mean)## estimadoresN1 estimadoresN2
## 0.10187500 0.09986111sapply(resu_segundoP, sd)## estimadoresN1 estimadoresN2
## 0.03114307 0.03313704g21=ggplot(resu_segundoP,aes(x=estimadoresN1))+geom_histogram(fill = "lightgreen",col = 1)+geom_density(col = 2,lwd=0.8)+xlab(paste("Muestra P1 n=",i))+theme_bw()
g22=ggplot(resu_segundoP,aes(x=estimadoresN2))+geom_histogram(fill = "lightgreen",col = 1)+geom_density(col = 2,lwd=0.8)+xlab(paste("Muestra P2 n=",j))+theme_bw()
ggarrange(g21,g22,ncol = 2, nrow = 1)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
DIFERRENCIA DE ESTIMADORES
difer=estimadoresN1-estimadoresN2
mean(difer)## [1] 0.002013889c. 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?.
estimadoresN1=sapply(rep(500,500), plan_aleat_N1)
estimadoresN2=sapply(rep(500,500), plan_aleat_N2)
difer=estimadoresN1-estimadoresN2
resu_segundoP=data.frame(estimadoresN1,estimadoresN2,difer)
sapply(resu_segundoP, mean)## estimadoresN1 estimadoresN2 difer
## 0.100308 0.098860 0.001448sapply(resu_segundoP, sd)## estimadoresN1 estimadoresN2 difer
## 0.009451974 0.011140405 0.014759505g23=ggplot(resu_segundoP,aes(x=estimadoresN1))+geom_histogram(fill = "lightgreen",col = 1)+geom_density(col = 2,lwd=0.8)+xlab("Muestra n=500 P1")+theme_bw()
g24=ggplot(resu_segundoP,aes(x=estimadoresN2))+geom_histogram(fill = "lightgreen",col = 1)+geom_density(col = 2,lwd=0.8)+xlab("Muestra n=500 P2")+theme_bw()
g25=ggplot(resu_segundoP,aes(x=difer))+geom_histogram(fill = "lightgreen",col = 1)+geom_density(col = 2,lwd=0.8)+xlab("Muestra n=500 P1-P2")+theme_bw()
ggarrange(g23,g24,g25,ncol = 3, nrow = 1)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
a=mean(difer)
b=sd(difer)
print(paste("La media de P1-P2 es= ",format(a, scientific = F)))## [1] "La media de P1-P2 es= 0.001448"print(paste("La desviacion estandar de P1-P2 es= ",format(b, scientific = F)))## [1] "La desviacion estandar de P1-P2 es= 0.0147595"Respuesta. En términos de simetría las dos poblaciones N1 y N2 conservan un comportamiento normal, y su media gira entorno al valor real del 10% en los dos caso, sin embargo se observa que en la muestra N1 con población de 1000 individuos la variabilidad de los datos es menor rondando el 0.9%, esto se puede deber principalmente a que al tener una menor población al momento de hacer las muestras es muy probable que el valor se acerque al real ya que las combinaciones que se pueden llegar a presentar son más pequeñas en comparación a la población N2 con 1500 individuos. Con relación a la diferencia de P1-P2 aunque los valores se acumulan entorno a 0, la media de este aunque se acerca nunca se vuelve 0, por más que se repita el ejercicio.
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_N1_5=sapply(rep(5,500), plan_aleat_N1)
est_N1_10=sapply(rep(10,500), plan_aleat_N1)
est_N1_15=sapply(rep(15,500), plan_aleat_N1)
est_N1_20=sapply(rep(20,500), plan_aleat_N1)
est_N1_30=sapply(rep(30,500), plan_aleat_N1)
est_N1_50=sapply(rep(50,500), plan_aleat_N1)
est_N1_60=sapply(rep(60,500), plan_aleat_N1)
est_N1_100=sapply(rep(100,500), plan_aleat_N1)
est_N1_200=sapply(rep(200,500), plan_aleat_N1)
est_N1_500=sapply(rep(500,500), plan_aleat_N1)
est_N2_5=sapply(rep(5,500), plan_aleat_N2)
est_N2_10=sapply(rep(10,500), plan_aleat_N2)
est_N2_15=sapply(rep(15,500), plan_aleat_N2)
est_N2_20=sapply(rep(20,500), plan_aleat_N2)
est_N2_30=sapply(rep(30,500), plan_aleat_N2)
est_N2_50=sapply(rep(50,500), plan_aleat_N2)
est_N2_60=sapply(rep(60,500), plan_aleat_N2)
est_N2_100=sapply(rep(100,500), plan_aleat_N2)
est_N2_200=sapply(rep(200,500), plan_aleat_N2)
est_N2_500=sapply(rep(500,500), plan_aleat_N2)
difer_5=est_N1_5 - est_N2_5
difer_10=est_N1_10 - est_N2_10
difer_15=est_N1_15 - est_N2_15
difer_20=est_N1_20 - est_N2_20
difer_30=est_N1_30 - est_N2_30
difer_50=est_N1_50 - est_N2_50
difer_60=est_N1_60 - est_N2_60
difer_100=est_N1_100 - est_N2_100
difer_200=est_N1_200 - est_N2_200
difer_500=est_N1_500 - est_N2_500
res_difer=data.frame(difer_5,difer_10,difer_15,difer_20,difer_30,difer_50,difer_60,difer_100,difer_200,difer_500)
sapply(res_difer, mean)## difer_5 difer_10 difer_15 difer_20 difer_30
## -0.0156000000 -0.0118000000 -0.0005333333 -0.0081000000 0.0033333333
## difer_50 difer_60 difer_100 difer_200 difer_500
## 0.0005200000 -0.0040333333 0.0035800000 -0.0024900000 -0.0001840000sapply(res_difer, sd)## difer_5 difer_10 difer_15 difer_20 difer_30 difer_50 difer_60
## 0.18127980 0.13856103 0.10685554 0.09411219 0.07651098 0.05539917 0.05271370
## difer_100 difer_200 difer_500
## 0.03961467 0.02752936 0.01424092sapply(res_difer,shapiro.test)## difer_5 difer_10
## statistic 0.8938953 0.9534241
## p.value 4.213882e-18 1.864966e-11
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## difer_15 difer_20
## statistic 0.9626701 0.9716328
## p.value 5.916738e-10 2.938903e-08
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## difer_30 difer_50
## statistic 0.9795659 0.9861968
## p.value 1.793712e-06 0.0001117251
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## difer_60 difer_100
## statistic 0.9889242 0.9933508
## p.value 0.0007841781 0.02648674
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## difer_200 difer_500
## statistic 0.9914292 0.9953602
## p.value 0.005457943 0.1422684
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"qqnorm(difer_5,main='Grafico QQ de Normalidad estimador 5 P1-P2',col="lightcoral")
qqline(difer_5)
qqnorm(difer_50,main='Grafico QQ de Normalidad estimador 50 P1-P2',col="lightcoral")
qqline(difer_50)
qqnorm(difer_100,main='Grafico QQ de Normalidad estimador 100 P1-P2',col="lightcoral")
qqline(difer_100)
qqnorm(difer_500,main='Grafico QQ de Normalidad estimador 500 P1-P2',col="lightcoral")
qqline(difer_500)
Respuesta. Si existe una relación evidente cuando se realiza el aumento de los estimadores y la normalidad que se presenta en la diferencia de P1-P2, por ejemplo, entre más aumenta el estimador en este caso al contar con la misma población de individuos enfermos el valor de la media obtenido se acerca cada vez mas a 0 sin embargo no llega a ser ese valor, en términos de variabilidad esta también disminuye pasando aproximadamente del 20% con un estimador de 5 muestras a menos del 2% con un estimador de 500 muestras, esto también se refleja en las pruebas de normalidad en donde solo serian aceptados los estimadores de 200 y 500 con un valor p sobre el 5%.
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 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%)?.
PROPORCION DE ENFERMOS DIFERENTES (P1=0.1 y P2=0.15) CON MISMO TAMANO N
lotePlantasN1=c(rep("sanas",900),rep("enfermas",100))
lotePlantasN2=c(rep("sanas",850),rep("enfermas",150))
plan_aleat_N1=function(i){
muestra=sample(lotePlantasN1, size=i)
return(sum(muestra=="enfermas")/i)
}
plan_aleat_N2=function(j){
muestra=sample(lotePlantasN2, size=j)
return(sum(muestra=="enfermas")/j)
}
i=j=72
plan_aleat_N1(i)## [1] 0.04166667plan_aleat_N2(j)## [1] 0.1666667estimadoresN1=sapply(rep(i,200), plan_aleat_N1)
estimadoresN2=sapply(rep(j,200), plan_aleat_N2)
resu_segundoP=data.frame(estimadoresN1,estimadoresN2)
sapply(resu_segundoP, mean)## estimadoresN1 estimadoresN2
## 0.09694444 0.14701389sapply(resu_segundoP, sd)## estimadoresN1 estimadoresN2
## 0.03430331 0.04149830g21=ggplot(resu_segundoP,aes(x=estimadoresN1))+geom_histogram(fill = "lightgreen",col = 1)+geom_density(col = 2,lwd=0.8)+xlab(paste("Muestra P1-0.1 n=",i))+theme_bw()
g22=ggplot(resu_segundoP,aes(x=estimadoresN2))+geom_histogram(fill = "lightgreen",col = 1)+geom_density(col = 2,lwd=0.8)+xlab(paste("Muestra P2-0.15 n=",j))+theme_bw()
ggarrange(g21,g22,ncol = 2, nrow = 1)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
estimadoresN1=sapply(rep(500,500), plan_aleat_N1)
estimadoresN2=sapply(rep(500,500), plan_aleat_N2)
difer=estimadoresN1-estimadoresN2
resu_segundoP=data.frame(estimadoresN1,estimadoresN2,difer)
sapply(resu_segundoP, mean)## estimadoresN1 estimadoresN2 difer
## 0.100068 0.150868 -0.050800sapply(resu_segundoP, sd)## estimadoresN1 estimadoresN2 difer
## 0.009616461 0.011177856 0.014317891g23=ggplot(resu_segundoP,aes(x=estimadoresN1))+geom_histogram(fill = "lightgreen",col = 1)+geom_density(col = 2,lwd=0.8)+xlab("Muestra n=500 P1")+theme_bw()
g24=ggplot(resu_segundoP,aes(x=estimadoresN2))+geom_histogram(fill = "lightgreen",col = 1)+geom_density(col = 2,lwd=0.8)+xlab("Muestra n=500 P2")+theme_bw()
g25=ggplot(resu_segundoP,aes(x=difer))+geom_histogram(fill = "lightgreen",col = 1)+geom_density(col = 2,lwd=0.8)+xlab("Muestra n=500 P1-P2")+theme_bw()
ggarrange(g23,g24,g25,ncol = 3, nrow = 1)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
a=mean(difer)
b=sd(difer)
print(paste("La media de P1-P2 es= ",format(a, scientific = F)))## [1] "La media de P1-P2 es= -0.0508"print(paste("La desviacion estandar de P1-P2 es= ",format(b, scientific = F)))## [1] "La desviacion estandar de P1-P2 es= 0.01431789"est_N1_5=sapply(rep(5,500), plan_aleat_N1)
est_N1_10=sapply(rep(10,500), plan_aleat_N1)
est_N1_15=sapply(rep(15,500), plan_aleat_N1)
est_N1_20=sapply(rep(20,500), plan_aleat_N1)
est_N1_30=sapply(rep(30,500), plan_aleat_N1)
est_N1_50=sapply(rep(50,500), plan_aleat_N1)
est_N1_60=sapply(rep(60,500), plan_aleat_N1)
est_N1_100=sapply(rep(100,500), plan_aleat_N1)
est_N1_200=sapply(rep(200,500), plan_aleat_N1)
est_N1_500=sapply(rep(500,500), plan_aleat_N1)
est_N2_5=sapply(rep(5,500), plan_aleat_N2)
est_N2_10=sapply(rep(10,500), plan_aleat_N2)
est_N2_15=sapply(rep(15,500), plan_aleat_N2)
est_N2_20=sapply(rep(20,500), plan_aleat_N2)
est_N2_30=sapply(rep(30,500), plan_aleat_N2)
est_N2_50=sapply(rep(50,500), plan_aleat_N2)
est_N2_60=sapply(rep(60,500), plan_aleat_N2)
est_N2_100=sapply(rep(100,500), plan_aleat_N2)
est_N2_200=sapply(rep(200,500), plan_aleat_N2)
est_N2_500=sapply(rep(500,500), plan_aleat_N2)
difer_5=est_N1_5 - est_N2_5
difer_10=est_N1_10 - est_N2_10
difer_15=est_N1_15 - est_N2_15
difer_20=est_N1_20 - est_N2_20
difer_30=est_N1_30 - est_N2_30
difer_50=est_N1_50 - est_N2_50
difer_60=est_N1_60 - est_N2_60
difer_100=est_N1_100 - est_N2_100
difer_200=est_N1_200 - est_N2_200
difer_500=est_N1_500 - est_N2_500
res_difer=data.frame(difer_5,difer_10,difer_15,difer_20,difer_30,difer_50,difer_60,difer_100,difer_200,difer_500)
sapply(res_difer, mean)## difer_5 difer_10 difer_15 difer_20 difer_30 difer_50
## -0.04680000 -0.04960000 -0.05173333 -0.05540000 -0.05273333 -0.05436000
## difer_60 difer_100 difer_200 difer_500
## -0.05286667 -0.04780000 -0.05220000 -0.05144000sapply(res_difer, sd)## difer_5 difer_10 difer_15 difer_20 difer_30 difer_50 difer_60
## 0.21142175 0.14373991 0.11699922 0.09960262 0.08245885 0.06387135 0.06056407
## difer_100 difer_200 difer_500
## 0.04517124 0.02963641 0.01475105sapply(res_difer,shapiro.test)## difer_5 difer_10
## statistic 0.9228634 0.9572425
## p.value 2.495263e-15 7.333358e-11
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## difer_15 difer_20
## statistic 0.9698454 0.9756492
## p.value 1.27834e-08 2.146419e-07
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## difer_30 difer_50
## statistic 0.9836708 0.9902963
## p.value 2.111883e-05 0.002228866
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## difer_60 difer_100
## statistic 0.9890993 0.992207
## p.value 0.000893785 0.01025991
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"
## difer_200 difer_500
## statistic 0.9948089 0.9968206
## p.value 0.09011493 0.435835
## method "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
## data.name "X[[i]]" "X[[i]]"qqnorm(difer_5,main='Grafico QQ de Normalidad estimador 5 P1-P2',col="lightcoral")
qqline(difer_5)
qqnorm(difer_50,main='Grafico QQ de Normalidad estimador 50 P1-P2',col="lightcoral")
qqline(difer_50)
qqnorm(difer_100,main='Grafico QQ de Normalidad estimador 100 P1-P2',col="lightcoral")
qqline(difer_100)
qqnorm(difer_500,main='Grafico QQ de Normalidad estimador 500 P1-P2',col="lightcoral")
qqline(difer_500)
Conclusion. Finalmente contrastando los dos escenarios aunque sucede el mismo comportamiento del porcentaje que se acerca cada vez al valor real de la diferencia y la variabilidad que cada vez es menor cuando aumentan los estimadores, si se observa una diferencia considerable a la media de los datos, para este segundo escenario los datos no se acumulan entorno al valor 0, si no por el contrario se acercan al valor -0.05 que representa un -5%, lo cual tiene lógica debido a que es la misma población, pero el tratamiento en N2 tuvo un resultado 5% mas de enfermos sobre N1, precisamente la diferencia que muestra la media.