Taller de simulación en R
German Avila Rodriguez
2022-09-2022
Actividad 1 - Datos Rotación
El Teorema del Límite Central es uno de los más importantes en la inferencia estadística y habla sobre la convergencia de los estimadores como la proporción muestral a la distribución normal. Algunos autores afirman que esta aproximación es bastante buena a partir del umbral n>30.
- Realice una simulación en la cual genere una población de N=1000 (Lote) y además que el porcentaje de individuos (plantas) enfermas sea del 50%. 1000 ceros o unos generados aleatoriamente con la misma probabilidad, el 1 representará a las plantas enfermas
set.seed(1)
Lote=sample(0:1,1000, replace=T)- 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.
MuestraLote=function(n){
muestra=sample(Lote,n, replace=T)
return(sum(muestra)/n)
} Muestra de 50% plantas del lote enfermas
sum(MuestraLote(50))/50## [1] 0.0116El estimador de la proporción para un tamaño de la muestra de 50 es de 0.58
Muestra del estimador de la proporción de un lote de 50 plantas repetido 500 veces
C. Repita el escenario anterior (b) 500 veces y analice los resultados en cuanto al comportamiento de los 500 estimadores. ¿Qué tan simétricos son los datos?, ¿Son sesgados y qué pasa en cuanto a variabilidad?
porcentajes_muestra=sapply(rep(50,500), MuestraLote)
hist(porcentajes_muestra,freq=FALSE,col="lightsalmon",main="Histograma",sub="Datos simulados de una N(50,500)")
abline(v=0.5,col="red",lwd=3)
Como se puede apreciar en la figura anterior, la media calculada para la muestra corresponde a la limitación que se tiene para la población, es decir, 50% de plantas enfermas.
Muestra de la prueba de tamaños de muestras n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500
- 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).
res=data.frame(p_m5,p_m10,p_m15,p_m30,p_m50,p_m60,p_m100,p_m200,p_m500)
boxplot(res, col = "lightsalmon",main="Boxplot de Diferentes muestras",outbg = "green",xlab = "Muestras",outbg = "blue")
abline(h=0.5,col="red",lwd=3)
Los diferentes tamaños de muestras tienden a oscilar a un valor de 0.5 pero entre más grande la muestra el margen de error es mucho menor
##
## Shapiro-Wilk normality test
##
## data: porcentajes_muestra
## W = 0.92574, p-value = 5.135e-15## [1] 0.4868
##
## Shapiro-Wilk normality test
##
## data: porcentajes_muestra
## W = 0.96074, p-value = 2.759e-10## [1] 0.5094
##
## Shapiro-Wilk normality test
##
## data: porcentajes_muestra
## W = 0.97214, p-value = 3.748e-08## [1] 0.4837333
##
## Shapiro-Wilk normality test
##
## data: porcentajes_muestra
## W = 0.98608, p-value = 0.000103## [1] 0.4991333
##
## Shapiro-Wilk normality test
##
## data: porcentajes_muestra
## W = 0.98926, p-value = 0.001007## [1] 0.49828
##
## Shapiro-Wilk normality test
##
## data: porcentajes_muestra
## W = 0.98917, p-value = 0.0009435## [1] 0.4971333
##
## Shapiro-Wilk normality test
##
## data: porcentajes_muestra
## W = 0.9944, p-value = 0.06404## [1] 0.49618
##
## Shapiro-Wilk normality test
##
## data: porcentajes_muestra
## W = 0.99531, p-value = 0.1368## [1] 0.49688
##
## Shapiro-Wilk normality test
##
## data: porcentajes_muestra
## W = 0.99525, p-value = 0.1298## [1] 0.498908Muestra N |
Estimador de la media |
P.valor |
|---|---|---|
| 5 | 0.48 | 5.135e-15 |
| 10 | 0.50 | 2.759e-10 |
| 15 | 0.48 | 3.748e-08 |
| 30 | 0.49 | 0.000103 |
| 50 | 0.49 | 0.001007 |
| 60 | 0.49 | 0.0009435 |
| 100 | 0.49 | 0.06404 |
| 200 | 0.49 | 0.1368 |
| 500 | 0.49 | 0.1298 |
De acuerdo a las muestras pequeñas en el gráfico QQ se observa un mayor desajuste con la linea recta en los extremos, conservando una buena meidad de tendencia centra que a medida que el valor de n crece se obtiene con mejor precisión una linea recta, y también, se presenta un incremento del valor p de la prueba Shapiro-Wilk en la que se va aceptando la hipotesis nula de normalidad.
Con un número de muestras de 60 sobre la población, se observa una gráfica mucho más continua y uniforme al valor de la media de la población. El valor de P en este punto ha mejorado ampliamente comparándolo con la primera muestra de 5.
Muestra de 10% plantas del lote enfermas
Se presenta una muestra con el 10% de plantas enfermas, evidenciando una muestra en el Histograma y con el grafico de cuantiles
MuestraLote_90=function(n){
muestra=sample(Lote_90,n, replace=T)
return(sum(muestra)/n)
} 
##
## Shapiro-Wilk normality test
##
## data: porcentajes_muestra
## W = 0.98647, p-value = 0.0001347## [1] 0.10196Muestra N |
Estimador de Media |
P. Valor |
|---|---|---|
| 5 | 0.0844 | < 2.2e-16 |
| 10 | 0.0934 | < 2.2e-16 |
| 15 | 0.0952 | < 2.2e-16 |
| 30 | 0.1022 | 0,000000001917 |
| 50 | 0.10324 | 0,000003579 |
| 60 | 0.1004667 | 0,001895 |
| 100 | 0.10188 | 0.0002779 |
| 200 | 0.10014 | 0.01764 |
| 500 | 0.099584 | 0.3138 |
#ciclo para tamaño de muestra 10%
par(mfrow=c(1,6))
n<-c(5, 10, 15, 30, 50, 60, 100, 200, 500)
for (i in 1:9) {
porcentajes_muestra=sapply(rep(n[i],500), MuestraLote_10)
hist(porcentajes_muestra,freq=FALSE,col="lightsalmon",main = paste("Histograma"))
abline(v=0.1,col="red",lwd=3)
qqnorm(porcentajes_muestra, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales",las=1,main=paste("Gráfico QQ"))
qqline(porcentajes_muestra, col = "red")
# print(shapiro.test(porcentajes_muestra))
# print(mean(porcentajes_muestra))
}
Muestra de 90% plantas del lote enfermas
Se presenta una muestra con el 90% de plantas enfermas, evidenciando una muestra en el Histograma y con el grafico de cuantiles
##
## Shapiro-Wilk normality test
##
## data: porcentajes_muestra
## W = 0.98425, p-value = 3.063e-05## [1] 0.90076Muestra N |
Estimador de Media |
P. Valor |
|---|---|---|
| 5 | 0.902 | < 2.2e-16 |
| 10 | 0.899 | < 2.2e-16 |
| 15 | 0.8998667 | < 2.2e-16 |
| 30 | 0.9034 | 0,00000001171 |
| 50 | 0.90156 | 0,000005005 |
| 60 | 0.9012667 | 0,0000005141 |
| 100 | 0.90074 | 0.0005241 |
| 200 | 0.90095 | 0.01766 |
| 500 | 0.899668 | 0.05123 |
#ciclo para tamaño de muestra 90%
Se desarrolla un ciclo el cual permite evidenciar el comportamiento de la muestra en un rango de 500 en donde la tendencia es que entre menor la muestra mayor el P.valor sosteniendo se el estimador sobre el 90%
par(mfrow=c(1,6))
n<-c(5, 10, 15, 30, 50, 60, 100, 200, 500)
for (i in 1:9) {
porcentajes_muestra=sapply(rep(n[i],500), MuestraLote_90)
hist(porcentajes_muestra,freq=FALSE,col="lightsalmon",main = paste("Histograma"))
abline(v=0.1,col="red",lwd=3)
qqnorm(porcentajes_muestra, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales",las=1,main=paste("Gráfico QQ"))
qqline(porcentajes_muestra, col = "red")
# print(shapiro.test(porcentajes_muestra))
# print(mean(porcentajes_muestra))
}
## Conclusión
Los valores muestras pequeños de no arrojan simetría alrededor de la media de la población. Mientras las muestras sean más pequeñas el estimador tendrá menos valores continuos, como se presenta en las graficas de igual manera pasa con el valor de P, en la medida que las muestras sean más pequeñas este valor se aleja mucho más del valor esperado, de las simulaciones se presentan valores de hasta con numeros exponenciales lo que definitivamente rechazaría cualquier hipótesis.
De igual manera que en la medida que aumentan las muestras el P.Valor se acerca más a un valor óptimo que permitiría aceptar la hipótesis planteada.
Como se evidencia en los diferentes graficos, los resultados con muestras más pequeñas se dispersan alrededor de la línea teórica de la normal; mientras que cuando las muestras son mayores estos valores se concentran sobre la línea teórica de la normal.
Recomendación de acuerdo a los valores de muestra muy pequeños el resultado del p-Value rechazaría cualquier hipótesis, por lo que es pertinente estimar muy bien los tamaños de las muestras y evitar resultados sesgados.
Actividad 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). Lotes de población 1000 y 1500 respectivamente.
LoteC1000 = sample(0:1, 1000, replace=TRUE, prob = c(90,10))
LoteC1500 = sample(0:1, 1500, replace=TRUE, prob = c(90,10))- 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.
MuestraLoteCar=function(n,Lote){
return(sample(Lote,n, replace=T))
}
p1=sum(MuestraLoteCar(40,LoteC1000))/40
p2=sum(MuestraLoteCar(40,LoteC1500))/40
p1-p2## [1] -0.125- 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?
simulNtimesCp<-vector("numeric",500)
for (i in 1:500) {
p1=sum(MuestraLoteCar(40,LoteC1000))/40
p2=sum(MuestraLoteCar(40,LoteC1500))/40
simulNtimesCp[i]<-(p1-p2)
}
simulNtimesCp <- as.numeric(simulNtimesCp) # turn into numbers
hist(simulNtimesCp,density = 50, col="lightsalmon", main = "Histograma de diferencias p1-p2",
xlab = "Estimador", ylab = "Frecuencia")
line = mean(simulNtimesCp)
abline(v=line, col="red",lwd=3)
Se presenta que hay simetría en los estimadores. No hay un sesgo de la media igual 0. tomando las poblaciones de diferente tamaño las muestras que se extraen de ellas de igual tamaño no varian respecto a la proporción de plantas enfermas.
- Realice los puntos b y c para tamaños de muestra n1=n2=5, 10, 15, 20, 30, 50, 60, 100, 200, 500. Y compare los resultados de los estimadores (p1-p2) en cuanto a la normalidad. También analice el comportamiento de las diferencias y evalúe. ¿Considera que es más probable concluir que existen diferencias entre los tratamientos con muestras grandes que pequeñas, es decir, cuál considera usted que es el efecto del tamaño de muestra en el caso de la comparación de proporciones?
par(mfrow=c(1,6))
funcionSimulCp=function(n){
SimulNtimesCp<-vector("numeric",500)
for (i in 1:500) {
p1=sum(MuestraLoteCar(n,LoteC1000))/n
p2=sum(MuestraLoteCar(n,LoteC1500))/n
simulNtimesCp[i]<-p1-p2
}
simulNtimesCp <- as.numeric(simulNtimesCp) # turn into numbers
hist(simulNtimesCp,density = 50, col="lightsalmon", main = paste("Histograma p1-p2 para n=", n),
xlab = "Estimador", ylab = "Frecuencia")
line = mean(simulNtimesCp)
abline(v=line, col="red",lwd=3)
qqnorm(simulNtimesCp, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales",las=1,main=paste("Gráfico QQ n=",n))
qqline(simulNtimesCp, col = "red")
# print(paste("Resumen Estadístico n=",n))
# print(summary(simulNtimesCp))
# print(shapiro.test(simulNtimesCp))
}
vecMuestraCp <- c(5,10,15,20,30,50,60,100,200,500)
for (i in vecMuestraCp) {
funcionSimulCp(i)
}
Se identifica que a medida que incrementa el tamaño de la muestra la diferencia de la proporción tiene una leve tendencia hacía la izquierda, es posible pensar que si se hace el análisis solo con un tamaño de la muestra, el efecto es causa de la aleatoriedad, de igual manera los resultados de la diferencia presentarían una forma de distribución normal, centrada con media del parámetro en cero (0).
A medida que aumentamos en tamaño de la muestra (n=200, 500) se evidencia que la diferencia entre estas poblaciones es pequeña ~0.025 esto se puede interpretar que las poblaciones tienen mucha probabilidad de ser iguales.
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%)?
Lote2C1000 = sample(0:1, 1000, replace=TRUE, prob = c(90,10))
Lote2C1500 = sample(0:1, 1500, replace=TRUE, prob = c(85,15))par(mfrow=c(1,6))
funcionSimul2Cp=function(n){
simulNtimes2Cp<-vector("numeric",500)
for (i in 1:500) {
p1=sum(MuestraLoteCar(n,Lote2C1000))/n
p2=sum(MuestraLoteCar(n,Lote2C1500))/n
simulNtimes2Cp[i]<-p1-p2
}
simulNtimes2Cp <- as.numeric(simulNtimes2Cp) # turn into numbers
hist(simulNtimes2Cp,density = 50, main = paste("Histograma p1-p2 para n=", n), col="lightsalmon",
xlab = "Estimador", ylab = "Frecuencia")
line = mean(simulNtimes2Cp)
abline(v=line, col="red",lwd=3)
qqnorm(simulNtimes2Cp, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales",las=1,main=paste("Gráfico QQ n=",n))
qqline(simulNtimes2Cp, col = "red")
# print(paste("Resumen Estadístico n=",n))
# print(summary(simulNtimes2Cp))
# print(shapiro.test(simulNtimes2Cp))
}
vecMuestra2Cp <- c(5,10,15,20,30,50,60,100,200,500)
for (i in vecMuestra2Cp) {
funcionSimul2Cp(i)
}
Se concluye que hay un sesgo notable con respecto a la media 0, pero se conserva la simetría alrededor de la media muestral.
En el primer escenario, el mayor sesgo respecto a la media fue para el tamaño de muestra igual a 500, aún así fue muy leve. En el segundo escenario se observa un valor de media aproximado a -0.06 en todos los casos, consecuencia de tomar un valor de proporción menor de plantas enfermas en un lote que en otro. Variar el tamaño de la población y el valor de la proporción de plantas enfermas hace que poder concluir hacerca del factor que hace que la media permanezca más estable pueda tornarse ambiguo, pero considerando que en el escenario uno la tendencia de la media de la diferencia de las proporciones era levemente negativa, podríamos sospechar que en el escenario 2 esta tendencia de un valor de media estable se deba a que se tomaron proporciones diferentes de plantas enfermas.
Actividad 3.
Con base a los artículos “Statistical Errors: P values, the gold standard of statistical validity, are not as reliable as many scientists assume” & “Statisticians issue warning on P values: Statement aims to halt missteps in the quest for certainty” escriba un resumen (máximo 2 páginas) sobre ambos artículos e incluya en este sus opiniones en cuanto al uso del valor p como criterio de decisión en inferencia estadística.
Con base en el artículo iniciamos de acuerod a la especificación de cuáles son los valores esperados o estimados para un óptimo resultado de los cálculos del p-value; es de relevancia para los resultados estadísticamente significantes que los valores de p se encuentren entre 0.01 y 0.05, mientras más pequeño sea este valor las probabilidades de rechazar la hipótesis nula.
también cabe aclarar que el entorno en el cual se da el uso de p-value por parte del estadístico Ronald Fisher se indica que dicho resultado no representa una prueba definitiva, un camino para validar la significancia de los resultados, de acuerdo a esto es menester que los resultados no solamente se deben basar en los datos obtenidos del valor de p, de acuerdo a que que podrían inducir a afirmaciones falsas. lo que nos llevaría a tomar los datos basados en una hipótesis especificada.
Por ende el solo valor de p no es un valor absoluto, este debe ser complementado con más evidencias para soportar una hipótesis, y es necesario demostrar que ese valor de p conlleva una probabilidad de falsos hallazgos, que a medida que el valor de p disminuye esa probabilidad de falsa alarma aumenta.
teniendo en cuenta al autor Geoff Cummin interpreta que la significancia no es un indicador práctico; este debería permitir cuestionarnos sobre ¿Cuándo efecto existe?, se plantea complementar los resultados con el tamaño de los efectos y los intervalos de confianza y alejarse un poco de la simple premisa si los resultados son o no significativos.
en conclusión el análisis exploratorio de datos es necesario orientarlo en diferentes etapas, realizando pequeños análisis exploratorios de los datos con la objetivo de generar hallazgos relevantes sin la necesidad de preocuparnos sobre las falsas hallazgos [1].
[1] R. Nuzzo. “Statistical errors: P values, the ‘gold standard’ of statistical validity, are not as reliable as many scientists assume.” Nature, vol. 506, no. 7487, pp. 150+, 13 Feb. 2014. Consultado Mar. 02, 2022. Disponible en: link.gale.com/apps/doc/A362064378/HRCA?u=anon~69369c4d&sid=googleScholar&xid=33e82427.