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$mer@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$mer@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]]

