Tappin73 - 1 year ago 57

R Question

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"
return(res)
}
```

Example:

```
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
```

This seems to work for `lme`

models too:

```
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
```