[R] F-Tests in generalized linear mixed models (GLMM)

[Björn Stollenwerk]

Hi!

I would like to perform an F-Test over more than one variable within a 
generalized mixed model with Gamma-distribution
and log-link function. For this purpose, I use the package mgcv.
Similar tests may be done using the function "anova", as for example in 
the case of a normal
distributed response. However, if I do so, the error message
"error in eval(expr, envir, enclos) : object "fixed" not found" occures.
Does anyone konw why, or how to fix the problem? To illustrate the 
problem, I send the output of a simulated example.
Thank you very much in advance.

Best regards, Björn

Example:

 > # simulation of data
 > n <- 300
 > x1 <- sample(c(T,F), n, replace=TRUE)
 > x2 <- rnorm(n)
 > random1 <- sample(c("level1","level2","level3"), n, replace=TRUE)
 > true.lp <- 5 + 1.1*x1 + 0.16 * x2
 > mu <- exp(true.lp)
 > sigma <- mu * 1
 > a <- mu^2/sigma^2
 > s <- sigma^2/mu
 > y <- rgamma(n, shape=a, scale=s)
 >
 > library(mgcv)
 >
 > # a mixed model without Gamma-distribution and without log-link works 
as follows:
 > glmm1 <- gamm(y ~ x1 + x2, random=list(random1 = ~1))
 > glmm2 <- gamm(y ~ 1, random=list(random1 = ~1))
 >
 > anova(glmm1$lme)
numDF denDF F-value p-value
X 3 295 103.4730 <.0001
 > anova(glmm2$lme, glmm1$lme)
Model df AIC BIC logLik Test L.Ratio p-value
glmm2$lme 1 3 4340.060 4351.172 -2167.030
glmm1$lme 2 5 4292.517 4311.036 -2141.258 1 vs 2 51.54367 <.0001
 >
 > # a linear model also works, though no p-value is reported
 > glm1 <- gam(y ~ x1 + x2)
 > glm2 <- gam(y ~ 1)
 > anova(glm1)

Family: gaussian
Link function: identity

Formula:
y ~ x1 + x2

Parametric Terms:
df F p-value
x1 1 45.58 7.69e-11
x2 1 13.96 0.000224

 > anova(glm2, glm1)
Analysis of Deviance Table

Model 1: y ~ 1
Model 2: y ~ x1 + x2
Resid. Df Resid. Dev Df Deviance
1 299 33024943
2 297 27811536 2 5213407
 >
 > # general linear models (GLM) with Gamma and log-link don't work
 > glm1.gamma <- gam(y ~ x1 + x2, family=Gamma(link="log"))
 > glm2.gamma <- gam(y ~ 1, family=Gamma(link="log"))
 > anova(glm1.gamma)

Family: Gamma
Link function: log

Formula:
y ~ x1 + x2

Parametric Terms:
df F p-value
x1 1 59.98 1.52e-13
x2 1 16.06 7.78e-05

 > anova(glm2.gamma, glm1.gamma)
Analysis of Deviance Table

Model 1: y ~ 1
Model 2: y ~ x1 + x2
Resid. Df Resid. Dev Df Deviance
1 299 413.59
2 297 343.90 2 69.69
 >
 > # neither do general linear mixed models (GLMM)
 >
 > glm1.gamma <- gamm(y ~ x1 + x2, random=list(random1 = ~1), 
family=Gamma(link="log"))

Maximum number of PQL iterations: 20
iteration 1
 > glm2.gamma <- gamm(y ~ 1, random=list(random1 = ~1), 
family=Gamma(link="log"))

Maximum number of PQL iterations: 20
iteration 1
 > summary(glm1.gamma$lme)
Linear mixed-effects model fit by maximum likelihood
Data: data
AIC BIC logLik
847.722 866.241 -418.861

Random effects:
Formula: ~1 | random1
(Intercept) Residual
StdDev: 2.954058e-05 0.9775214

Variance function:
Structure: fixed weights
Formula: ~invwt
Fixed effects: list(fixed)
Value Std.Error DF t-value p-value
X(Intercept) 5.066376 0.08363392 295 60.57801 0e+00
Xx1TRUE 0.884486 0.11421762 295 7.74387 0e+00
Xx2 0.234537 0.05851689 295 4.00802 1e-04
Correlation:
X(Int) X1TRUE
Xx1TRUE -0.733
Xx2 -0.008 0.085

Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.0207671 -0.6911364 -0.2899184 0.3665161 4.9603830

Number of Observations: 300
Number of Groups: 3
 > summary(glm1.gamma$gam)

Family: Gamma
Link function: log

Formula:
y ~ x1 + x2

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.06638 0.08363 60.578 < 2e-16 ***
x1TRUE 0.88449 0.11422 7.744 1.53e-13 ***
x2 0.23454 0.05852 4.008 7.75e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


R-sq.(adj) = 0.171 Scale est. = 0.95555 n = 300
 >
 > anova(glm1.gamma$lme)
numDF denDF F-value p-value
X 3 295 3187.192 <.0001
 > anova(glm2.gamma$lme, glm1.gamma$lme)
Fehler in eval(expr, envir, enclos) : objekt "fixed" nicht gefunden
 >