Dear R users
I have a question related to the interpretation of results based on GAMMs using
Simon Woods package gamm4.
I have repeated measurements (hours24) of subjects (vpnr) and one factor with
three levels (pred). The outcome (dv) is binary.
In the first model I'd like to test for differences among factor levels (main
effects only):
gamm.11<-gamm4(dv ~ pred +s(hours24), random = ~
(1|vpnr),data=sdata,family=binomial)
In the second model I'd like to test whether the smooths vary across factor
levels (interaction):
gamm.12<-gamm4(dv ~ pred +s(hours24) +s(hours24,by=pred), random = ~
(1|vpnr),data=sdata,family=binomial)
Part of the output for both models is shown below.
My questions are:
1) Which is the best way to test whether smooths actually vary across factor
levels?
Is thought that one way would be to compare the two models, e.g. by calculating
differences between the two deviances and degrees of freedom (as long as this
has a chi-square distribution). If so this would serve as a kind of overall
test for interaction which would be fine.
Would this be suitable or is there a better way to accomplish this?
2) It would be interesting to compare smooths of specific factor levels with
each other.
I thought that I would get the answer to this question by correctly
interpreting the coefficients of the second model. I've assumed that in the
output below (at the very end) the (approximate) significance of the first
smooth term s(hours24) tests for the smooth of the averaged values (across all
three factor levels) against a horizontal line, whereas the three additional
smooths (s(hours24):pred1 etc.) test for the deviation from this average to the
respective factor level. If so how would I have to set up a model in which I
could directly compare specific factors with each other?
I've also tried to omit the smooth term s(hours24):
gamm.12<-gamm4(dv ~ pred +s(hours24,by=pred), random = ~
(1|vpnr),data=sdata,family=binomial)
However this would obviously test each of the smooths aginst a horizontal line
which is not what I want.
Any help is greatly appreciated, thanks!
Andrea
Output of 1st model gamm.11:
----------------------------------------
logLik(gamm.11$mer);deviance(gamm.11$mer);attributes(summary(gamm.11$mer))$AICtab[1];gamm.11$...@deviance["disc"]
'log Lik.' -49054.65 (df=5)
ML
98109.3
AIC
98119.3
disc
97600
> summary(gamm.11$mer);summary(gamm.11$gam)
Generalized linear mixed model fit by the Laplace approximation
AIC BIC logLik deviance
98119 98167 -49055 98109
Random effects:
Groups Name Variance Std.Dev.
vpnr (Intercept) 0.12965 0.36007
Xr.1 s(hours24) 1291.42444 35.93639
Number of obs: 97920, groups: vpnr, 114; Xr.1, 8
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
X(Intercept) 0.3713345 0.0644469 5.762 8.32e-09 ***
Xpred2 -0.0575848 0.0865231 -0.666 0.506
Xpred3 0.0003748 0.0869543 0.004 0.997
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.3713345 0.0149858 24.779 < 2e-16 ***
pred2 -0.0575848 0.0197488 -2.916 0.00355 **
pred3 0.0003748 0.0198803 0.019 0.98496
Output of 2nd model gamm.12:
-----------------------------------------
logLik(gamm.12$mer);deviance(gamm.12$mer);attributes(summary(gamm.12$mer))$AICtab[1];gamm.12$...@deviance["disc"]
'log Lik.' -48714.29 (df=8)
ML
97428.59
AIC
97444.59
disc
96840.12
> summary(gamm.12$mer);summary(gamm.12$gam)
Generalized linear mixed model fit by the Laplace approximation
AIC BIC logLik deviance
97445 97521 -48714 97429
Random effects:
Groups Name Variance Std.Dev.
vpnr (Intercept) 0.12777 0.35744
Xr.4 s(hours24):pred3 262.64433 16.20631
Xr.3 s(hours24):pred2 0.00000 0.00000
Xr.2 s(hours24):pred1 422.37197 20.55169
Xr.1 s(hours24) 1377.63864 37.11655
Number of obs: 97920, groups: vpnr, 114; Xr.4, 8; Xr.3, 8; Xr.2, 8; Xr.1, 8
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
X(Intercept) 0.33991 0.06402 5.310 1.10e-07 ***
Xpred2 -0.04912 0.08607 -0.571 0.568
Xpred3 0.11460 0.08691 1.319 0.187
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.33991 0.01504 22.598 < 2e-16 ***
pred2 -0.04912 0.02036 -2.412 0.0159 *
pred3 0.11460 0.02217 5.170 2.34e-07 ***
Approximate significance of smooth terms:
edf Ref.df Chi.sq p-value
s(hours24) 7.943e+00 7.943e+00 9030.4 <2e-16 ***
s(hours24):pred1 7.598e+00 7.598e+00 139.1 <2e-16 ***
s(hours24):pred2 5.000e-12 5.000e-12 0.0 NA
s(hours24):pred3 7.367e+00 7.367e+00 350.7 <2e-16 ***
[[alternative HTML version deleted]]
______________________________________________
[email protected] 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.
