31 July 2014

Looking at the predictors to water quality: using GAM



Dear friends,

Continuing my analysis on water quality, on this post I try to look at the predictors that control water quality. I am still using GAM with mgcv package. We can use (x,y) coordinate as one of the predictor, as tensor function using "te()" command. Character-type column can also be used as the predictor without "s()" command. In the following code, I am putting the:
  • EC (electro-conductivity) and concentrations of CO3, HCO3, CO2, Cl, SO4, NO2, NO3, Fe (in mg/L) as predictant
  • Spatial distribution "te(x,y)" + geology "(data$aq) column" and + elevation "s(elevation)" as predictors
One of the result from "plot(gam.., pages=1)" command will look something like the above image. I will elaborate more on the results and the workaround later on. The following is the sample code.


 

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
#---
# title : MGCV package tryout
# author: Dasapta Erwin Irawan^1 and Farzina Akter^2
# affiliation^1: Institut Teknologi Bandung (Indonesia)
# affiliation^2: University of Sydney (Australia)
# date  : 01 Aug 2014
#---

##########################
##### GAM ANALYSIS #######
##########################

# load library and data
require("mgcv")
data <- read.csv("alldata23.csv")
group1 <- data[,c("x", "y", "ec", "elv", "aq", 
                  "ph", "hard", "tds", "temp", 
                  "eh", "Q")]
group2 <- data[,c("x", "y", "ec", "elv", "aq", "Ca", "Mg", 
                  "Fe", "Mn", "K", "Na")]
group3 = data[,c("x", "y", "ec", "elv", "aq", "CO3","HCO3",
                 "CO2","Cl","SO4","NO2",
                 "NO3","SiO2")]


################## FAMILY = GAUSSIAN #####################

# smoothing=cubic shrinkage version ########################
k1<-10
bsm<-"cs"
Gcr10 <- gam(ec ~ te(x, y, k=k1, bs=bsm) +
            s(elv, k=k1, bs=bsm) + 
            (data$aq),
            data=group1)       
plot(Gcr10, pages=1)
gam.check(Gcr10)
summary(Gcr10)

Gcr11 <- gam(CO3 ~ te(x, y, k=k1, bs=bsm) +
            s(elv, k=k1, bs=bsm) + 
            (data$aq),
            data=group3) 
plot(Gcr11, pages=1)
gam.check(Gcr11)               
summary(Gcr11)

Gcr12 <- gam(HCO3 ~ te(x, y, k=k1, bs=bsm) +
            s(elv, k=k1, bs=bsm) + 
            (data$aq),
            data=group3) 
plot(Gcr12, pages=1)
gam.check(Gcr12)               
summary(Gcr12)

Gcr13 <- gam(CO2 ~ te(x, y, k=k1, bs=bsm) +
            s(elv, k=k1, bs=bsm) + 
            (data$aq),
            data=group3) 
plot(Gcr13, pages=1)
gam.check(Gcr13)               
summary(Gcr13)

Gcr14 <- gam(Cl ~ te(x, y, k=k1, bs=bsm) +
            s(elv, k=k1, bs=bsm) + 
            (data$aq),
            data=group3) 
plot(Gcr14, pages=1)
gam.check(Gcr14)               
summary(Gcr14)

Gcr15 <- gam(SO4 ~ te(x, y, k=k1, bs=bsm) +
            s(elv, k=k1, bs=bsm) + 
            (data$aq),
            data=group3) 
plot(Gcr15, pages=1)
gam.check(Gcr15)               
summary(Gcr15)

Gcr16 <- gam(NO2 ~ te(x, y, k=k1, bs=bsm) +
            s(elv, k=k1, bs=bsm) + 
            (data$aq),
            data=group3) 
plot(Gcr16, pages=1)
gam.check(Gcr16)               
summary(Gcr16)

Gcr17 <- gam(NO3 ~ te(x, y, k=k1, bs=bsm) +
            s(elv, k=k1, bs=bsm) + 
            (data$aq),
            data=group3) 
plot(Gcr17, pages=1)
gam.check(Gcr17)               
summary(Gcr17)

Gcr18 <- gam(Fe ~ te(x, y, k=k1, bs=bsm) +
            s(elv, k=k1, bs=bsm) + 
            (data$aq),
            data=group2) 
plot(Gcr18, pages=1)
gam.check(Gcr18)               
summary(Gcr18)

AIC(Gcr10, Gcr11, Gcr12)
AIC(Gcr13, Gcr14, Gcr15)
AIC(Gcr16, Gcr17, Gcr18)
Posted by at
Labels: , , ,

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)