back to first page..

First method

​http://tolstoy.newcastle.edu.au/R/e2/help/07/03/11485.html
​http://www.math.yorku.ca/Who/Faculty/Monette/S-news/2284.html

Test avec annee:as.factor(uga) il faut utiliser la library mgcv car ne fonctionne pas avec gam
Mêmes résultats avec annee:uga[[BR]]

library(mgcv)
gam1<-gam(d~annee:as.factor(uga)+s(elev_mean),data=ers[ers$ot_effectif>0,],family=Gamma(link=log))
summary(gam1)
Family: Gamma 
Link function: log 

Formula:
d ~ annee:as.factor(uga) + s(elev_mean)

Parametric coefficients:
                                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)                           53.643721   7.348534   7.300 3.72e-13 ***
annee:as.factor(uga)Adour             -0.026323   0.003675  -7.162 1.01e-12 ***
annee:as.factor(uga)ArtoisPicardie    -0.026567   0.003668  -7.243 5.64e-13 ***
annee:as.factor(uga)Bretagne          -0.026196   0.003682  -7.114 1.42e-12 ***
annee:as.factor(uga)Corse             -0.025447   0.003674  -6.926 5.32e-12 ***
annee:as.factor(uga)Garonne           -0.026514   0.003677  -7.212 7.06e-13 ***
annee:as.factor(uga)Loire             -0.026645   0.003676  -7.248 5.43e-13 ***
annee:as.factor(uga)Rhin              -0.026609   0.003695  -7.202 7.55e-13 ***
annee:as.factor(uga)RhoneMediterranee -0.025929   0.003677  -7.052 2.20e-12 ***
annee:as.factor(uga)SeineNormandie    -0.026052   0.003680  -7.080 1.81e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Approximate significance of smooth terms:
               edf Ref.df     F p-value    
s(elev_mean) 7.947  8.691 135.7  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

R-sq.(adj) =  0.206   Deviance explained = 34.3%
GCV score = 1.5366  Scale est. = 1.527     n = 2861

> anova(gam1)

Family: Gamma 
Link function: log 

Formula:
d ~ annee:as.factor(uga) + s(elev_mean)

Parametric Terms:
                     df     F p-value
annee:as.factor(uga)  9 75.24  <2e-16

Approximate significance of smooth terms:
               edf Ref.df     F p-value
s(elev_mean) 7.947  8.691 135.7  <2e-16

Autre method

See in Mixed Effect Modes and Extensions in Ecology with R (2009), A. F.Zuur et al., Statistics for Biology and Health, Springer chapter 15 (page367-368)
See in ​http://www.highstat.com/book2.htm for the code.

#   Mixed Effects Models and Extensions in Ecology with R (2009)
#    Zuur, Ieno, Walker, Saveliev and Smith.    Springer
#    This file was produced by Alain Zuur (highstat@highstat.com)
#    www.highstat.com

#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.




library(AED); data(Ythan)
library(lattice)
Birds <- as.vector(as.matrix(Ythan[,2:8]))
X <- as.vector(as.matrix(Ythan[,9:15]))
YX14 <- c(Birds,X)
Year14 <- rep(Ythan$Year,14)
N <- length(Ythan$Year)
ID14<-factor(rep(names(Ythan[,2:15]),each = N),
     levels=c("wheat",
    "barley","oats","cattle","sheep","pigs",
    "nitrate","Oystercatcher","Turnstone","Curlew",
    "BartailedGodwit","Redshank","Knot","Dunlin"))


xyplot(YX14~Year14|ID14,xlab="Time",ylab="Variable",
   layout=c(4,4),
   scales=list(alternating=T,x=list(relation = "same"),
               y = list(relation = "free")),
   panel=function(x,y){
     panel.xyplot(x,y,col=1)
     panel.grid(h=-1,v=2)
     panel.loess(x,y,col=1,span=0.5)
     I2<-is.na(y)
     panel.text(x[I2],min(y,na.rm=T),'|',cex=0.5)  })


##########  Section 15.3 Estimation of trends


Birds7 <- as.vector(as.matrix(Ythan[,2:8]))
BirdNames<-c("Oystercatcher","Turnstone",
   "Curlew","BartailedGodwit","Redshank","Knot",
   "Dunlin")
ID7<-factor(rep(BirdNames, each = N), levels = BirdNames)
Year7<-rep(Ythan$Year,7)




Oyst.01 <- as.numeric(ID7=="Oystercatcher")
Turn.01 <- as.numeric(ID7=="Turnstone")
Curl.01 <- as.numeric(ID7=="Curlew")
Bart.01 <- as.numeric(ID7=="BartailedGodwit")
Reds.01 <- as.numeric(ID7=="Redshank")
Knot.01 <- as.numeric(ID7=="Knot")
Dunl.01 <- as.numeric(ID7=="Dunlin")

#For older R versions (2.7.0 and before), use:
f7 <- formula(Birds7 ~ ID7 +
    s(Year7, by = Oyst.01, bs = "cr") +
    s(Year7, by = Turn.01, bs = "cr") +
    s(Year7, by = Curl.01, bs = "cr") +
    s(Year7, by = Bart.01, bs = "cr") +
    s(Year7, by = Reds.01, bs = "cr") +
    s(Year7, by = Knot.01, bs = "cr") +
    s(Year7, by = Dunl.01, bs = "cr"))
    
    
#For later R versions, use:
f7 <- formula(Birds7 ~ ID7 + s(Year7, by = ID7, bs = "cr"))



library(mgcv);library(nlme)
lmc<-lmeControl(niterEM=5000,msMaxIter=1000)

M0<-gamm(f7,
      control=lmc,method="REML",
      weights=varIdent(form=~1|ID7))

AIC(M0$lme)


p1<-plot(M0$lme,resid(.,type="n")~fitted(.),abline=0,col=1)
p2<-plot(M0$lme,resid(.,type="n")~Year7,abline=0,col=1,xlab="Year")
p3<-plot(M0$lme,resid(.,type="n")~fitted(.)|ID7,
         abline=0,col=1,par.strip.text=list(cex = 0.75))


print(p1, position = c(0,0,1,1), split=c(1,1,2,2), more=TRUE)
print(p2, position = c(0,0,1,1), split=c(2,1,2,2), more=TRUE)
print(p3, position = c(0,0,2,1), split=c(1,2,2,2), more=FALSE)


plot(M0$lme,resid(.,type="n")~Year7|ID7,abline=0,col=1,xlab="Year")
########################################################