Caso 1.

A lotes con distintas cantidades de moscas (tot) se les aplicó 3 tipos de veneno (veneno) a distintas dosis (dosis) y se contabilizó el número de moscas muertas (muertos). Se intenta conocer cuál es el veneno más efectivo. Dado que la variable respuesta es una proporción (moscas muertas / moscas totales) se propone un modelo lineal generalizado con estructura de errores binomial para datos agregados.

ven <- read.table("veneno.txt", header = TRUE)

# construcción de la variable respuesta
rta <- cbind(ven$muertos, ven$tot-ven$muertos)

#modelo
vfit <- glm(rta ~ veneno + dosis, data = ven, family = binomial(logit))

#significancia según el estadístico de Wald
summary(vfit)

#análisis de la devianza 
anova(vfit, test = "Chisq")

#examen gráfico de los residuos
layout(matrix(1:4, 2, 2))
plot(vfit)
layout(1)

#examen sobre sobredispersión
vfit2 <- glm(rta ~ veneno + dosis, data = ven, family = quasibinomial(logit))
summary(vfit2)

#examen sobre la pertinencia del enlace
LP <- vfit$linear.predictors^2
vfit3 <- glm(rta ~ veneno + dosis + LP, data = ven, family = binomial(logit))
summary(vfit3)

#INTERPRETACIÓN DE PARÁMETROS
exp(vfit$coeff)

# parámetros faltantes:

ven$veneno2 <- relevel(ven$veneno, "R")
reor.vfit <- glm(rta ~ veneno2 + dosis, family = binomial(logit), data = ven)
summary(reor.vfit)
exp(reor.vfit$coeff)
1/exp(reor.vfit$coeff)

Caso 2.

El archivo uta.txt contiene datos de un estudio donde se examinó la presencia de lagartijas del género Uta en 19 islas de Baja California. Se desea probar si la presencia de lagartijas (Uta: 0 = ausentes, 1= presentes) depende de la relación perímetro/área de las islas (PA.ratio).
Dado que los datos son de naturaleza binaria, se propone aplicar un modelo lineal generalizado con estructura de errores binomial, para datos no agregados (regresión logística).

La respuesta es en realidad una categoría (presencia/ausencia), por lo que este es también un problema de clasificación, donde podemos utilizar el modelo para predecir una variable categórica.

library(caret)
library(ggplot2)
library(ROCR)

datos <- read.table("uta.txt", header = TRUE)

# modelo
fit <- glm(Uta ~ PA.ratio, data = datos, family = binomial(logit))

# parámetros y su significancia según el estadístico de Wald
summary(fit)

# interpretación de parámetros.
exp(fit$coefficients)

# notar cómo se vuelve más interpretable el intercepto al 
# centrar la variable.
datos$PA.ratio.2 <- datos$PA.ratio - mean(datos$PA.ratio, na.rm = T)

fit2 <- glm(Uta ~ PA.ratio.2, data = datos, family = binomial(logit))
summary(fit2)
exp(fit2$coefficients)

# análisis de la devianza 
anova(fit2, test = "Chisq")

# examen gráfico de los residuos (?)
layout(matrix(1:4, 2, 2))
plot(fit2)
layout(1)

# examen sobre el enlace
LP <- fit2$linear.predictors^2
fit3 <- glm(Uta ~ PA.ratio.2 + LP, data = datos, family = binomial(logit))
summary(fit3)

#gráfico con predict
X <- seq(-20, 64, 0.5)
Y <- predict(fit2, data.frame(PA.ratio.2 = X), type = "response")
plot(Uta ~ PA.ratio.2, data = datos, xlab = "perímetro/área", 
     ylab = "presencia de Uta")
points(X, Y, type = "l")

# gráfico en ggplot2
g1 <- ggplot(datos, aes(x = PA.ratio.2, y = Uta)) + geom_point() +
  stat_smooth(method = "glm", color = "green", se = TRUE,
                method.args = list(family=binomial))
g1

#matriz de confusión
obs <- datos$Uta
esp <- as.numeric(predict(fit2, type="response") > 0.5)
table(esp, obs)

confusionMatrix(table(esp, obs)) # caret

# curva ROC sobre los datos de entrenamiento (!)
# (pocos datos para que sea informativo, sólo como ejemplo)
pr <- prediction(esp, datos$Uta)
prf <- performance(pr, measure = "tpr", x.measure = "fpr")
plot(prf)

Ejercicios

  1. El archivo budworm.txt contiene los resultados de un experimento donde a lotes de 20 gusanos del tabaco (Number) machos o hembras (Gender) se les aplicó un veneno en distintas concentraciones (Dose) y se examinó cuántos se morían (Killed). ¿Existe un efecto de la dosis de veneno sobre la mortalidad de los gusanos? ¿Existe un efecto del sexo? ¿Modifica el sexo del gusano su respuesta a dosis crecientes de veneno? ¿A qué dosis se alcanza la Dosis Letal 50, para cada sexo? Fuente: Dunn PK, Smyth GK. 2018. Generalized Linear Models With Examples in R. Springer. (a través del paquete GLMsData)

  2. El archivo toxo.txt muestra el número de individuos afectados por toxoplasmosis (Infected), el número de individuos examinados (Sampled) y la precipitación anual en 34 ciudades de El Salvador. ¿Hay un efecto de la lluvia sobre la ocurrencia de toxoplasmosis? Examine los gráficos de diagnóstico, ¿es un modelo lineal lo más adecuado? Fuente: Dunn PK, Smyth GK. 2018. Generalized Linear Models With Examples in R. Springer. (a través del paquete GLMsData)

  3. Supongamos una enfermedad que tiene baja incidencia (1 en 10000). Para realizar un estudio sobre los factores ambientales que predisponen a esa enfermedad se localizaron 20 individuos afectados en una ciudad y se tomaron como control 50 individuos elegidos de forma aleatoria en la misma ciudad. ¿Cuáles son las consecuencias sobre el modelos estadístico de este diseño?

  4. En una revisión reciente de modelos logísticos Ives y Garland afirman que datos (o simulaciones) donde más de 7/8 de las simulaciones sean 0 o 1 “contienen poca información, y un investigador prudente no debería analizarlos en primer lugar”. Explique los motivos. (Ayuda: utilice/visualice tablas de confusión).

  5. ¿Tiene usted datos de naturaleza Binomial? (n > 50). Discuta con sus compañeres de mesa, prestando especial atención a: a) el proceso ue originaron los datos, b) distinga la naturaleza agregada o no agregada de los datos. c) Presente sus datos con gráficos y análisis. d) En caso de datos de naturaleza binaria, divida el set de datos en dos para realizar matrices de confusión sobre datos de entrenamiento y datos de prueba.