Estadística Computacional

8.3 Simulación para evaluar el ajuste del modelo

La simulación de datos falsos (fake data) o datos replicados se puede usar para explorar las implicaciones del modelo ajustado. En este ejemplo buscamos comparar datos simulados bajo el modelo ajustado con los datos observados.

8.3.0.1 Ejemplo: modelo Poisson con sobreabundancia de ceros

Se busca estudiar el efecto de pesticidas en el control de cucarachas en departamentos urbanos. Se realiza un experimento deonde se dividen los departamentos en:

grupo de tratamiento (160 deptos.) y
grupo de control (104 deptos.).

En cada departamento se mide el número de cucarachas atrapadas \(y_i\) en un conjunto de trampas. Distintos departamentos tuvieron trampas un número distinto de días, y denotamos por \(u_i\) el número de días-trampa en el i-ésimo departamento. Se propone el siguiente modelo:

\[y_i \sim Poisson(u_iexp(X\beta))\] donde X reprersenta variables explicativas (predictores), en este caso, consisten en el nivel de cucarachas antes del tratamiento (roach1), una variable binaria indicando si se aplica insecticida en el departamento (treatment) y una variable binaria indicando si el edificio es de personas mayor (senior). En R el modelo se ajusta como sigue:

library(arm)
roachdata <- read.csv("data/roachdata.csv", stringsAsFactors = FALSE)
#> Warning in file(file, "rt"): cannot open file 'data/roachdata.csv': No such
#> file or directory
#> Error in file(file, "rt"): cannot open the connection
glm_1 <- glm(y ~ roach1 + treatment + senior, family = poisson, 
  offset = log(exposure2), data = roachdata)
#> Error in is.data.frame(data): object 'roachdata' not found
display(glm_1)
#> Error in display(glm_1): object 'glm_1' not found

¿Qué tan bien se ajusta el modelo a los datos? Para responder esta pregunta simularemos del modelo.

X <- model.matrix(~ roach1 + treatment + senior, family = poisson, 
  data = roachdata)
#> Error in terms.formula(object, data = data): object 'roachdata' not found
n <- nrow(X)
simula_modelo <- function(n_sims = 19, ajuste){
    # simulamos los coeficientes del modelo
    betas <- coef(sim(ajuste, n_sims))
    # calculamos ui*exp(Xb)
    y_hat <- roachdata$exposure2 * exp(X %*% t(betas))
    # creamos una lista con las y_hat de cada simulación
    y_hat_list <- split(y_hat, rep(1:ncol(y_hat), each = nrow(y_hat)))
    # simulamos observaciones
    y_sims <- map_df(y_hat_list, ~rpois(n, .))
    y_sims_df <- bind_cols(X = 1:n, y_sims) %>% 
        gather(sim, y, -X) 
    # código para esconder los datos
    codigo <- sample(1:(n_sims + 1), n_sims + 1)
    sims_datos <- y_sims_df %>% 
        bind_rows(dplyr::select(roachdata, X, y)) %>% 
        mutate(sample = rep(codigo, each = n))
    list(sims_datos = sims_datos, y_sims_df = y_sims_df,
        codigo = codigo[n_sims + 1])
}

sim_1 <- simula_modelo(n_sims = 9, glm_1)
#> Error in sim(ajuste, n_sims): object 'glm_1' not found
ggplot(sim_1$sims_datos, aes(x = y)) +
    geom_histogram(binwidth = 3) + 
    xlim(0, 40) +
    facet_wrap(~ sample, nrow = 2)
#> Error in ggplot(sim_1$sims_datos, aes(x = y)): object 'sim_1' not found
# los datos están en sim_1$codigo

¿En que se diferencían los datos observados de los simulados?

Comparemos el número de ceros de los datos observados y el primer conjunto de datos simulados:

mean(roachdata$y == 0)
#> Error in mean(roachdata$y == 0): object 'roachdata' not found
mean(sim_1$sims_datos$y == 0)
#> Error in mean(sim_1$sims_datos$y == 0): object 'sim_1' not found

Vemos que el 36% de los datos observados hay ceros mientras que en los datos replicados el porcentaje de ceros es cercano a cero.

Podemos pensar en la proporción de ceros como una estadística de prueba, simulamos 1000 conjuntos de datos y calculamos la proporción de ceros:

sim_2 <- simula_modelo(1000, glm_1)
#> Error in sim(ajuste, n_sims): object 'glm_1' not found
# calculamos el porcentaje de ceros en cada conjunto simulado
sims_p_ceros <- sim_2$y_sims_df %>% 
    group_by(sim) %>% 
    summarise(p_ceros = mean(y == 0))
#> Error in eval(lhs, parent, parent): object 'sim_2' not found
min(sims_p_ceros$p_ceros)
#> Error in eval(expr, envir, enclos): object 'sims_p_ceros' not found
max(sims_p_ceros$p_ceros)
#> Error in eval(expr, envir, enclos): object 'sims_p_ceros' not found

Vemos que en el porcentaje de ceros varía entre 0 y 0.008, todos ellos mucho menores a la estadística de prueba \(0.36\).

Ahora veamos que ocurre si ajustamos un modelo Poisson con sobredispersión.

glm_2 <- glm(y ~ roach1 + treatment + senior, family = quasipoisson, 
  offset = log(exposure2), data = roachdata)
#> Error in is.data.frame(data): object 'roachdata' not found
display(glm_2)
#> Error in display(glm_2): object 'glm_2' not found

simula_modelo <- function(n_sims = 19, ajuste){
    # simulamos los coeficientes del modelo
    betas <- coef(sim(ajuste, n_sims))
    # calculamos ui*exp(Xb)
    y_hat <- roachdata$exposure2 * exp(X %*% t(betas))
    # creamos una lista con las y_hat de cada simulación
    y_hat_list <- split(y_hat, rep(1:ncol(y_hat), each = nrow(y_hat)))
    # simulamos observaciones
    y_sims <- map_df(y_hat_list, ~rnegbin(n, ., . / (65.4 - 1)))
    y_sims_df <- bind_cols(X = 1:n, y_sims) %>% 
        gather(sim, y, -X) 
    # código para esconder los datos
    codigo <- sample(1:(n_sims + 1), n_sims + 1)
    sims_datos <- y_sims_df %>% 
        bind_rows(dplyr::select(roachdata, X, y)) %>% 
        mutate(sample = rep(codigo, each = n))
    list(sims_datos = sims_datos, y_sims_df = y_sims_df, codigo = codigo[n_sims + 1])
}

sim_2 <- simula_modelo(n_sims = 9, glm_2)
#> Error in sim(ajuste, n_sims): object 'glm_2' not found
ggplot(sim_2$sims_datos, aes(x = y)) +
    geom_histogram(binwidth = 3) + 
    xlim(0, 40) +
    facet_wrap(~ sample, nrow = 2)
#> Error in ggplot(sim_2$sims_datos, aes(x = y)): object 'sim_2' not found

El panel que corresponde a los datos es:

sim_2$codigo
#> Error in eval(expr, envir, enclos): object 'sim_2' not found

Podemos comparar el número de ceros.

sim_1000 <- simula_modelo(1000, glm_2)
#> Error in sim(ajuste, n_sims): object 'glm_2' not found
# calculamos el porcentaje de ceros en cada conjunto simulado
sims_p_ceros <- sim_1000$y_sims_df %>% 
    group_by(sim) %>% 
    summarise(p_ceros = mean(y == 0))
#> Error in eval(lhs, parent, parent): object 'sim_1000' not found
mean(sims_p_ceros$p_ceros >= 0.36)
#> Error in mean(sims_p_ceros$p_ceros >= 0.36): object 'sims_p_ceros' not found

En este caso el 19% de los datos muestran una proporción de ceros al menos tan alta como la observada.

La simulación de datos falsos no debe ser la única herramienta para evaluar el ajuste de un modelo; sin embargo, es una herramienta útil que nos puede ayudar a detectar desajustes y en caso de revelarlos nos da indicios de porque falla el modelo.