Tappin73 Tappin73 - 1 year ago 57
R Question

Hausman's specification test for "glmer" from lme4

First of all remark that I´m learning R lenguage so this question could be basic, but I don´t find anything that can help me.

I would make a "Fixed/Random Models for Gegeneralized linear model" (family="binomial"), because I have a data base where observations come from a population and there are a grouping structure. Then I use the function "glmer" form de "lme4" package, also I have read than can be used the glmmPQL function from library "MASS" (Faraway,2006).

My problem arises when I want to justify the use of random versus fixed model using the Hausman's test (Greene,2012), I don't find a specific function that allows me to do this similar to the "phtest" feature package "plm" test.

It may be a difficult or basic question (depends on how you look at it) but all help will be appreciated or how to justify using the random model!

Thanks and sorry for the format mistakes.

Answer Source

This is a straightforward tweaking of the plm::phtest function. I've commented on the only lines of code that I actually changed. USE AT YOUR OWN RISK and if at all possible test it against the results from plm::phtest. I have just adapted the code, not thought about whether it's really doing the right thing!

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"


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)

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

This seems to work for lme models too:

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

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