Inferencia estadística y simulación

1.El Teorema del Límite Central es uno de los más importantes en la inferencia estadística y habla sobre la convergencia de los estimadores como la proporción muestral a la distribución normal. Algunos autores afirman que esta aproximación es bastante buena a partir del umbral n>30.

  1. Realice una simulación en la cual genere una población de N=1000 (Lote) y además que el porcentaje de individuos (plantas) enfermas sea del 50%.
lote=c(rep("Enfermo",500),rep("Sana",500))
  1. Genere una función que permita obtener una muestra aleatoria de la población y calcule el estimador de la proporción muestral para un tamaño de muestra dado n.
calc_p_gorro = function(n){
muestra = sample(lote, size = n)
p_gorro = sum(muestra == "Enfermo")/n
return(p_gorro)
}

calc_p_gorro(n=100)
## [1] 0.44
  1. 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?
posibles_p_gorro = sapply(rep(100, 500), calc_p_gorro)
summary(posibles_p_gorro)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3500  0.4700  0.5000  0.5016  0.5300  0.6500
hist(posibles_p_gorro)
line = mean(posibles_p_gorro)
abline(v=0.50, col="red", lwd=4)

De acuerdo con los indicadores calculados mediante la función Summary, vemos como la media y la mediana tienen un valor cercano, que permite inferir que los datos analizados no presentan tanta dispersión, pero la distribución noe s completamente simética debido a la que se puede observar en el histograma posibles_p_gorro.

  1. 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).
n = 5
require(car)
## Loading required package: car
## Warning: package 'car' was built under R version 4.1.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.1.3
posibles_p_gorro_5 = sapply(rep(5, 500), calc_p_gorro)
par(mfrow=c(1,3))
hist(posibles_p_gorro_5, las=1, ylab = "Frecuencia", main = "", col = "blue")
plot(density(posibles_p_gorro_5), las=1, ylab = "Densidad", main = "", col ="blue")
qqPlot(posibles_p_gorro_5, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales",las=1,main="")

## [1]  7 20
shapiro.test(posibles_p_gorro_5)
## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro_5
## W = 0.92621, p-value = 5.798e-15
n = 10
require(car)
posibles_p_gorro_10 = sapply(rep(10, 500), calc_p_gorro)
par(mfrow=c(1,3))
hist(posibles_p_gorro_10, las=1, ylab = "Frecuencia", main = "", col = "blue")
plot(density(posibles_p_gorro_10), las=1, ylab = "Densidad", main = "", col ="blue")
qqPlot(posibles_p_gorro_10, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales",las=1,main="")

## [1]  81 153
shapiro.test(posibles_p_gorro_10)
## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro_10
## W = 0.9607, p-value = 2.711e-10
n = 15
require(car)
posibles_p_gorro_15 = sapply(rep(15, 500), calc_p_gorro)
par(mfrow=c(1,3))
hist(posibles_p_gorro_15, las=1, ylab = "Frecuencia", main = "", col = "blue")
plot(density(posibles_p_gorro_15), las=1, ylab = "Densidad", main = "", col ="blue")
qqPlot(posibles_p_gorro_15, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales",las=1,main="")

## [1]  99 234
shapiro.test(posibles_p_gorro_15)
## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro_15
## W = 0.97074, p-value = 1.932e-08
n = 20
require(car)
posibles_p_gorro_20 = sapply(rep(20, 500), calc_p_gorro)
par(mfrow=c(1,3))
hist(posibles_p_gorro_20, las=1, ylab = "Frecuencia", main = "", col = "blue")
plot(density(posibles_p_gorro_20), las=1, ylab = "Densidad", main = "", col ="blue")
qqPlot(posibles_p_gorro_20, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales",las=1,main="")

## [1] 479  76
shapiro.test(posibles_p_gorro_20)
## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro_20
## W = 0.98103, p-value = 4.202e-06
n = 30
require(car)
posibles_p_gorro_30 = sapply(rep(30, 500), calc_p_gorro)
par(mfrow=c(1,3))
hist(posibles_p_gorro_30, las=1, ylab = "Frecuencia", main = "", col = "blue")
plot(density(posibles_p_gorro_30), las=1, ylab = "Densidad", main = "", col ="blue")
qqPlot(posibles_p_gorro_30, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales",las=1,main="")

## [1] 121 308
shapiro.test(posibles_p_gorro_30)
## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro_30
## W = 0.97999, p-value = 2.292e-06
n = 50
require(car)
posibles_p_gorro_50 = sapply(rep(50, 500), calc_p_gorro)
par(mfrow=c(1,3))
hist(posibles_p_gorro_50, las=1, ylab = "Frecuencia", main = "", col = "blue")
plot(density(posibles_p_gorro_50), las=1, ylab = "Densidad", main = "", col ="blue")
qqPlot(posibles_p_gorro_50, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales",las=1,main="")

## [1] 155 219
shapiro.test(posibles_p_gorro_50)
## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro_50
## W = 0.99033, p-value = 0.002289
n = 60
require(car)
posibles_p_gorro_60 = sapply(rep(60, 500), calc_p_gorro)
par(mfrow=c(1,3))
hist(posibles_p_gorro_60, las=1, ylab = "Frecuencia", main = "", col = "blue")
plot(density(posibles_p_gorro_60), las=1, ylab = "Densidad", main = "", col ="blue")
qqPlot(posibles_p_gorro_60, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales",las=1,main="")

## [1]  24 378
shapiro.test(posibles_p_gorro_60)
## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro_60
## W = 0.99187, p-value = 0.007812
n = 100
require(car)
posibles_p_gorro_100 = sapply(rep(100, 500), calc_p_gorro)
par(mfrow=c(1,3))
hist(posibles_p_gorro_100, las=1, ylab = "Frecuencia", main = "", col = "blue")
plot(density(posibles_p_gorro_100), las=1, ylab = "Densidad", main = "", col ="blue")
qqPlot(posibles_p_gorro_100, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales",las=1,main="")

## [1] 179  12
shapiro.test(posibles_p_gorro_100)
## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro_100
## W = 0.99493, p-value = 0.09988
n = 200
require(car)
posibles_p_gorro_200 = sapply(rep(200, 500), calc_p_gorro)
par(mfrow=c(1,3))
hist(posibles_p_gorro_200, las=1, ylab = "Frecuencia", main = "", col = "blue")
plot(density(posibles_p_gorro_200), las=1, ylab = "Densidad", main = "", col ="blue")
qqPlot(posibles_p_gorro_200, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales",las=1,main="")

## [1] 460 317
shapiro.test(posibles_p_gorro_200)
## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro_200
## W = 0.99573, p-value = 0.192
n = 500
require(car)
posibles_p_gorro_500 = sapply(rep(500, 500), calc_p_gorro)
par(mfrow=c(1,3))
hist(posibles_p_gorro_500, las=1, ylab = "Frecuencia", main = "", col = "blue")
plot(density(posibles_p_gorro_500), las=1, ylab = "Densidad", main = "", col ="blue")
qqPlot(posibles_p_gorro_500, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales",las=1,main="")

## [1] 202 214
shapiro.test(posibles_p_gorro_500)
## 
##  Shapiro-Wilk normality test
## 
## data:  posibles_p_gorro_500
## W = 0.99439, p-value = 0.06326
Análisis

En la mayoria de las muestras aleatorias realizadas, teniendo en cuenta la forma de las gráficos, se puede inferir que se presenta una distribución normal y la mayoria de los puntos de las muestras se encuentra dentro de los intervalos de confianza.

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.

  1. 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 elcontrol de una plaga presente en ambos al momento inicial. Para ello utilizará como criterio de desempeño el tratamiento que menor % de plantas enfermaspresente 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ónen 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).
lote1=c(rep("Enfermo",100),rep("Sana",900))
lote2=c(rep("Enfermo",150),rep("Sana",1350))
  1. 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_diferencia=function(m){

lote1sample=sample(lote1,m)
p1=sum(lote1sample=="Enfermo")/m

lote2sample=sample(lote2,m)
p2=sum(lote2sample=="Enfermo")/m

diferenciap=p1-p2
return(diferenciap)
}
calc_diferencia(m=250)
## [1] 0
  1. 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?
simula_lote=sapply(rep(250,500),calc_diferencia)
summary(simula_lote)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -0.080000 -0.016000  0.000000  0.001864  0.020000  0.076000
hist(simula_lote)
line=mean(simula_lote)
abline(v=0.0, col="red", lwd=4)

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

De acuerdo con la lectura, el valor p fue considerado en su momento el estandar como criterio de decisión sobre si la inferencia a cerca de un conjunto de datos era correcta, y fue un factor determinante en sus casi 9 decadas desde su implementación. Sin embargo, con el paso del tiempo y nuevos estudios relacionados con la estadísitica, se ha identificado que existen otros factores que tambien se deben tener en cuenta para garantizar un correcto análisis de los datos y no basarse unicamente en un solo calculo que puede conllevar a conclusiones erroneos que difieran de la realidad de la situación analizada.

El objetivo del valor p es que permite un indicador apriori sobre los resultados que son productos del azar o de aquellos resultados que son estadisticamente significativos, mediante un único calculo. Lo anaterio conlleva a que la relevancia pese sobre un único argumento.El autor plantea que el valor p puede ser un indcador de apoyo para la comprensión del conjunto de datos, pero que en adición, es necesario utilizar criterios adicionales de soporte como por ejemplo los intervalos de confianza o métodos bayesianos, entre otros. La correlación de estas difrentes medidas e indicadores permiten tener un mejor sustento para garantizar

Para concluir, al momento de realizar un análisis exploratorio de los datos es necesario identificar el contexto en el cual se generan los datos a analizar y contar con el apoyo de un experto de la situación que se esta analizando, el cual nos puede brindar orientaciones sobre si los indicadores o medidas que se realizan pueden ser posibles, y en adición utilizar diferentes indicadores estadísticos que permitan una mejor comrpensión de los datos y de la situación de estudio.