Regresión binomial "logit" frente a "cloglog"

Question

Estoy utilizando una regresión binomial con un factor categórico con 9 niveles (denominado "grupo de tratamiento") y tamaños de muestra en cada grupo de 1-8. Un grupo de tratamiento tiene todos los casos positivos (es decir, 1) - y esto crea un problema de estimación con la función "estándar" glm() en R causado por la "separación perfecta" para ese nivel de tratamiento. Por lo tanto, estoy utilizando bayesglm del paquete arm.

Mi pregunta es que el enlace de identidad por defecto es "lógico" pero he leído que "cloglog" o (Log-Log complementario) se utiliza frecuentemente cuando la probabilidad de un evento es muy pequeña o muy grande. Por lo tanto, dado que mi modelo muestra una separación perfecta en un grupo de tratamiento, la probabilidad del evento es muy grande y debería utilizar "cloglog". El uso de cloglog me da un resultado significativo para el grupo de tratamiento con separación perfecta, mientras que "logit" no lo hace.

¿Estoy justificado en el uso de "cloglog" o hay una manera de ver mis resultados y estar seguro de qué enlace es mejor?

f1<- bayesglm(response~ treat.group,family=binomial(link="logit"), data=df)

f2 <- bayesglm(response~ treat.group,family=binomial(link="cloglog"), data=df)

(Cuadro de datos más abajo)

{ structure(list(response = c(0L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L), treat.group = c("pctc", "phth", "phth", "phtl", "pltc", "pcth", "pltl", "phtc", "pctl", "phtc", "pcth", "pctl", "pctc", "phtl", "plth", "pltc", "phtc", "pcth", "phtl", "plth", "pctl", "pltc", "phtl", "pctc", "pcth", "pltc", "phtc", "phtl", "phtc", "pctl", "pctc", "pcth", "phth", "pctc", "phtl", "pcth", "phth", "phtc", "pcth", "phth")), .Names = c("response", "treat.group"), row.names = c(NA, -40L), class = c("tbl_df", "tbl", "data.frame"), na.action = structure(c(1L, 4L, 5L, 7L, 15L, 21L, 23L, 24L, 27L, 29L, 33L, 37L, 39L, 48L, 50L, 53L, 54L, 55L, 56L, 57L, 58L, 59L, 60L, 62L, 63L, 65L, 66L, 67L, 68L, 70L, 71L, 72L, 73L), .Names = c("1", "4", "5", "7", "15", "21", "23", "24", "27", "29", "33", "37", "39", "48", "50", "53", "54", "55", "56", "57", "58", "59", "60", "62", "63", "65", "66", "67", "68", "70", "71", "72", "73"), class = "omit")) }

Answer 1

Si ajusta dos modelos, uno con logit y otro con cloglog, debe informar de los resultados de ambos, y también llevar a cabo algún tipo de técnica de comparación de modelos e informar de los resultados de la misma.

En cuanto a los modelos, se trata de una situación ideal para utilizar modelos multinivel bayesianos (véase este documento [PDF] de Gelman). Podemos agrupar la información entre los grupos para informar de las estimaciones de los grupos con tamaños de muestra más pequeños y de los grupos con resultados aparentemente extremos (como en phth y pltl).

Ajuste de modelos multinivel con el brms es bastante sencillo. Los dos modelos aquí, uno con un enlace logit y otro con un enlace cloglog, se ajustarían como sigue:

fit1 <- brm(
    response ~ (1 | treat.group),
    family = bernoulli,  # logit link is the default
    data = df,
    iter = 2e4,
    warmup = 2e3,
    control = list(adapt_delta = 0.999)
)

fit2 <- brm(
    response ~ (1 | treat.group),
    family = bernoulli(link = "cloglog"),
    data = df,
    iter = 2e4,
    warmup = 2e3,
    control = list(adapt_delta = 0.999)
)

Cada una de ellas generará 72000 muestras de los modelos, a las que se puede acceder con as.data.frame(fit) . Tras aplicar el logit inverso y el cloglog inverso a las muestras, obtenemos las siguientes estimaciones para las probabilidades deseadas:

group    prob using logit    sd       89% probability interval
--------------------------------------------------------------
pctc     0.405               0.169    0.12 to 0.66    
pcth     0.562               0.135    0.34 to 0.78    
pctl     0.438               0.170    0.15 to 0.7    
phtc     0.604               0.141    0.39 to 0.84    
phth     0.729               0.162    0.52 to 1    
phtl     0.529               0.142    0.3 to 0.76    
pltc     0.532               0.157    0.28 to 0.79    
plth     0.539               0.180    0.24 to 0.84    
pltl     0.623               0.193    0.39 to 1

group    prob using cloglog    sd       89% probability interval
----------------------------------------------------------------
pctc     0.395                 0.162    0.12 to 0.64    
pcth     0.545                 0.136    0.32 to 0.76    
pctl     0.422                 0.164    0.14 to 0.67    
phtc     0.589                 0.144    0.36 to 0.83    
phth     0.760                 0.182    0.52 to 1    
phtl     0.511                 0.142    0.28 to 0.74    
pltc     0.512                 0.158    0.25 to 0.77    
plth     0.517                 0.181    0.21 to 0.81    
pltl     0.627                 0.213    0.37 to 1

En general, el modelo que utiliza el enlace cloglog indujo una menor contracción de las probabilidades. Cabe destacar que los intervalos de probabilidad del 89% para phth son los mismos en los dos modelos.

Le site brms admite la comparación de modelos mediante PSIS-LOO . Utilizando el comando

loo(fit1, fit2, reloo = TRUE)

tarda un poco pero finalmente devuelve el siguiente resultado:

            LOOIC   SE
fit1        59.19 2.95
fit2        58.59 2.91
fit1 - fit2  0.60 0.46

Aquí se estima que el modelo con un LOOIC más bajo tiene una mejor precisión de predicción fuera de la muestra. El modelo con el vínculo clog tiene un LOOIC marginalmente inferior, y el error estándar de la diferencia, 0,46, es pequeño, por lo que concluiríamos que el modelo que utiliza el vínculo clog es marginalmente mejor que el que utiliza el vínculo logit.