################
#Cesar
#2018-3-23
#GAM model
################
#安装mgcv包，安装一次即可
install.packages("mgcv")
#载入
library(mgcv)
#设置空间，放到自己想放的文件夹中，原始数据也放到该文件夹下
setwd("/home/cesar/LXF/3/")
#读取数据
dat<-read.csv("data.csv",header = T)
# attach 数据表，便于打参数
attach(dat)
#初步看一下各变量之间的关系，可以忽略
pairs(dat[,2:ncol(dat)],lower.panel = NULL)
#开始建模，选择变量
modgam1<-gam(density~s(av_tem),data = dat)
#输出，下同
sink("density~av-tem.txt")
summary(modgam1)
sink()
modgam2<-gam(density~s(av_tem)+s(flow),data = dat)
sink("density~av-tem+flow.txt")
summary(modgam2)
sink()
modgam3<-gam(density~s(av_tem)+s(flow)+s(water_level),data = dat)
sink("density~av-tem+flow+waterlevel.txt")
summary(modgam3)
sink()
modgam4<-gam(density~s(av_tem)+s(flow)+s(water_level)+s(av_tran),data=dat)
sink("density~av-tem+flow+waterlevel+av_tran.txt")
summary(modgam4)
sink()
#计算AIC
sink("AIC.txt")
AIC(modgam1,modgam2,modgam3,modgam4)
sink()
#计算anova
sink("anova.txt")
anova.gam(modgam1,modgam2,modgam3,modgam4,test = "Chisq")
sink()
##输出图形
dev.off()
plot(modgam1, all.terms=F, residuals=T, pch=20)
par(mfrow=c(1,2))
plot(modgam2, all.terms=F, residuals=T, pch=20)
par(mfrow=c(2,2))
plot(modgam3, all.terms=F, residuals=T, pch=20)
par(mfrow=c(2,2))
plot(modgam4, all.terms=F, residuals=T, pch=20)
#解除attach
detach(dat)
data.csv数据示例
summary结果示例
AIC结果示例
anova结果示例
作图示例
