Introducción

En esta clase vamos a hablar de trampas estadísitcas que consciente o inconscientemente podemos cometer en nuestro análisis, generando resultados irreplicables. Nos vamos a enfocar en las consecuencias de estas trampas y cómo evitamos caer en ellas.

Figura del artículo de Baker (2016).

Error tipo I

Empecemos con un pequeño repaso del (principal) problema que queremos evitar. Queremos evitar concluir que hay un mecanismo/patrón biológico, cuando en realidad no lo hay

Veamos un ejemplo de la vida real para distinguir error tipo I de error tipo II

¿Qué tan común debe ser el error tipo I?

Depende de \(\alpha\). Si alpha = 0.05

error1 <-
  function(# promedio de x
    m1 = 0,
    # desviacion de x
    sd1 = 1,
    # n de x
    n1 = 1000,
    # promedio de y
    m2 = m1,
    # desviacion de y
    sd2 = sd1,
    # n de y
    n2 = n1,
    # n de simulaciones
    N,
    # alpha
    a = 0.05) {
    
    # Crear un vector para los valores de P
    V <- NULL
    
    for (i in 1:N) {
      # Simular dos variables independientes
      x <- rnorm(n = n1, mean = m1, sd = sd1)
      y <- rnorm(n = n2, mean = m2, sd = sd2)
      
      # x explica y?
      M <- lm(y ~ x)
      S <- summary(M)
      
      # extraer valor de P
      P <- S$coefficients[2, 4]
      
      # agregarlo al vector
      V <- append(x = V, values = P)
    }
    
    # proporcion de resultados significativos
    Err <- length(which(V < a)) / N
    
    # graficar el histograma, marcando el valor de a
    H <- hist(V,
              main = "",
              xlab = "Valor de P",
              ylab = "Frequencia")
    abline(v = a)
    
    cat("Tasa de error tipo 1\n")
    return(Err)
    return(H)
  }

error1(a = 0.05, N = 1000)

## Tasa de error tipo 1
## [1] 0.048

 

Ejercicio 1

  1. Corra la función error1() con dos diferentes valores de alpha
  2. Corra la función error1() con dos diferentes valores de N

 

Ahora veamos algunas cosas que pueden hacer que este supuesto no se cumpla…

Trampa 1: Comparaciones indirectas

Esta trampa consiste en hacer una interpretación acerca la diferencia entre dos tratamientos o categorías, cuando en realidad lo que hicimos fue investigar cada una por aparte.

Veámoslo en la práctica. La interpretación sería que el efecto de X sobre Y es diferente entre las categorías a y b.

trampa1 <-
  function(# promedio de x
    m1 = 0,
    # desviacion de x
    sd1 = 1,
    # n de x
    n1 = 1000,
    # promedio del error de y
    m2 = m1,
    # desviacion del error de y
    sd2 = sd1 / 2,
    # efecto de x sobre y, por defecto un valor modesto
    m = 0.05,
    # valor de y cuando x = 0
    t = 0,
    # numero de simulaciones
    N,
    #alpha
    a = 0.05) {
    
    # Crear un vector para los valores de P
    # interacccion treat:X
    V.i <- NULL
    # efecto de x solo en treat = a
    V.a <- NULL
    #efecto de x solo en treat = b
    V.b <- NULL
    
    for (i in 1:N) {
      # simular situacion en que y depende de x
      x <- rnorm(n = n1, mean = m1, sd = sd1)
      y <- t + x * m  + rnorm(n = n1, mean = m2, sd = sd2)
      
      # asignar los valores a los dos tratamientos
      dat <-
        data.frame(treat = rep(c("a", "b"), each = n1 / 2),
                   X = c(x),
                   Y = rep(y))
      
      # el modelo correcto
      M <- lm(Y ~ X * treat, data = dat)
      S <- summary(M)
      
      # los modelos que se hicieron
      dat.a <- dat[1:(n1 / 2), ]
      M.a <-  lm(Y ~ X, dat = dat.a)
      S.a <- summary(M.a)
      
      dat.b <- dat[(n1 / 2 + 1):n1, ]
      M.b <-  lm(Y ~ X, data = dat.b)
      S.b <- summary(M.b)
      
      # extraer valor de P de la interaccion
      P.i <- S$coefficients[4, 4]
      P.a <- S.a$coefficients[2, 4]
      P.b <- S.b$coefficients[2, 4]
      
      # agregarlo al vector
      V.i <- append(x = V.i, values = P.i)
      V.a <- append(x = V.a, values = P.a)
      V.b <- append(x = V.b, values = P.b)
    }
    
    # Cuantos son significativos?
    # error tipo 1 de la interaccion
    Err.i <- length(which(V.i < a)) / N
    
    # proporcion de N en la hay significancia en un solo tratamiento
    Err.t <-
      length(which(V.a < a & V.b >= a | V.a >= a & V.b < a)) / N
    
    Err <-
      data.frame("Tipo" = c("interaccion", "trampa"),
                 "Tasa" = c(Err.i, Err.t))
    cat("Tasas de error \n")
    return(Err)
  }

trampa1(N=1000)
## Tasas de error
Tipo Tasa
interaccion 0.043
trampa 0.476

Ejercicio 2

  1. Investigue si este error depende del tamaño del efecto de X sobre Y (m)
  2. Investigue si este error depende del tamaño de muestra (n1)
  3. Investigue si este error depende del error de medición en Y (sd2)

 

Yo soy muy inteligente

y no cometería un error así, porque es obvio que hay que incluir la interacción en el modelo. Veamos un caso más complejo y más real.

Por ejemplo, en este diseño experimental, estamos viendo el efecto de contaminación (Ct), y competencia (Cp) en una respuesta fisiológica al estrés (R) y en tres especies de pecesitos (E).

 

Ejercicio 3

  1. En este caso, lo correcto es hacer R ~ Ct * Cp para cada especie, o R ~ Ct * Cp * E?
  2. Si su modelo es R ~ Ct * Cp * E, y la interacción Ct:Cp:E es significativa, cómo la interpreta?

 

Trampa 2: Inflación de N

Sigamos con el estudio del ejercicio anterior, pero por ahora vamos a asumir que sólo hay una especie y una densidad:

Tenemos 30 tanques por tratamiento pcada uno con 9 peces, es decir 270 peces por tratamiento cada uno con una medida de estrés (E) asociada.

¿Cuál es la unidad de muestreo en este estudio?

¿Cuáles unidades podrían considerarse independientes entre sí?

Veamos que pasa si pensamos que los peces son unidades de muestreo independientes. Vamos a simular una situación en la que la contaminación aumenta el estrés pero el efecto es modualdo por condiciones específicas en cada tanque.

require(ggplot2)
require(lme4)
trampa2 <- function(# efecto promedio de Ct
  m1 = 1,
  # desviacion de Ct
  sd1 = 1,
  # n de peces por tanque
  n1 = 9,
  # n de tanques
  n2 = 30,
  # efecto promedio de los tanques (Tnq)
  m2 = 0,
  # desviacion de Tnq -> indica cuánto se diferencian los tanques entre sí
  sd2 = 1,
  # promedio del error de E
  m3 = 0,
  # desviacion del error de E -> indica la variación entre peces que no se debe al tanque
  sd3 = sd2 / 2,
  # valor de E cuando Ct = 0
  t0 = 0,
  # numero de simulaciones
  N,
  #alpha
  a = 0.05) {
  # crear vectores para guardar los resutados
  # efecto estimado en 3 modelos
  V.m1 <- NULL
  V.m2 <- NULL
  V.m3 <- NULL
  
  # error estandar del efecto en 3 modelos
  V.s1 <- NULL
  V.s2 <- NULL
  V.s3 <- NULL
  
  # simulemos
  for (i in 1:N) {
    # un mismo efecto para todos los peces en un tanque
    Tnq.0 <- rep(rnorm(n = n2 , mean = m2, sd = sd2), each = n1)
    # Sin Ct E depende solo del tanque y el pez
    E.0 <- t0 + Tnq.0 + rnorm(n = n1 * n2 ,
                              mean = m3,
                              sd = sd3)
    
    # efecto de las contaminacion
    Ct <- rnorm(n = n2 * n1,
                mean = m1,
                sd = sd1)
    # un mismo efecto adicional para todos los peces en un tanque
    Tnq.1 <- rep(rnorm(n = n2 , mean = m2, sd = sd2), each = n1)
    # estres depende del pez, el tanque y la contaminacion
    E.1 <- t0 +  Ct + Tnq.1 + rnorm(n = n1 * n2 ,
                                    mean = m3,
                                    sd = sd3)
    # juntamos todo
    dat1 <- data.frame(
      treat = rep(c("0" , "1"), each = n1 * n2),
      E = c(E.0, E.1),
      Tnq = rep(seq(1:(n2 * 2)), each = n1)
    )
    # modelo con peces como unidades de muestreo
    M.1 <- lm(E ~ treat, data = dat1)
    S.1 <- summary(M.1)
    
    # y si sacamos el promedio del nivel de estres por tanque
    dat2 <- aggregate(dat1$E,
                      by = list (dat1$treat, dat1$Tnq),
                      FUN = mean)
    colnames(dat2) <- c("treat", "Tnq", "E")
    
    # modelo con tanques como unidades de muestreo
    M.2 <- lm(E ~ treat, data = dat2)
    S.2 <- summary(M.2)
    
    # modelo mixto
    M.3 <- lmer(E ~ treat + (1|Tnq), data = dat1)
    S.3 <- summary(M.3)
    
    # guardamos efectos y errores estandar en los vectores
    V.m1 <- append(V.m1, S.1$coefficients[2, 1])
    V.s1 <- append(V.s1, S.1$coefficients[2, 2])
    
    V.m2 <- append(V.m2, S.2$coefficients[2, 1])
    V.s2 <- append(V.s2, S.2$coefficients[2, 2])
    
    V.m3 <- append(V.m3, S.3$coefficients[2, 1])
    V.s3 <- append(V.s3, S.3$coefficients[2, 2])
  }
  
  # juntamos todo todo
  err.df <-
    data.frame(
      modelo = rep(rep(
        c("n.total", "n.tanques", "mixto"), each = N
      ), 2),
      valor = c(V.m1, V.m2, V.m3, V.s1, V.s2, V.s3),
      par = rep(c("promedio", "se"), each = N * 3)
    )
  
  # graficamos los histogramas de cada modelo para el efecto y el error
  p <-
    ggplot(data = err.df, aes(x = valor, colour = modelo, fill = modelo)) + geom_histogram(
      bins = 50,
      position = "stack",
      alpha = 0.5,
      lwd = 0.2
    ) +
    theme_minimal(base_size = 16) + facet_wrap(~ par, scales = "free") +
    scale_colour_manual(values = c("#BDB76B", "#32CD32", "#00416A")) +
    scale_fill_manual(values = c("#BDB76B", "#32CD32", "#00416A"))
  return(p)
}

trampa2(m1 = 1, sd2 = 2, sd3 = 1, N = 1000)

Ejercicio 4

  1. ¿Cuál es la consecuencia de aumentar sd2? ¿Por qué?
  2. ¿Cuáles son los grados de libertad en el modelo 1 y en el modelos 2?

 

Trampa 3: Tamaño de muestra pequeño

En biología nos formamos pensando que un tamaño de muestra prqueño no es deseable. Pero ¿por qué? Tal vez tenemos la intuición de que el poder estadístico (la capacidad de detectar un efecto cuando el efecto existe) está relacionadon con nuestro tamaño de muestra. Hasta ahí vamos bien.

Tal vez también sabemos que un mayor tamaño de muestra nos ayuda a obtener estimados más precisos. Y si comobinamos estas dos afirmaciones podemos concluir que si un efecto es pequeño (i.e. cercano a 0) vamos a necesitar un mayor tamaño de muestra para poder deterctarlo.

Creo que hasta aquí seguimos bien. El problema es que a veces derivamos una conclusión más, que sería algo así:

Esto quiere decir que si yo ya detecté un efecto con un tamaño de muestra pequeño, el efecto real debe ser grande, ya que si hubiera sido pequeño no lo hubiera podido detectar

Este razonamiento nos llevaría a no desconfiar de un resultado significativo obtenido con un tamaño de muestra pequeño y tal vez nos llevaría incluso a creer todavía más en la importancia de este efecto.

¡Veamos el problema!

trampa3 <-
  function(
    # promedio de x
    m1 = 0,
    # desviacion de x
    sd1 = 1,
    # promedio de y
    m2 = m1,
    # desviacion de y
    sd2 = sd1,
     # n de muestra pequenia
    np = 10,
    # n de de muestra grande 
    ng = 30,
    # n de simulaciones
    N,
    # alpha
    a = 0.05) {
    # Crear vectores para los valores de P, R y beta
    V.P  <- NULL
    V.R <- NULL
    V.b <- NULL
    
    # n pequenio
    for (i in 1:N) {
      
      # Simulamos dos variables independientes
      x.p <- rnorm(n = np, mean = m1, sd = sd1)
      y.p <- rnorm(n = np, mean = m2, sd = sd2)
      
      # x explica y?
      M.p <- lm(y.p ~ x.p)
      S.p <- summary(M.p)
      
      # extraemos valor de P
      P.p <- S.p$coefficients[2, 4]
      
      # extraemos R
      R.p <- S.p$adj.r.squared
      
      # extraemos pendiente
      beta.p <- S.p$coefficients[2, 2]
      
      # agregamos los valores a los vectores correspondientes
      V.P <- append(x = V.P, values = P.p)
      V.R <- append(x = V.R, values = R.p)
      V.b <- append(x = V.b, values = beta.p)
    }
    
    # n grande
    for (i in 1:N) {
      
      # Simulamos dos variables independientes
      x.g <- rnorm(n = ng, mean = m1, sd = sd1)
      y.g <- rnorm(n = ng, mean = m2, sd = sd2)
      
      # x explica y?
      M.g <- lm(y.g ~ x.g)
      S.g <- summary(M.g)
      
      # extraemos valor de P
      P.g <- S.g$coefficients[2, 4]
      
      # extraemos R
      R.g <- S.g$adj.r.squared
      
      # extraemos pendiente
      beta.g <- S.g$coefficients[2, 2]
      
      # aagregamos los valores a los vectores correspondientes
      V.P <- append(x = V.P, values = P.g)
      V.R <- append(x = V.R, values = R.g)
      V.b <- append(x = V.b, values = beta.g)
    }
    
    # juntamos todo
    dat <-
      data.frame(
        n = rep(c("pq", "gd"), each = N),
        P = V.P,
        R = V.R,
        beta = V.b
      )
    
    # proporcion de resultados significativos
    Err.p <- length(which(dat$n == "pq" & dat$P < a)) / N
    Err.g <- length(which(dat$n == "gd" & dat$P < a)) / N
    
    # R promedio de falsos positivos
    R.p <- mean(dat$R[which(dat$n == "pq" & dat$P < a)])
    R.g <- mean(dat$R[which(dat$n == "gd" & dat$P < a)])
    
    # beta promedio de falsos positivos
    b.p <- mean(dat$beta[which(dat$n == "pq" & dat$P < a)])
    b.g <- mean(dat$beta[which(dat$n == "gd" & dat$P < a)])
    
    # cuadro para imprimir
    err.df <-
      data.frame(
        par = c("error.1", "R.promedio", "beta.promedio"),
        pq = c(Err.p, R.p, b.p),
        gd = c(Err.g, R.g, b.g)
      )
    
    # agreguemos una columna indicando si la prueba es significativa
    for (i in 1:(2 * N)) {
      ifelse(dat$P[i] >= a, dat$S[i] <-
               "N.S.", dat$S[i] <- "Significativo")
    }
    
    # graficamos los histogramas de R para cada tamanio de muestra y en pruebas N.S. o significativas
    p <-
      ggplot(data = dat, aes(x = R, colour = n, fill = n)) + geom_histogram(
        bins = 50,
        position = "identity",
        alpha = 0.5,
        lwd = 0.2
      ) +
      theme_minimal(base_size = 16) + facet_wrap( ~ S, scales = "free") +
      scale_colour_manual(values = c("#BDB76B", "#32CD32")) +
      scale_fill_manual(values = c("#BDB76B", "#32CD32"))
    p
    
    # resultado
    cat("Comparacion entre tamanios de muestra\n")
    out <- list(err.df, p)
    return(out)
    
  }

trampa3(N = 1000)
## Comparacion entre tamanios de muestra
## [[1]]
##             par      pq      gd
## 1       error.1 0.04500 0.05100
## 2    R.promedio 0.45509 0.16142
## 3 beta.promedio 0.26912 0.17323
## 
## [[2]]

Moraleja

Los falsos postivos con tamaños de muestra pequeños típicamente son efectos grandes.

Ejercicio 5

  1. ¿Qué pasa si aumenta el tamaño de muestra grande a 100? ¿Por qué?
  2. ¿Qué esperaría con respecto a la distribución de valores de R si el tamaño de muestra fuera infinito?

 

Discusión

¿Y entonces cuál es un tamaño de muestra apropiado?

 

Trampa 4: Comparaciones múltiples

Comic: Randall Munroe, xkcd.com / CC BY-NC (https://creativecommons.org/licenses/by-nc/2.5/)

Al inicio de la clase vimos que alpha representa la probabilidad de un efecto espúreo, es decir, un efecto que detectamos estadísticamente aún cuando las variables son independientes. La idea aquí es que cuando hacemos más pruebas con una misma variable o variables del mismo set de datos la probabilidad de que al menos una sea espúrea aumenta.

Veamos cómo.

trampa4 <-
  function(# promedio de x
    m1 = 0,
    # desviacion de x
    sd1 = 1,
    # tamaño de muestra
    n1 = 100,
    # numero de pruebas
    ntest = 5,
    # n de simulaciones
    N,
    # alpha
    a = 0.05) {
    # crear un vector para la siginficancia de la prueba
    V.P  <- NULL
    
    # simulemos
    for (i in 1:N) {
      # variable de respuesta
      y <- rnorm(n = n1, mean = m1, sd = sd1)
      
      # vector para los valores de p
      V.p <- NULL
      
      # un analisis con cada predictor
      for (j in 1:ntest) {
        M <-  lm (y ~ rnorm(n = n1, mean = m1, sd = sd1))
        S <- summary(M)
        V.p <- append(x = V.p, values = S$coefficients[2, 4])
      }
      
      # registremos si hay al menos un resultado significativo
      ifelse(test = any(V.p < a),
             yes = V.P[i] <- 1,
             no = V.P[i] <- 0)
    }
    
    # errir tipo 1 con n comparaciones
    err <- length(which(V.P == 1)) / N
    
    cat("Tasa de error tipo 1\n")
    return(err)
  }
      
trampa4(ntest = 20, N = 1000)
## Tasa de error tipo 1
## [1] 0.635

Ejercicio 6

  1. ¿Qué pasa si aumenta el número de tests realizados? ¿Por qué?
  2. ¿Qué pasa si reduce el valor de alpha? ¿Por qué?

 

Dos preguntas sobresalientes

Rara vez colectamos un set de datos completo para hacer un único modelo estadístico. Si esta es su realidad ud. tiene al menos dos preguntas.

¿Cómo puedo corregir este error?

Existen en realidad varias maneras de lidear con este error y hace unos años la discusión de cómo hacerlo era intensa. Primero, si el estudio es del tipo de “prueba de hipótesis” podría limitarse mucho el error haciendo un número razonable (pequeño) de análisis que respondan la pregunta sin probar todo contra todo.

Si el análisis es más exploratorio o por alguna razón se requiren muchas pruebas (e.g. análsis de expresión genética diferencial en los que hay tantas pruebas como genes). necesitamos algún método para corregir por esta inflación del error.

La literatura es extensa y aquí no vamos a entrar en detalles. Sólo quiero hablar de lo básico. Uno de los primeros métodos es la corrección de Boferroni. Esta corrección consiste en simplemente dividir alpha entre el número de pruebas y es abrumadoramente conservadora. Con esto me refiero a que hace que el error tipo II se dispare.

Otro problema es que la corrección de Bonferroni no reconoce el número de comparaciones significativas. Bonferroni no distingue entre la probabilidad de que al menos una de n comparaciones sea espúrea vs. la probabilidad de que dos o tres o cuatro de las n comparaciones sean simultáneamente espúreas. Los últimos tres escenarios son cada uno menos probable que el anterior.

Dados estos problemas (entre otros) y la existencia de métodos más apropiados, como por ejemplo la tasa de descubrimientos falsos de Benjamini-Hochberg y la alternativa de enfocarnos en tamaños de efecto en lugar de valores de P, es más que suficiente para enterrar a Bonferroni.

¿Qué es “un mismo set de datos”?

Si nos preocupa inflar el error tipo I al hacer muchas comparaciones con el mismo set de datos. No es trivial preguntarse ¿qué es un mismo set de data? ¿serán todos los datos usados en un mismo artículo? ¿todos los datos colectados en un mismo proyecto (que produce n artículos)? ¿Todos los datos colectados en una misma población de estudio o en un mismo laboratorio?

Trampa 5: Reporte selectivo

Esta trampa consiste en correr una serie de análisis, pero mencionar sólo los análisis que dieron resultados significativos, como si los demás análisis nunca se nos hubieran ocurrido. En una forma extrema, esta trampa implica modificar y repetir algún análisis hasta que sea significativo (p-hacking).

Ahora, todxs sabemos que esto es un no no, pero ¿por qué? ¿Cuál es el problema realmente de reportar sólo los análisis que dieron significativos?

Veamos cómo esta trampa se relaciona con la anterior.

trampa5 <-
  function(
    # promedio de x
    m1 = 0,
    # desviacion de x
    sd1 = 1,
    # tamaño de muestra
    n1 = 100,
    # n de simulaciones
    N,
    # alpha
    a = 0.05) {
    # Crear vectores para los valores de P, R y beta
    V.E  <- NULL
    
    # simulemos
    for (i in 1:N) {
      # variable de respuesta
      y <- rnorm(n = n1, mean = m1, sd = sd1)
      #vector con para guardar valores de P
      V.P <- NULL
      j <- 1
      
      # hagamos pruebas con predictores independientes, hasta que alguno de significativo
      repeat {
        M <-  lm (y ~ rnorm(n = n1, mean = m1, sd = sd1))
        S <- summary(M)
        V.P <- append(x = V.P, values = S$coefficients[2, 4])
        if (V.P[j] >= a){
                j <- j +1} else {break}
      }
    len.P <- length(V.P)  
     V.E <- append (x = V.E, values = len.P)  
    }
  
  # grafiquemos el histograma
   hist(V.E,
              main = "",
              xlab = "Número de pruebas para llegar a p < 0.05",
              ylab = "Frequencia")
}
      
trampa5(N = 1000)

Al hacer p-hacking estamos inflando el error tipo 1 porque hay más pruebas sobre la(s) misma(s) variable(s) de respuesta que las que estamos reportando.

Ejercicio 7

  1. ¿Qué otras maneras de hacer p-hacking existen?
  2. ¿Cómo podemos hacer para evitarlas?

 

Trampa 6: Especificación del modelo

Comic: Randall Munroe, xkcd.com / CC BY-NC (https://creativecommons.org/licenses/by-nc/2.5/)

Aún si tomamos todas las precauciones para evitar las trampas anteriores, nuestras inferencias pueden ser salvajemente incorrectas si nuestro modelo es una representación muy pobre de la realidad. Con esto hay que decir que todos los modelos son representaciones incompletas de la realidad. Esa es la idea. Un modelo es útil si podemos simplificar la realidad y aún así descubrir algo verdadero sobre la la naturaleza. Para citar a Goerge E. P. Box “todos los modelos son esencialmente incorrectos, pero algunos son útiles.”

La pregunta más difícil es por supuesto, qué tan incorrecto tiene que ser un modelo para dejar de ser útil. Una pregunta más a nuestro alcance por ahora es por qué es importante evitar un modelo muy incorrecto, cuando hay otro modelo menos incorrecto disponible. Por supuesto, las maneras en que un modelo puede ser incorrecto son infinitas. Entonces aquí vamos a ver sólo un ejemplo.

Nos vamos a enfocar en situaciones en las que érroneamente asumimos una distribución normal en los errores. Recordemos que errores son simplemente las desviaciones entre valores observados y un estimado (como un promedio o una línea de ajuste). Digamos que colectamos datos sobre el número de parásitos en una población de peces. En la vida real es común que muchos peces tengan cero o pocos parásitos, ya que han evolucionado adptaciones para resistirlos, mientras que algunos pocos individuos (tal vez los que ya estaban en una mala condición) tienen muchos parásitos.

Aquí queremos investigar si el número de parásitos aumenta con la temperatura del agua. Asumamos que la mayoría de los lagos tienen una temperatura intermedia, algunos son relativamente fríos y todavía menos son muy fríos. Igual con las temperaturas calientes, cuánto más caliente menos común. Veamos.

require(gridExtra)
  require(DAAG)

trampa6 <- function(
  # intercepto
  a = 0,
  # coeficiente lineal
  b = 2,
  # tamanio de muestra
  n = 50,
  # valor en x para la prediccion
  P = 0.5,
  # numero de simulaciones
  N = 1000) {
  
  # vectores para deviance
  V.D1 <- NULL
  V.D2 <- NULL
  
  # vectores para valores predichos
  V.P1 <- NULL
  V.P2 <- NULL
  
  # simulemos
  for (i in 1:N) {
    # x puede tomar cualquier valor entre 0 y 1
    x <- runif(n = n, min = 0, max = 1)
    # usamos una funcion del paquete DAAG para simular la respuesta de acuerdo con Poisson 
    dat <- poissonsim(x = x, a = a, b = b)
    
    # ahora analizamos los resultados como si los errores fueran normales
    M.1 <- glm(y ~ x, family = "gaussian", data = dat)
    S.1 <- summary.glm(M.1)
    
    # deviance del modelo con el predictor
    V.D1 <- append (x = V.D1, values = S.1$deviance)
    
    # numero de parasitos esperados cuando x = P
    P1 <-
      predict.glm(object = M.1,
                  newdata = data.frame(x = P),
                  type = "response")
    V.P1 <- append (x = V.P1, values = P1)
    
    # ahora analizamos con la estructura de error que corresponde
    M.2 <- glm(y ~ x, family = "poisson", data = dat)
    S.2 <- summary(M.2)
    
    # deviance del modelo con el predictor
    V.D2 <- append (x = V.D2, values = S.2$deviance)
    
    # numero de parasitos esperados cuando x = P
    P2 <-
      predict.glm(object = M.2,
                  newdata = data.frame(x = P),
                  type = "response")
    V.P2 <- append (x = V.P2, values = P2)
  }
  
  # juntamos todo
  err.df <-
    data.frame(
      modelo = rep(c("normal", "poisson"), each = N),
      Deviance = c(V.D1, V.D2),
      Prediction = c(V.P1, V.P2)
    )
 
  #graficamos los datos en la escala original
  p1 <-
    ggplot(data = dat, aes(x = x, y = y)) + geom_point(size = 3,
                                                       color =  "#32CD32",
                                                       alpha = 0.5) +
    theme_minimal(base_size = 16)
  
  # graficamos los histogramas de cada modelo para "deviance" y los valores predichos
  p2 <-
    ggplot(data = err.df, aes(x = Deviance, colour = modelo, fill = modelo)) + geom_histogram(
      bins = 50,
      position = "identity",
      alpha = 0.5,
      lwd = 0.2
    ) +
    theme_minimal(base_size = 16) +
    scale_colour_manual(values = c("#BDB76B", "#32CD32")) +
    scale_fill_manual(values = c("#BDB76B", "#32CD32")) +
    geom_vline(
      xintercept = n - 2,
      size = 1,
      lty = 2,
      colour = "grey20"
    )
  
  p3 <-
    ggplot(data = err.df, aes(x = Prediction, colour = modelo, fill = modelo)) + geom_histogram(
      bins = 50,
      position = "identity",
      alpha = 0.5,
      lwd = 0.2
    ) +
    theme_minimal(base_size = 16) +
    scale_colour_manual(values = c("#BDB76B", "#32CD32")) +
    scale_fill_manual(values = c("#BDB76B", "#32CD32")) +
    xlab("Predicción para x = P")
  
  grid.arrange(p1, p2, p3, ncol = 1)
}

trampa6(N = 1000, b = 5)

Ejercicio 8

  1. ¿Qué pasa al aumentar b? ¿Por qúe?
  2. ¿Qué pasa cuando b = 0? ¿Por qué?

 

Moraleja

Tres preguntas de que preocuparse son:

  1. ¿Cuáles son los supuestos que hace mi análisis?

  2. ¿Qué razones tengo para creer que mis datos no violan esos supuestos?

  3. ¿Qué tan robusto es el modelo a la violación de distintos supuestos?

sessionInfo()
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=es_ES.UTF-8      
##  [2] LC_NUMERIC=C              
##  [3] LC_TIME=es_CR.UTF-8       
##  [4] LC_COLLATE=es_ES.UTF-8    
##  [5] LC_MONETARY=es_CR.UTF-8   
##  [6] LC_MESSAGES=es_ES.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8      
##  [8] LC_NAME=C                 
##  [9] LC_ADDRESS=C              
## [10] LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8
## [12] LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets 
## [6] methods   base     
## 
## other attached packages:
##  [1] DAAG_1.25.4         gridExtra_2.3      
##  [3] lme4_1.1-29         Matrix_1.3-4       
##  [5] ggplot2_3.3.6       plyr_1.8.7         
##  [7] rgbif_3.7.3         taxize_0.9.100     
##  [9] tufte_0.12          rticles_0.24       
## [11] revealjs_0.9        rmdformats_1.0.4   
## [13] kableExtra_1.3.4    rmarkdown_2.14     
## [15] sketchy_1.0.2       remotes_2.4.2      
## [17] leaflet_2.1.1       knitr_1.39         
## [19] xaringanExtra_0.7.0 emo_0.0.0.9000     
## [21] cowsay_0.8.0       
## 
## loaded via a namespace (and not attached):
##   [1] uuid_1.1-0          workflowr_1.7.0    
##   [3] systemfonts_1.0.4   lazyeval_0.2.2     
##   [5] splines_4.1.1       crosstalk_1.2.0    
##   [7] usethis_2.1.6       urltools_1.7.3     
##   [9] digest_0.6.29       foreach_1.5.2      
##  [11] htmltools_0.5.3     rsconnect_0.8.26   
##  [13] fansi_1.0.3         magrittr_2.0.3     
##  [15] memoise_2.0.1       vertical_0.1.0.0000
##  [17] svglite_2.1.0       prettyunits_1.1.1  
##  [19] jpeg_0.1-9          colorspace_2.0-3   
##  [21] rvest_1.0.2         rbibutils_2.2.9    
##  [23] xfun_0.31           dplyr_1.0.9        
##  [25] callr_3.7.0         crayon_1.5.1       
##  [27] jsonlite_1.8.0      zoo_1.8-10         
##  [29] iterators_1.0.14    ape_5.6-2          
##  [31] glue_1.6.2          gtable_0.3.0       
##  [33] emmeans_1.7.4-1     webshot_0.5.3      
##  [35] pkgbuild_1.3.1      scales_1.2.0       
##  [37] oai_0.3.2           mvtnorm_1.1-3      
##  [39] solrium_1.2.0       DBI_1.1.3          
##  [41] Rcpp_1.0.9          viridisLite_0.4.0  
##  [43] xtable_1.8-4        bold_1.2.0         
##  [45] clisymbols_1.2.0    datawizard_0.4.1   
##  [47] htmlwidgets_1.5.4   httr_1.4.3         
##  [49] papaja_0.1.1        RColorBrewer_1.1-3 
##  [51] ellipsis_0.3.2      farver_2.1.1       
##  [53] pkgconfig_2.0.3     reshape_0.8.9      
##  [55] sass_0.4.1          utf8_1.2.2         
##  [57] here_1.0.1          conditionz_0.1.0   
##  [59] crul_1.2.0          labeling_0.4.2     
##  [61] tidyselect_1.1.2    rlang_1.0.4        
##  [63] later_1.3.0         effectsize_0.7.0   
##  [65] munsell_0.5.0       tools_4.1.1        
##  [67] cachem_1.0.6        fortunes_1.5-4     
##  [69] cli_3.3.0           generics_0.1.2     
##  [71] tinylabels_0.2.3    devtools_2.4.3     
##  [73] evaluate_0.15       stringr_1.4.0      
##  [75] fastmap_1.1.0       yaml_2.3.5         
##  [77] processx_3.6.1      fs_1.5.2           
##  [79] purrr_0.3.4         ritis_1.0.0        
##  [81] packrat_0.8.0       rrtools_0.1.5      
##  [83] nlme_3.1-152        whisker_0.4        
##  [85] formatR_1.12        xml2_1.3.3         
##  [87] brio_1.1.3          compiler_4.1.1     
##  [89] rstudioapi_0.13     png_0.1-7          
##  [91] curl_4.3.2          testthat_3.1.4     
##  [93] tibble_3.1.8        bslib_0.3.1        
##  [95] stringi_1.7.8       highr_0.9          
##  [97] ps_1.7.1            parameters_0.18.1  
##  [99] desc_1.4.1          lattice_0.20-44    
## [101] nloptr_2.0.3        vctrs_0.4.1        
## [103] pillar_1.8.0        lifecycle_1.0.1    
## [105] rmsfact_0.0.3       Rdpack_2.4         
## [107] triebeard_0.3.0     jquerylib_0.1.4    
## [109] estimability_1.3    data.table_1.14.2  
## [111] insight_0.17.1      httpuv_1.6.5       
## [113] latticeExtra_0.6-29 R6_2.5.1           
## [115] bookdown_0.27       promises_1.2.0.1   
## [117] sessioninfo_1.2.2   codetools_0.2-18   
## [119] boot_1.3-28         MASS_7.3-54        
## [121] assertthat_0.2.1    pkgload_1.2.4      
## [123] rprojroot_2.0.3     withr_2.5.0        
## [125] httpcode_0.3.0      bayestestR_0.12.1  
## [127] parallel_4.1.1      grid_4.1.1         
## [129] minqa_1.2.4         coda_0.19-4        
## [131] git2r_0.30.1        getPass_0.2-2      
## [133] lubridate_1.8.0