Качество подгонки для логит-модели с фиксированным эффектом с использованием пакета «bife»

Пусть ДГП будет

n <- 1000
x <- rnorm(n)
id <- rep(1:2, each = n / 2)
y <- 1 * (rnorm(n) > 0)

так что мы будем в рамках нулевой гипотезы. Как говорится в ?bife, когда нет коррекции смещения, все то же самое, что и в glm, кроме скорости. Итак, начнем с glm.

modGLM <- glm(y ~ 1 + x + factor(id), family = binomial())
modGLM0 <- glm(y ~ 1, family = binomial())

Один из способов выполнить тест LR с

library(lmtest)
lrtest(modGLM0, modGLM)
# Likelihood ratio test
#
# Model 1: y ~ 1
# Model 2: y ~ 1 + x + factor(id)
#   #Df  LogLik Df  Chisq Pr(>Chisq)
# 1   1 -692.70                     
# 2   3 -692.29  2 0.8063     0.6682

Но мы также можем сделать это вручную,

1 - pchisq(c((-2 * logLik(modGLM0)) - (-2 * logLik(modGLM))),
           modGLM0$df.residual - modGLM$df.residual)
# [1] 0.6682207

Теперь приступим к bife.

library(bife)
modBife <- bife(y ~ x | id)
modBife0 <- bife(y ~ 1 | id)

Здесь modBife — полная спецификация, а modBife0 — только с фиксированными эффектами. Для удобства пусть

logLik.bife <- function(object, ...) object$logl_info$loglik

для извлечения логарифмического правдоподобия. Тогда мы можем сравнить modBife0 с modBife, как в

1 - pchisq((-2 * logLik(modBife0)) - (-2 * logLik(modBife)), length(modBife$par$beta))
# [1] 1

в то время как modGLM0 и modBife можно сравнить, запустив

1 - pchisq(c((-2 * logLik(modGLM0)) - (-2 * logLik(modBife))), 
           length(modBife$par$beta) + length(unique(id)) - 1)
# [1] 0.6682207

что дает тот же результат, что и раньше, хотя с bife у нас по умолчанию есть коррекция смещения.

Наконец, в качестве бонуса, мы можем смоделировать данные и убедиться, что тест работает так, как предполагалось. 1000 итераций ниже показывают, что оба теста (поскольку два теста одинаковы) действительно отклоняются так часто, как должны при нулевом значении.