Спецификационный тест Хаусмана для "glmer" от lme4

Question

Я хочу сделать "Фиксированные / случайные модели для обобщенной линейной модели" (family="binomial"), потому что у меня есть база данных, где наблюдения происходят из популяции, и существует структура группировки. Тогда я использую функцию glmer от lme4 пакет, также я прочитал, что я могу использовать glmmPQL функция из библиотеки MASS (Далекий,2006).

Моя проблема возникает, когда я хочу обосновать использование случайной или фиксированной модели с помощью теста Хаусмана (Грин,2012), я не нахожу конкретной функции, которая позволяет мне делать это подобно phtest тест, включенный в пакет plm,

Как обосновать использование случайной модели?

Answer 1

Это прямая настройка plm::phtest функция. Я прокомментировал только те строки кода, которые я на самом деле изменил. ИСПОЛЬЗУЙТЕ НА СВОЙ СОБСТВЕННЫЙ РИСК и, если возможно, проверяйте его на соответствие результатам plm::phtest, Я только что адаптировал код, не задумываясь о том, действительно ли он работает правильно!

phtest_glmer <- function (glmerMod, glmMod, ...)  {  ## changed function call
    coef.wi <- coef(glmMod)
    coef.re <- fixef(glmerMod)  ## changed coef() to fixef() for glmer
    vcov.wi <- vcov(glmMod)
    vcov.re <- vcov(glmerMod)
    names.wi <- names(coef.wi)
    names.re <- names(coef.re)
    coef.h <- names.re[names.re %in% names.wi]
    dbeta <- coef.wi[coef.h] - coef.re[coef.h]
    df <- length(dbeta)
    dvcov <- vcov.re[coef.h, coef.h] - vcov.wi[coef.h, coef.h]
    stat <- abs(t(dbeta) %*% as.matrix(solve(dvcov)) %*% dbeta)  ## added as.matrix()
    pval <- pchisq(stat, df = df, lower.tail = FALSE)
    names(stat) <- "chisq"
    parameter <- df
    names(parameter) <- "df"
    alternative <- "one model is inconsistent"
    res <- list(statistic = stat, p.value = pval, parameter = parameter, 
        method = "Hausman Test",  alternative = alternative,
                data.name=deparse(getCall(glmerMod)$data))  ## changed
    class(res) <- "htest"
    return(res)
}

Пример:

library(lme4)
gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
                   data = cbpp, family = binomial)
gm0 <- glm(cbind(incidence, size - incidence) ~ period +  herd,
                   data = cbpp, family = binomial)

phtest_glmer(gm1,gm0)
##  Hausman Test
## data:  cbpp
## chisq = 10.2747, df = 4, p-value = 0.03605
## alternative hypothesis: one model is inconsistent

Кажется, это работает для lme модели тоже:

library("nlme")
fm1 <- lme(distance ~ age, data = Orthodont) # random is ~ age
fm0 <- lm(distance ~ age*Subject, data = Orthodont)
phtest_glmer(fm1,fm0)

## Hausman Test 
## data:  Orthodont
## chisq = 0, df = 2, p-value = 1
## alternative hypothesis: one model is inconsistent