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

Literal A

Realice una simulación en la cual genere una población de N=1000 (Lote) y además que el porcentaje de individuos (plantas) enfermas sea del 50%.

Solución

Se genera la simulación empleando el siguiente código:

simulacion=sample(c(rep("Enferma",500),rep("Sana",500)))

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

Solución

Se genera la función empleando el siguiente código:

funcion=function(n){
  b=sample(simulacion,size=n)
  p_circunflejo=sum(b=="Enferma")/n
return(p_circunflejo)
}
funcion(8)
## [1] 0.625

Literal 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?

Solución

Se genera la función que toma como variable el tamaño de muestra deseado \(n\), crea un vector con \(n\) registros, y a través de un ciclo con contador \(i\) que a través de 1000 iteraciones calcula 1000 muestras de tamaño \(n\) de la variable, para luego hallar el estimador \(\hat{p}\), calculado como la suma de los valores en que la planta de la muestra está enferma, sobre el número total de observaciones \(n\), que corresponde al porcentaje de plantas enfermas en la muestra, y finalmente guarda el resultado de cada iteración en la posición i del vector anteriormente definido; el código empleado es el siguiente, y se ejecuta con un tamaño de muestra \(n\) igual a 10:

ciclo=function(n){
vector=c(1:n)
for (i in 1:1000){
  b=sample(simulacion,size=n)
  p_circunflejo=sum(b=="Enferma")/n
vector[i]=p_circunflejo
}
return(vector)}
coef_var = function(vector) {
  sd(vector) / mean(vector)
}

n=10
Ejecucion=(ciclo(n))

Ahora, se calculan diferentes indicadores de estadística descriptiva, para concluir sobre los resultados

library(psych)
promedio=mean(Ejecucion)
desviacion=sd(Ejecucion)
coeficiente_variacion=coef_var(Ejecucion)
coeficiente_asimetria=skew(Ejecucion,type=1)
kurtosis=kurtosi(Ejecucion)
hist(Ejecucion,main=paste("Histograma del valor de p circunflejo para muestras de tamaño",n),ylab="Frecuencia",xlab="Valor de p circunflejo")

Datos=data.frame(promedio,desviacion,coeficiente_variacion,coeficiente_asimetria,kurtosis)
knitr::kable(Datos, align = "c", caption = "Medidas Relevantes" )
Medidas Relevantes
promedio desviacion coeficiente_variacion coeficiente_asimetria kurtosis
0.4981 0.1588446 0.3189011 -0.0588013 -0.1572936

Para un valor de muestra \(n\) igual a 10, se logra observar un valor promedio del estimador \(\hat{p}\), bastante cercano al parámetro \(p\), generalmente entre 0.48 y 0.52, por lo cual no sería un resultado sesgado, sin embargo, el valor de la desviación estándar fue en el mayor de los casos superior a 0.15, dando a entender un rango entre 0.35 y 0.65 para el valor de \(\hat{p}\), y con un coeficiente de variación \(cv\) regularmente superior a 30%, lo que indica un resultado muy heterógeneo.

Literal 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).

Solución

Para la repetición solicitada del punto b, se decide crear un dataframe que contiene dos columnas: la primera corresponde el número de escenarios dados (1 a 10), y la segunda al tamaño de muestra \(n\) de 5 hasta un \(n\) de 500, y se adicionn dos escenarios para un totla de 12, uno con \(n= 350\), y otro con \(n= 1000\), este último correspondiente a la población completa. Por medio de un ciclo for, con un contador que va desde 1 hasta el total de escenarios (12), se ejecuta la función creada en el literal A que calcula el valor de \(\hat{p}\) para cada escenario, teniendo en cuenta que este toma la muestra de tamañao \(n\) una única vez. El código empleado se muestra a continuación:

N_Escenarios=12
Ejecuciones=data.frame(1:N_Escenarios,c(5,10,15,20,30,50,60,100,200,350,500,1000))
colnames(Ejecuciones)=c("Escenario","n")
Ejecuciones
Escenario n
1 5
2 10
3 15
4 20
5 30
6 50
7 60
8 100
9 200
10 350
11 500
12 1000 
for (i in 1:N_Escenarios){
  funcion(Ejecuciones$n[i])
  print(funcion(Ejecuciones$n[i]))
}
## [1] 0.2
## [1] 0.6
## [1] 0.5333333
## [1] 0.45
## [1] 0.4333333
## [1] 0.56
## [1] 0.45
## [1] 0.56
## [1] 0.54
## [1] 0.4885714
## [1] 0.538
## [1] 0.5

Para repetir el punto c,se genera un ciclo for que ejecuta el ciclo creado en el punto C, calculando en 1000 iteraciones el valor de \(\hat{p}\) de acuerdo al tamaño de muestra de cada escenario. Adicionalmente, se opta por crear 4 vectores con la misma dimensión del número de escenarios, cuya posición \(i\) es luego reemplazada por el ciclo for para almacenar 4 medidas de estadística descriptiva: Promedio, Desviación Estándar, Coeficiente de Variación y Coeficiente de Asimetría; estas finalmente se incluyen dentro de un dataframe con el fin de permitir la comparación posterior de los resultados de cada escenario.

#Definición de vectores para contener resultados.

Promedios=c(1:N_Escenarios)
Desviaciones=c(1:N_Escenarios)
Coeficientes_Variacion=c(1:N_Escenarios)
Coeficientes_Asimetria=c(1:N_Escenarios)

#Ejecución del ciclo for para cada escenario.

Conf2x5 = matrix(c(1:12), nrow=3, byrow=FALSE)
layout(Conf2x5)
#layout.show(12)
for (i in 1:N_Escenarios){
  Ejecucion=(ciclo(Ejecuciones$n[i]))
Promedios[i]=mean(Ejecucion)
Desviaciones[i]=sd(Ejecucion)
Coeficientes_Variacion[i]=coef_var(Ejecucion)
Coeficientes_Asimetria[i]=skew(Ejecucion,type=1)
hist(Ejecucion,main=paste("n=",Ejecuciones$n[i]),ylab="Frecuencia",xlab="Valor de p")
abline(v=Promedios[i],col="red",lwd=2)

}

Conf2x5 = matrix(c(1:12), nrow=3, byrow=FALSE)
for (i in 1:N_Escenarios){
  Ejecucion=(ciclo(Ejecuciones$n[i]))
qqnorm(Ejecucion,xlab="Cuantiles Teóricos",ylab="Cuant. Muestrales",main=paste("n=",Ejecuciones$n[i]))
qqline(Ejecucion, col="red", lwd=1)

}

Se puede apreciar en los histogramas, que para mayor valor de \(n\), la gráfica se asemeja más a una distribución normal; además, los límites inferior y superior de \(\hat{p}\) se vuelven más cercanos al parámetro \(p\) a medida que se aumenta el valor de \(n\). Además, en todas las ejecuciones, los gráficos de qq-normalidad dejan en evidencia el comportamiento de normalidad de las simulaciones. También se presenta una tabla con distintos cálculos relevantes para cada escenario:

Datos_Generales=data.frame(Ejecuciones,Promedios,Desviaciones,Coeficientes_Variacion,Coeficientes_Asimetria)
knitr::kable(Datos_Generales, align = "c", caption = "Medidas Relevantes para 10 escenarios con 1000 repeticiones de tamaño de muestra variable" )
Medidas Relevantes para 10 escenarios con 1000 repeticiones de tamaño de muestra variable
Escenario n Promedios Desviaciones Coeficientes_Variacion Coeficientes_Asimetria
1 5 0.4930000 0.2266327 0.4597011 0.0103136
2 10 0.5079000 0.1582803 0.3116367 0.1094494
3 15 0.5004000 0.1272199 0.2542365 0.0426607
4 20 0.4962000 0.1152556 0.2322764 0.0245884
5 30 0.4956000 0.0917251 0.1850789 0.0062422
6 50 0.5036800 0.0730684 0.1450692 -0.0762910
7 60 0.5015667 0.0623858 0.1243820 0.0196682
8 100 0.5002500 0.0457792 0.0915126 0.0183801
9 200 0.4995800 0.0310232 0.0620985 0.0505494
10 350 0.5000686 0.0213418 0.0426777 0.0565639
11 500 0.5002620 0.0164803 0.0329433 -0.0207981
12 1000 0.5000000 0.0000000 0.0000000 NaN

Finalmente, se trazan dos diagramas de dispersión con curva de tendencia suavizada, donde se evidencia que el tamaño de muestra \(n\) es inversamente proporcional a la desviación estándar de \(\hat{p}\), así como al coeficiente de variación, y su comportamiento es bastante similar teniendo en consideración el comportamiento previamente observado del valor promedio de \(\hat{p}\) (El sesgo es bastante bajo, es decir, se acerca bastante al parámetro \(p\)).

cor(Datos_Generales$n,Datos_Generales$Desviaciones)
## [1] -0.6654371
cor(Datos_Generales$n,Datos_Generales$Coeficientes_Variacion)
## [1] -0.6624728
require(ggplot2)
require(plotly)

Dispersion1=ggplot(Datos_Generales,aes(x=n,y=Coeficientes_Variacion)) + geom_point() + geom_smooth(method="loess")
ggplotly(Dispersion1)
Dispersion2=ggplot(Datos_Generales,aes(x=n,y=Desviaciones)) + geom_point() + geom_smooth(method="loess")
ggplotly(Dispersion2)

Literal E

Repita toda la simulación (puntos a – d) pero ahora con lotes con 10% y 90% de plantas enfermas. Concluya todo el ejercicio.

Lotes de 10% Enfermas

Generación del Lote
simulacion10=sample(c(rep("Enferma",100),rep("Sana",900)))
Función para una muestra de tamaño n
funcion10=function(n10){
  b10=sample(simulacion10,size=n10)
  p_circunflejo10=sum(b10=="Enferma")/n10
return(p_circunflejo10)
}
funcion10(12)
## [1] 0.1666667
1000 repeticiones con n = 12
ciclo10=function(n10){
vector10=c(1:n10)
for (i in 1:1000){
  b10=sample(simulacion10,size=n10)
  p_circunflejo10=sum(b10=="Enferma")/n10
vector10[i]=p_circunflejo10
}
return(vector10)}
coef_var10 = function(vector10) {
  sd(vector10) / mean(vector10)
}

n10=12
Ejecucion10=(ciclo10(n10))
library(psych)
promedio10=mean(Ejecucion10)
desviacion10=sd(Ejecucion10)
coeficiente_variacion10=coef_var10(Ejecucion10)
coeficiente_asimetria10=skew(Ejecucion10,type=1)
kurtosis10=kurtosi(Ejecucion10)
hist(Ejecucion10,main=paste("Histograma del valor de p circunflejo para muestras de tamaño",n10),ylab="Frecuencia",xlab="Valor de p circunflejo")

Datos10=data.frame(promedio10,desviacion10,coeficiente_variacion10,coeficiente_asimetria10,kurtosis10)
knitr::kable(Datos10, align = "c", caption = "Medidas Relevantes" )
Medidas Relevantes
promedio10 desviacion10 coeficiente_variacion10 coeficiente_asimetria10 kurtosis10
0.102 0.0841974 0.8254646 0.7063316 0.4205837
1000 repeticiones por cada valor de n
N_Escenarios10=12
Ejecuciones10=data.frame(1:N_Escenarios10,c(5,10,15,20,30,50,60,100,200,350,500,1000))
colnames(Ejecuciones10)=c("Escenario","n")
Ejecuciones10
Escenario n
1 5
2 10
3 15
4 20
5 30
6 50
7 60
8 100
9 200
10 350
11 500
12 1000 
for (i in 1:N_Escenarios10){
  funcion10(Ejecuciones10$n[i])
  print(funcion10(Ejecuciones10$n[i]))
}
## [1] 0
## [1] 0.1
## [1] 0
## [1] 0.2
## [1] 0.06666667
## [1] 0.12
## [1] 0.1333333
## [1] 0.09
## [1] 0.095
## [1] 0.08
## [1] 0.094
## [1] 0.1
#Definición de vectores para contener resultados.

Promedios10=c(1:N_Escenarios10)
Desviaciones10=c(1:N_Escenarios10)
Coeficientes_Variacion10=c(1:N_Escenarios10)
Coeficientes_Asimetria10=c(1:N_Escenarios10)

#Ejecución del ciclo for para cada escenario.

Conf2x5 = matrix(c(1:12), nrow=3, byrow=FALSE)
layout(Conf2x5)
#layout.show(12)
for (i in 1:N_Escenarios10){
  Ejecucion10=(ciclo10(Ejecuciones10$n[i]))
Promedios10[i]=mean(Ejecucion10)
Desviaciones10[i]=sd(Ejecucion10)
Coeficientes_Variacion10[i]=coef_var10(Ejecucion10)
Coeficientes_Asimetria10[i]=skew(Ejecucion10,type=1)
hist(Ejecucion10,main=paste("n=",Ejecuciones10$n[i]),ylab="Frecuencia",xlab="Valor de p")
}

Conf2x5 = matrix(c(1:12), nrow=3, byrow=FALSE)
for (i in 1:N_Escenarios){
  Ejecucion10=(ciclo10(Ejecuciones$n[i]))
qqnorm(Ejecucion10,xlab="Cuantiles Teóricos",ylab="Cuant. Muestrales",main=paste("n=",Ejecuciones$n[i]))
qqline(Ejecucion10, col="red", lwd=1)

}

Datos_Generales10=data.frame(Ejecuciones10,Promedios10,Desviaciones10,Coeficientes_Variacion10,Coeficientes_Asimetria10)
knitr::kable(Datos_Generales10, align = "c", caption = "Medidas Relevantes para 10 escenarios con 1000 repeticiones de tamaño de muestra variable" )
Medidas Relevantes para 10 escenarios con 1000 repeticiones de tamaño de muestra variable
Escenario n Promedios10 Desviaciones10 Coeficientes_Variacion10 Coeficientes_Asimetria10
1 5 0.1020000 0.1348116 1.3216826 1.2389783
2 10 0.0996000 0.0947038 0.9508413 0.8781349
3 15 0.0996667 0.0771525 0.7741053 0.5549195
4 20 0.1000000 0.0668166 0.6681665 0.5111578
5 30 0.0996667 0.0540425 0.5422320 0.4053367
6 50 0.1008800 0.0429543 0.4257959 0.3111297
7 60 0.0987333 0.0366024 0.3707197 0.4382467
8 100 0.1000100 0.0281548 0.2815198 0.2598761
9 200 0.1001400 0.0185479 0.1852196 0.1502272
10 350 0.0996171 0.0133548 0.1340611 0.2509273
11 500 0.1001080 0.0095699 0.0955962 0.0170250
12 1000 0.1000000 0.0000000 0.0000000 NaN 
cor(Datos_Generales10$n,Datos_Generales10$Desviaciones10)
## [1] -0.6645037
cor(Datos_Generales10$n,Datos_Generales10$Coeficientes_Variacion10)
## [1] -0.669693
require(ggplot2)
require(plotly)

Dispersion1_10=ggplot(Datos_Generales10,aes(x=n,y=Coeficientes_Variacion10)) + geom_point() + geom_smooth(method="loess")
ggplotly(Dispersion1_10)
Dispersion2_10=ggplot(Datos_Generales10,aes(x=n,y=Desviaciones10)) + geom_point() + geom_smooth(method="loess")
ggplotly(Dispersion2_10)

Lotes de 90% Enfermas

Generación del Lote
simulacion10=sample(c(rep("Enferma",900),rep("Sana",100)))
Función para una muestra de tamaño n
funcion10=function(n10){
  b10=sample(simulacion10,size=n10)
  p_circunflejo10=sum(b10=="Enferma")/n10
return(p_circunflejo10)
}
funcion10(12)
## [1] 0.9166667
1000 repeticiones con n = 12
ciclo10=function(n10){
vector10=c(1:n10)
for (i in 1:1000){
  b10=sample(simulacion10,size=n10)
  p_circunflejo10=sum(b10=="Enferma")/n10
vector10[i]=p_circunflejo10
}
return(vector10)}
coef_var10 = function(vector10) {
  sd(vector10) / mean(vector10)
}

n10=12
Ejecucion10=(ciclo10(n10))
library(psych)
promedio10=mean(Ejecucion10)
desviacion10=sd(Ejecucion10)
coeficiente_variacion10=coef_var10(Ejecucion10)
coeficiente_asimetria10=skew(Ejecucion10,type=1)
kurtosis10=kurtosi(Ejecucion10)
hist(Ejecucion10,main=paste("Histograma del valor de p circunflejo para muestras de tamaño",n10),ylab="Frecuencia",xlab="Valor de p circunflejo")

Datos10=data.frame(promedio10,desviacion10,coeficiente_variacion10,coeficiente_asimetria10,kurtosis10)
knitr::kable(Datos10, align = "c", caption = "Medidas Relevantes" )
Medidas Relevantes
promedio10 desviacion10 coeficiente_variacion10 coeficiente_asimetria10 kurtosis10
0.8974167 0.0866314 0.0965342 -0.7046527 0.2313587
1000 repeticiones por cada valor de n
N_Escenarios10=12
Ejecuciones10=data.frame(1:N_Escenarios10,c(5,10,15,20,30,50,60,100,200,350,500,1000))
colnames(Ejecuciones10)=c("Escenario","n")
Ejecuciones10
Escenario n
1 5
2 10
3 15
4 20
5 30
6 50
7 60
8 100
9 200
10 350
11 500
12 1000 
for (i in 1:N_Escenarios10){
  funcion10(Ejecuciones10$n[i])
  print(funcion10(Ejecuciones10$n[i]))
}
## [1] 0.8
## [1] 0.8
## [1] 1
## [1] 0.85
## [1] 0.8333333
## [1] 0.88
## [1] 0.9
## [1] 0.92
## [1] 0.885
## [1] 0.9
## [1] 0.9
## [1] 0.9
#Definición de vectores para contener resultados.

Promedios10=c(1:N_Escenarios10)
Desviaciones10=c(1:N_Escenarios10)
Coeficientes_Variacion10=c(1:N_Escenarios10)
Coeficientes_Asimetria10=c(1:N_Escenarios10)

#Ejecución del ciclo for para cada escenario.

Conf2x5 = matrix(c(1:12), nrow=3, byrow=FALSE)
layout(Conf2x5)
#layout.show(12)
for (i in 1:N_Escenarios10){
  Ejecucion10=(ciclo10(Ejecuciones10$n[i]))
Promedios10[i]=mean(Ejecucion10)
Desviaciones10[i]=sd(Ejecucion10)
Coeficientes_Variacion10[i]=coef_var10(Ejecucion10)
Coeficientes_Asimetria10[i]=skew(Ejecucion10,type=1)
hist(Ejecucion10,main=paste("n=",Ejecuciones10$n[i]),ylab="Frecuencia",xlab="Valor de p")
}

Conf2x5 = matrix(c(1:12), nrow=3, byrow=FALSE)
for (i in 1:N_Escenarios){
  Ejecucion10=(ciclo10(Ejecuciones$n[i]))
qqnorm(Ejecucion10,xlab="Cuantiles Teóricos",ylab="Cuant. Muestrales",main=paste("n=",Ejecuciones$n[i]))
qqline(Ejecucion10, col="red", lwd=1)

}

Datos_Generales10=data.frame(Ejecuciones10,Promedios10,Desviaciones10,Coeficientes_Variacion10,Coeficientes_Asimetria10)
knitr::kable(Datos_Generales10, align = "c", caption = "Medidas Relevantes para 10 escenarios con 1000 repeticiones de tamaño de muestra variable" )
Medidas Relevantes para 10 escenarios con 1000 repeticiones de tamaño de muestra variable
Escenario n Promedios10 Desviaciones10 Coeficientes_Variacion10 Coeficientes_Asimetria10
1 5 0.9002000 0.1357143 0.1507602 -1.2440616
2 10 0.8984000 0.0944794 0.1051641 -0.7661765
3 15 0.9022000 0.0766589 0.0849689 -0.6024209
4 20 0.8985000 0.0676929 0.0753399 -0.6841419
5 30 0.9014000 0.0548629 0.0608641 -0.5038397
6 50 0.8969200 0.0403205 0.0449544 -0.3391170
7 60 0.8998667 0.0360810 0.0400960 -0.3829591
8 100 0.9005200 0.0288229 0.0320070 -0.0951453
9 200 0.8991650 0.0189707 0.0210981 -0.2056489
10 350 0.9003771 0.0128612 0.0142843 -0.0917307
11 500 0.8998620 0.0097931 0.0108829 0.1542631
12 1000 0.9000000 0.0000000 0.0000000 NaN 
cor(Datos_Generales10$n,Datos_Generales10$Desviaciones10)
## [1] -0.6609705
cor(Datos_Generales10$n,Datos_Generales10$Coeficientes_Variacion10)
## [1] -0.6610883
require(ggplot2)
require(plotly)

Dispersion1_10=ggplot(Datos_Generales10,aes(x=n,y=Coeficientes_Variacion10)) + geom_point() + geom_smooth(method="loess")
ggplotly(Dispersion1_10)
Dispersion2_10=ggplot(Datos_Generales10,aes(x=n,y=Desviaciones10)) + geom_point() + geom_smooth(method="loess")
ggplotly(Dispersion2_10)

Conclusiones

En los 3 escenarios simulados, se observa que el valor del estimador \(\hat{p}\) es insesgado con respecto a \({p}\) de cada escenario, es decir, es bastante cercano o idéntico, sin importar el tamaño de la muestra \({n}\); no obstante, no ocurre lo mismo con la varibilidad, ya que el tamaño de muestra \({n}\) se observa inversamente proporcional a la desviación estándar de \(\hat{p}\), así como al coeficiente de variación.

Es por ello que se puede apreciar en los distintos conjuntos de histogramas, que para mayor valor de \({n}\), la gráfica se asemeja más a una distribución normal; además, los límites inferior y superior de \(\hat{p}\) se vuelven más cercanos al parámetro \(p\) a medida que se aumenta el valor de \(n\), en la medida que la variación disminuye.

Por otra parte, es interesante observar que los coeficientes de asimetría para el escenario con porcentaje de plantas uniforme (\(p=0.5\)) posee valores positivos y negativos indistintamente, mientras que en el escenario de 10% enfermas los coeficientes son principalmente positivos, y en el escenario de 90% enfermas estos son principalmente negativos, lo cual también se puede percibir en el desplazamiento de las colas de los histogramas para cada escenario.

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

Literal 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).

Solución

Se generan dos Lotes con las condiciones dadas empleando la función rep para repetir el porcentaje dado el número de plantas enfermas en el lote, y la misma función para el complemento correspondiente a plantas sanas:

Lote1 = sample(c(rep("Enferma",1000*0.1),rep("Sana",1000*0.9)))
Lote2 = sample(c(rep("Enferma",1500*0.1),rep("Sana",1500*0.9)))

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

Solución

Se crea una función que toma como variable entrada el valor de \(n\) para generar dos muestras del mismo tamaño en cada lote, y encontrar la proporción de plantas enfermas en lote dado con dicho tamaño de muestra; luego, calcula la diferencia entre las proporciones resultantes.

muestra=function(n){
  muestra1=sample(Lote1,size=n)
  p_circunflejo_1=sum(muestra1=="Enferma")/n
  muestra2=sample(Lote2,size=n)
  p_circunflejo_2=sum(muestra2=="Enferma")/n
  diferencia=p_circunflejo_1-p_circunflejo_2
return(diferencia)
}

muestra(100)
## [1] -0.05

Literal C

Repita el escenario anterior (b) 500 veces y analice los resultados en cuanto alcomportamiento de los 500 estimadores (diferencias p1-p2). ¿Qué tan simétricos son los datos?, ¿Son siempre cero las diferencias?

Solución

Se genera la función que toma como variable de entrada el tamaño de muestra deseado \(n\), crea dos vectores con \(n\) registros, y a través de un ciclo con contador \(i\) de 500 iteraciones calcula 500 muestras de tamaño \(n\) para cada lote, con su respectivo estimador \(\hat{p}\), calculado como la suma de los valores en que la planta de la muestra de cada lote está enferma, sobre el número total de observaciones \(n\), que corresponde al porcentaje de plantas enfermas en la muestra, luego almacena el resultado de cada iteración en la posición i de cada vector anteriormente definido, y finalmente crea un vector que calcula la diferencia entre los registros por pares de los vectores previamente definidos.

Se logra apreciar en el histograma que los datos tienden a una distribución normal, y de igual manera, al calcular el coeficiente de asimetría por medio de la función skewness, se observa que bastante cercano a 0, por lo cual la distribución es bastante simétrica; además, se plantea una prueba de hipótesis con el test de Shapiro-Wilks, y se logra concluir que, al ser el valor p bastante superior al nivel de confianza de 0.05, no se rechaza la hipótesis nula de que la distribución obtenida sea normal.

\(H_0\): La distribución es normal.

\(H_1\): La distribución no es normal.

Sin embargo, las diferencias no son siempre 0, y existen iteraciones en las cuales los valores están por encima o por debajo de 0, no obstante, estas se compensan y es por ello que el promedio de los valores es bastante cercano a 0.

repeticion=function(n){
vector1=c(1:n)
vector2=c(1:n)
for (i in 1:500){
  muestra1=sample(Lote1,size=n)
  p_circunflejo_1=sum(muestra1=="Enferma")/n
  muestra2=sample(Lote2,size=n)
  p_circunflejo_2=sum(muestra2=="Enferma")/n
vector1[i]=p_circunflejo_1
vector2[i]=p_circunflejo_2}

diferencia=c(vector1-vector2)
return(diferencia)}

resultado=repeticion(100)
hist(resultado,main="Histograma de la diferencia de p para 500 iteraciones con n=100",xlab="Diferencia de p",ylab="Frecuencia")

mean(resultado)
## [1] -0.00162
sd(resultado)
## [1] 0.04129147
library(moments)
skewness(resultado)
## [1] 0.1142026
shapiro.test(resultado)
## 
##  Shapiro-Wilk normality test
## 
## data:  resultado
## W = 0.99216, p-value = 0.009852

Literal 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 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?

Para esto, se emplea nuevamente una metodología similar a la empleada en el Enunciado 1, cambiando la función requerida como argumento de los cálculos y gráficos.

N_Escenarios=12
Ejecuciones=data.frame(1:N_Escenarios,c(5,10,15,20,30,50,60,100,200,350,500,1000))
colnames(Ejecuciones)=c("Escenario","n")
Ejecuciones
Escenario n
1 5
2 10
3 15
4 20
5 30
6 50
7 60
8 100
9 200
10 350
11 500
12 1000 
for (i in 1:N_Escenarios){
  muestra(Ejecuciones$n[i])
  print(muestra(Ejecuciones$n[i]))
}
## [1] 0.2
## [1] 0
## [1] 0
## [1] -0.05
## [1] -0.03333333
## [1] 0.02
## [1] 0
## [1] -0.04
## [1] -0.025
## [1] -0.02857143
## [1] 0
## [1] -0.003
#Definición de vectores para contener resultados.

Promedios=c(1:N_Escenarios)
Desviaciones=c(1:N_Escenarios)
Coeficientes_Asimetria=c(1:N_Escenarios)

#Ejecución del ciclo for para cada escenario.
library(psych)

Conf2x5 = matrix(c(1:12), nrow=3, byrow=FALSE)
layout(Conf2x5)
for (i in 1:N_Escenarios){
  Ejecucion=(repeticion(Ejecuciones$n[i]))
Promedios[i]=mean(Ejecucion)
Desviaciones[i]=sd(Ejecucion)
hist(Ejecucion,main=paste("n=",Ejecuciones$n[i]),ylab="Frecuencia",xlab="Diferencia de p")
abline(v=Promedios[i],col="red",lwd=2)
}

Conf2x5 = matrix(c(1:12), nrow=3, byrow=FALSE)
for (i in 1:N_Escenarios){
  Ejecucion=(repeticion(Ejecuciones$n[i]))
qqnorm(Ejecucion,xlab="Cuantiles Teóricos",ylab="Cuant. Muestrales",main=paste("n=",Ejecuciones$n[i]))
qqline(Ejecucion, col="red", lwd=1)

}

kurtosi(repeticion(1000))
## [1] 3.797166

Se puede apreciar en los histogramas, que para mayor valor de \(n\), la gráfica se asemeja más a una distribución normal; además, los límites inferior y superior de la diferencia de \(\hat{p}\) se vuelven más cercanos a 0 a medida que se aumenta el valor de \(n\). Igualmente, la gráfica de qq-normalidad muestra un comportamiento normal en la mayor parte de los escenarios, excepto para n=1000, ya que se puede apreciar en el último histograma, que hay un apuntamiento sustancial para el valor de diferencia de \(\hat{p}\) igual a 0, y una poca concentración de valores en las colas, lo cual se verifica al calcular el coeficiente de kurtosis mayor a 0 que deja en evidencia un comportamiento leptocúrtico de la distribución; esto sucede porque el tamaño total de la población del lote 1 es igual al tamaño de muestra para este escenario, por lo cual, la probabilidad de obtener valores de diferencia distintos a 0 es mucho menor que en los escenarios anteriores.

También se presenta una tabla con distintos cálculos relevantes para cada escenario:

Datos_Generales=data.frame(Ejecuciones,Promedios,Desviaciones)
knitr::kable(Datos_Generales, align = "c", caption = "Medidas Relevantes para 10 escenarios con 1000 repeticiones de tamaño de muestra variable" )
Medidas Relevantes para 10 escenarios con 1000 repeticiones de tamaño de muestra variable
Escenario n Promedios Desviaciones
1 5 -0.0032000 0.1915807
2 10 -0.0046000 0.1305102
3 15 0.0034667 0.1134330
4 20 0.0052000 0.0944500
5 30 0.0048667 0.0749429
6 50 0.0002400 0.0602062
7 60 0.0037667 0.0529337
8 100 -0.0013000 0.0408589
9 200 -0.0002700 0.0266533
10 350 -0.0018857 0.0193190
11 500 -0.0004760 0.0138988
12 1000 -0.0000320 0.0038897

Literal 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: diferenciade 5%)

Solución

Generación de Lote
Lote1 = sample(c(rep("Enferma",1000*0.1),rep("Sana",1000*0.9)))
Lote2 = sample(c(rep("Enferma",1500*0.15),rep("Sana",1500*0.85)))
Función para una muestra de tamaño n
muestra=function(n){
  muestra1=sample(Lote1,size=n)
  p_circunflejo_1=sum(muestra1=="Enferma")/n
  muestra2=sample(Lote2,size=n)
  p_circunflejo_2=sum(muestra2=="Enferma")/n
  diferencia=p_circunflejo_1-p_circunflejo_2
return(diferencia)
}

muestra(100)
## [1] 0.01
500 Iteraciones con n = 100
repeticion=function(n){
vector1=c(1:n)
vector2=c(1:n)
for (i in 1:500){
  muestra1=sample(Lote1,size=n)
  p_circunflejo_1=sum(muestra1=="Enferma")/n
  muestra2=sample(Lote2,size=n)
  p_circunflejo_2=sum(muestra2=="Enferma")/n
vector1[i]=p_circunflejo_1
vector2[i]=p_circunflejo_2}

diferencia=c(vector1-vector2)
return(diferencia)}

resultado=repeticion(100)
hist(resultado,main="Histograma de la diferencia de p para 500 iteraciones con n=100",xlab="Diferencia de p",ylab="Frecuencia")
abline(v=mean(resultado),col="red",lwd="4")

sd(resultado)
## [1] 0.04366908
library(moments)
skewness(resultado)
## [1] -0.06833376
shapiro.test(resultado)
## 
##  Shapiro-Wilk normality test
## 
## data:  resultado
## W = 0.99383, p-value = 0.03967
500 Iteraciones con n variable
N_Escenarios=12
Ejecuciones=data.frame(1:N_Escenarios,c(5,10,15,20,30,50,60,100,200,350,500,1000))
colnames(Ejecuciones)=c("Escenario","n")
Ejecuciones
Escenario n
1 5
2 10
3 15
4 20
5 30
6 50
7 60
8 100
9 200
10 350
11 500
12 1000 
for (i in 1:N_Escenarios){
  muestra(Ejecuciones$n[i])
  print(muestra(Ejecuciones$n[i]))
}
## [1] 0
## [1] 0
## [1] 0.1333333
## [1] -0.2
## [1] -0.1
## [1] -0.02
## [1] 0
## [1] -0.06
## [1] -0.11
## [1] -0.05428571
## [1] -0.028
## [1] -0.055
#Definición de vectores para contener resultados.

Promedios=c(1:N_Escenarios)
Desviaciones=c(1:N_Escenarios)
Coeficientes_Asimetria=c(1:N_Escenarios)

#Ejecución del ciclo for para cada escenario.
library(psych)

Conf2x5 = matrix(c(1:12), nrow=3, byrow=FALSE)
layout(Conf2x5)
for (i in 1:N_Escenarios){
  Ejecucion=(repeticion(Ejecuciones$n[i]))
Promedios[i]=mean(Ejecucion)
Desviaciones[i]=sd(Ejecucion)
hist(Ejecucion,main=paste("n=",Ejecuciones$n[i]),ylab="Frecuencia",xlab="Diferencia de p")
abline(v=Promedios[i],col="red",lwd=2)
}

Conf2x5 = matrix(c(1:12), nrow=3, byrow=FALSE)
for (i in 1:N_Escenarios){
  Ejecucion=(repeticion(Ejecuciones$n[i]))
qqnorm(Ejecucion,xlab="Cuantiles Teóricos",ylab="Cuant. Muestrales",main=paste("n=",Ejecuciones$n[i]))
qqline(Ejecucion, col="red", lwd=1)

}

kurtosi(repeticion(1000))
## [1] -1.88869

Enunciado 3