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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
