Taller de Simulación en R

Autor: Jaime Reinoso

PROBLEMA 2

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

Para facilidad de análisis, se generan 3 códigos distintos según cada caso.

Para este primer caso, se generará una población para el grupo 1 y grupo 2 con un 10% de plantas enfermas.

Luego se genera una función que dada una población y un n, obtiene para dicha población una muestra de n elementos y calcula la proporción de plantas enfermas (enfermas/total instancias) y retorna dicho estadístico.

Finalmente, se hace un for para obtener 500 muestras para cada uno de los diferentes valores de n definidos en vn. En cada ciclo, se obtiene una muestra de tamaño n de cada grupo, se calcula el respectivo estadístico de proporción de enfermas, y se calcula la diferencia entre p1 (estadístico calculado sobre la muestra que viene de la población 1) y p2 (idem, pero de la población 2). Se obtiene para este estadístico p1-p2 la prueba Shapiro para determinar si su distribución es normal. También se calcula el histograma correspondiente y el gráfico qq.

estadoPersona = c('Sana','Enferma')

#Definimos los porcentajes de enfermas del groupoA (varía) y el grupo B queda constante con 10%
#vProbEnfermaGrupoA va a ser probado posteriormente con valores 0.05, 0.10,  0.15, pasados por parámetro

#La cantidad de enfermas del grupo B permanecerá constante en 10%
probEnfermaGrupo2 = 0.1

#Se crea la población del grupo B
pobGrupo2 <- sample(estadoPersona, 
                   size=1500, 
                   replace=TRUE,     # se debe permitir pues si no se acaban los valores a escoger
                   prob=c(1-probEnfermaGrupo2, probEnfermaGrupo2)  # 90% sanas y 0% enfermas
                  )

probEnfermaGrupo1 = 0.1
#Se crea la población correspondiente al grupo A que varía según el parámetro de entrada.
pobGrupo1 <- sample(estadoPersona, 
                   size=1000, 
                   replace=TRUE,     # se debe permitir pues si no se acaban los valores a escoger
                   prob=c(1-probEnfermaGrupo1, probEnfermaGrupo1)  # 50% cada tipo, si probEnfermaGrupoA = 0.5
                  )

table(pobGrupo1)
## pobGrupo1
## Enferma    Sana 
##      98     902
table(pobGrupo2)
## pobGrupo2
## Enferma    Sana 
##     157    1343
#  definimnos una funcion, que dada una población y un n, se obtiene una muestra de tamaño n de esa población
# y retorna la proporción de enfermas
Muestra <- function(poblacion, n){
              muestra <- sample(poblacion, size=n, replace=FALSE)  # es un vector de caracteres
           
              #debemos contar cuantos plantas hay Enfermas
              y = length(muestra[muestra== 'Enferma'])
              
              # calculamos la proporción. 
              p = y/n
              return(p)
}
  
  
#Se define el vector de tamaño de muestra
vn<-c(5,10,15,20,30,50,60,100,200,500)   # vn es el vector que tiene los diferentes valores de n a probar
  
#n es el tamaño de muestra a probar
for (n in vn){
    
       t<-c()  # se crea el vector vacío
    
       for (i in 1:500) {      # Se hace un ciclo 500 veces, o sea, se obtienen 500 muestras, cada una de tamaño n.
      
         p1 <- Muestra(pobGrupo1, n)
         p2 <- Muestra(pobGrupo2, n)
      
         t<-c(t,p1-p2)   # Se captura el indicador de diferencia de proporción para ese tamaño de muestra (p1 - p2)
       }
    
      # Se obtiene la prueba shapiro para t para ver si el indicador pA-pB es normal
      pruebaST <-shapiro.test(t)
    
      #Se obtiene el histograma correspondiente.
      hist(t, main = paste(" Porc Enfermas 1: ", probEnfermaGrupo1, "n =", n, "Shap:",format(pruebaST$p.value, digits=4)))  
      qqnorm(t, main = paste(" Porc Enfermas 1: ", probEnfermaGrupo1, "n =", n, "Shap:",format(pruebaST$p.value, digits=4)))
    
} # end for n in vn

Como primera conclusión, se observa que a medida que n tiende a ser más grande, más claro es que la distribución del estadístico p1-p2 tiende a normal. Esto sigue lo esperado por el teorema del límite central visto en clase y las evidencias se observan tanto en la inspección visual del histograma resultando, el resultado de la prueba Shapiro y la inspección visual del diagrama qq.

Analisando el experimento en si (ambas poblaciones con igual % de plantas enfermas del 10%), a medida que n es más y más grande, más claro se observa que la distribución de n tiende a tener un centro en cero, indicando que p1-p2 tendería a cero en la población. Esto es completamente esperado, pues esta simulación parte del hecho que ambas poblaciones tienen un 10% de plantas enfermas, y por tanto p1-p2 va a estar centrado en cero, permitiendo al investigador concluir que no hay diferencias significativas entre ambos tratamientos.

Vamos a comparar ahora con el caso en que el grupo 1 presenta una mejora del 5% frente a enfermas (o sea, que el porcentaje de enfermas se reduce de 10% a 5%), mientras que el grupo 2 continúa con una proporción de plantas enfermas del 10%.

estadoPersona = c('Sana','Enferma')

#Definimos los porcentajes de enfermas del groupoA (varía) y el grupo B queda constante con 10%
#vProbEnfermaGrupoA va a ser probado posteriormente con valores 0.05, 0.10,  0.15, pasados por parámetro

#La cantidad de enfermas del grupo B permanecerá constante en 10%
probEnfermaGrupo2 = 0.1

#Se crea la población del grupo B
pobGrupo2 <- sample(estadoPersona, 
                   size=1500, 
                   replace=TRUE,     # se debe permitir pues si no se acaban los valores a escoger
                   prob=c(1-probEnfermaGrupo2, probEnfermaGrupo2)  # 90% sanas y 0% enfermas
                  )

probEnfermaGrupo1 = 0.05
#Se crea la población correspondiente al grupo A que varía según el parámetro de entrada.
pobGrupo1 <- sample(estadoPersona, 
                   size=1000, 
                   replace=TRUE,     # se debe permitir pues si no se acaban los valores a escoger
                   prob=c(1-probEnfermaGrupo1, probEnfermaGrupo1)  # 50% cada tipo, si probEnfermaGrupoA = 0.5
                  )

table(pobGrupo1)
## pobGrupo1
## Enferma    Sana 
##      56     944
table(pobGrupo2)
## pobGrupo2
## Enferma    Sana 
##     175    1325
#  definimnos una funcion, que dada una población y un n, se obtiene una muestra de tamaño n de esa población
# y retorna la proporción de enfermas
Muestra <- function(poblacion, n){
              muestra <- sample(poblacion, size=n, replace=FALSE)  # es un vector de caracteres
           
              #debemos contar cuantos plantas hay Enfermas
              y = length(muestra[muestra== 'Enferma'])
              
              # calculamos la proporción. 
              p = y/n
              return(p)
}
  
  
#Se define el vector de tamaño de muestra
vn<-c(5,10,15,20,30,50,60,100,200,500)   # vn es el vector que tiene los diferentes valores de n a probar
  
#n es el tamaño de muestra a probar
for (n in vn){
    
       t<-c()  # se crea el vector vacío
    
       for (i in 1:500) {      # Se hace un ciclo 500 veces, o sea, se obtienen 500 muestras, cada una de tamaño n.
      
         p1 <- Muestra(pobGrupo1, n)
         p2 <- Muestra(pobGrupo2, n)
      
         t<-c(t,p1-p2)   # Se captura el indicador de diferencia de proporción para ese tamaño de muestra (p1 - p2)
       }
    
      # Se obtiene la prueba shapiro para t para ver si el indicador pA-pB es normal
      pruebaST <-shapiro.test(t)
    
      #Se obtiene el histograma correspondiente.
      hist(t, main = paste(" Porc Enfermas 1: ", probEnfermaGrupo1, "n =", n, "Shap:",format(pruebaST$p.value, digits=4)))  
      qqnorm(t, main = paste(" Porc Enfermas 1: ", probEnfermaGrupo1, "n =", n, "Shap:",format(pruebaST$p.value, digits=4)))
    
} # end for n in vn

Primero que todo, a medida que n crece, se observa que la distribución de p1-p2 tiende a un valor negativo de -0.5. Esto es consistente con la simulación pues la proporción de plantas enfermas en el grupo 1 tiende a 0.05 y en el grupo 2 sigue en 0.1 p1-p2 por supuesto entonces va a tender a -0.5. El investigador entonces podría concluir que la evidencia parece señalar que el tratamiento 1 tiene un mejor resultado que el tratamiento 2. Esto por supuesto se hace más y más claro a medida que n crece.

Se sigue también cumpliendo lo señalado por el teorema del límite central, en donde la distribución de p1-p2 tiende a ser normal, particularmente a medida que n crece. Esto se confirma tanto por la inspección visual del histograma, el resultado de la prueba Shaphiro y también de la inspección visual del gráfico qq.

Ahora veamos a estudair el caso que ocurre cuando se aumenta la proporción de personas enfermas a 15% en el grupo 1 (pasa de 10% a 15%) y se deja el grupo 2 en 10% de proporción de plantas enfermas.

estadoPersona = c('Sana','Enferma')

#Definimos los porcentajes de enfermas del groupoA (varía) y el grupo B queda constante con 10%
#vProbEnfermaGrupoA va a ser probado posteriormente con valores 0.05, 0.10,  0.15, pasados por parámetro

#La cantidad de enfermas del grupo 2 permanecerá constante en 10%
probEnfermaGrupo2 = 0.1

#Se crea la población del grupo B
pobGrupo2 <- sample(estadoPersona, 
                   size=1500, 
                   replace=TRUE,     # se debe permitir pues si no se acaban los valores a escoger
                   prob=c(1-probEnfermaGrupo2, probEnfermaGrupo2)  # 90% sanas y 0% enfermas
                  )


# el grupo 1 pasa a una proporción de plantas enfermas del 15%
probEnfermaGrupo1 = 0.15

#Se crea la población correspondiente al grupo A que varía según el parámetro de entrada.
pobGrupo1 <- sample(estadoPersona, 
                   size=1000, 
                   replace=TRUE,     # se debe permitir pues si no se acaban los valores a escoger
                   prob=c(1-probEnfermaGrupo1, probEnfermaGrupo1)  # 50% cada tipo, si probEnfermaGrupoA = 0.5
                  )

table(pobGrupo1)
## pobGrupo1
## Enferma    Sana 
##     141     859
table(pobGrupo2)
## pobGrupo2
## Enferma    Sana 
##     158    1342
#  definimnos una funcion, que dada una población y un n, se obtiene una muestra de tamaño n de esa población
# y retorna la proporción de enfermas
Muestra <- function(poblacion, n){
              muestra <- sample(poblacion, size=n, replace=FALSE)  # es un vector de caracteres
           
              #debemos contar cuantos plantas hay Enfermas
              y = length(muestra[muestra== 'Enferma'])
              
              # calculamos la proporción. 
              p = y/n
              return(p)
}
  
  
#Se define el vector de tamaño de muestra
vn<-c(5,10,15,20,30,50,60,100,200,500)   # vn es el vector que tiene los diferentes valores de n a probar
  
#n es el tamaño de muestra a probar
for (n in vn){
    
       t<-c()  # se crea el vector vacío
    
       for (i in 1:500) {      # Se hace un ciclo 500 veces, o sea, se obtienen 500 muestras, cada una de tamaño n.
      
         p1 <- Muestra(pobGrupo1, n)
         p2 <- Muestra(pobGrupo2, n)
      
         t<-c(t,p1-p2)   # Se captura el indicador de diferencia de proporción para ese tamaño de muestra (p1 - p2)
       }
    
      # Se obtiene la prueba shapiro para t para ver si el indicador pA-pB es normal
      pruebaST <-shapiro.test(t)
    
      #Se obtiene el histograma correspondiente.
      hist(t, main = paste(" Porc Enfermas 1: ", probEnfermaGrupo1, "n =", n, "Shap:",format(pruebaST$p.value, digits=4)))  
      qqnorm(t, main = paste(" Porc Enfermas 1: ", probEnfermaGrupo1, "n =", n, "Shap:",format(pruebaST$p.value, digits=4)))
    
} # end for n in vn

En cuanto a la distribución obtenida, se sigue confirmando lo señalado por el teormea del límite central, en que a medida que n crece, la distribución tiende a normal, hecho que se confirma tanto por la inspección visual del histograma de p1-p2 obtenido, el resultado de Shapiro y también de la inspección visual del gráfico qq.

En cuanto al resultado del experimento como tal en donde la población 1 tiene un tasa de plantas enfermas del 15% y el grupo 2 del 10%, se observa que la distribución resultado de p1-p2 tiende a un valor positivo de 0.5 a medida que n crece (explicado por 0.15-0.1 > 0 y tendiente a 0.05). Esto permitiría al investigador señalar que la evidencia parece indicar que el tratamiento en el grupo 2 es más exitoso que en el grupo 1.