Changes between Version 9 and Version 10 of GAM with interaction

Timestamp:

Apr 8, 2011 11:24:12 AM (14 years ago)

Author:

celine

Comment:

--

GAM with interaction

@@ -61 +61 @@
 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)[[BR]]
 See in http://www.highstat.com/book2.htm  for the code. [[BR]]
+
+{{{
+#   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")
+########################################################
+}}}