来自 lme4 的 Hausman 对“glmer”的规范测试答案

【问题标题】：Hausman's specification test for "glmer" from lme4来自 lme4 的 Hausman 对“glmer”的规范测试
【发布时间】：2014-07-01 01:34:59
【问题描述】：

我想制作一个“广义线性模型的固定/随机模型”（family="binomial"），因为我有一个数据库，其中观察来自人群并且有一个分组结构。然后我使用lme4 包中的函数glmer，我还读到我可以使用库MASS 中的glmmPQL 函数（Faraway，2006）。

当我想使用 Hausman 检验 (Greene,2012) 证明使用随机模型与固定模型的合理性时，我的问题出现了，我没有找到一个特定的函数可以让我做到这一点，类似于 phtest 测试特色在包中plm。

如何证明使用随机模型的合理性？

【问题讨论】：

您能提供更多信息吗？我找到了Wikipedia page，但我不确定 b_0 和 b_1 代表什么：它们只是固定效应协变量（即不包括分组变量）吗？
感谢您的帮助！我很担心，因为我在 R 中找不到特定的函数，特别是在 lme4 包中执行这个测试，所以我怀疑它是否可以完成。或者，如果相反，我可以证明随机模型最适合数据（anova）。有一个类似的问题：stats.stackexchange.com/questions/58900/fixed-vs-random-effects
我不确定维基百科的解释，所以我给你看格林的定义。 “该检验是基于在无相关假设下，OLS、LSDV 和 FGLS 估计量是一致的，但 OLS 是低效的，而在备择假设下，LSDV 是一致的，但 FGLS 不是。因此，在原假设下，两个估计不应有系统的差异，并且可以基于差异进行测试”（Greene，2012）。
Greene 的示例：“我们从固定效应结果和 B_re 的前九个元素和 V_re（不包括常数项）中检索系数向量和估计渐近协方差矩阵 b_fe 和 V_fe”[I不知道为什么有九个第一元素。] H=(b_fe - B_re)' [V_fe - V_re]^-1 (b_fe - B_re)
anova 不能很好地工作，因为很难证明一个模型嵌套在另一个模型中的要求是正确的（这是应用似然比检验所必需的）。

标签： r random panel

【解决方案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

【讨论】：