Showing posts with label ubuntu. Show all posts

13 July 2014

#R: GAM exercise using mgcv package (the workflow)

mgcv1

Dear friends,

This post is related to:

Hydromad exercise

Variogram exercise

Pairs exercise

I am still finishing my example of GAM (Generalised Additive Model) using MGCV package (by MGCV package)(by Simon Wood). I've been applying some variations of parameter setting. Next step will be choosing which one is the best model. In this practice I collaborate my code with Farzina Akter (a PhD student from Univ. of Sydney).

Our work flow would be like this:

Dataset preparation
GAM fitting using Gaussian family.
GAM fitting using Gamma family
Choosing the best model using AIC and ANOVA
Predict GAM for new data

Model try outs

Our models try outs can be seen in this table.

GAM	Family	Smoothing	Code
gam1	Gaussian	no smoothing	`gam1<-gam(y ~ s(cov1) + s(cov2) + ... + s(cov-n),data=data.frame)`
gam2	Gaussian	thin plate (tp)	`gam2<-gam(y ~ s(cov1,bs="tp") + s(cov2,bs="tp") + ... + s(cov-n,bs="tp"),data=data.frame)`
gam3	Gaussian	thin shrinkage (ts)	`gam3<-gam(y ~ s(cov1,bs="ts") + s(cov2,bs="ts") + ... + s(cov-n,bs="ts"),data=data.frame)`
gam4	Gaussian	cubic reg spline (cr)	`gam4<-gam(y ~ s(cov1,bs="cr") + s(cov2,bs="cr") + ... + s(cov-n,bs="cr"),data=data.frame)`
gam5	Gaussian	cubic shrinkage spline (cs)	`gam5<-gam(y ~ s(cov1,bs="cs") + s(cov2,bs="cs") + ... + s(cov-n,bs="cs"),data=data.frame)`
gam6	Gaussian	cubic cyclic spline (cc)	`gam6<-gam(y ~ s(cov1,bs="cc") + s(cov2,bs="cc") + ... + s(cov-n,bs="cc"),data=data.frame)`
gam7	Gamma	no smoothing	`gam7<-gam(y ~ s(cov1) + s(cov2) + ... + s(cov-n),Gamma (link="log"),data=data.frame)`
gam8	Gamma	thin plate (tp)	`gam8<-gam(y ~ s(cov1,bs="tp") + s(cov2,bs="tp") + ... + s(cov-n,bs="tp"),Gamma (link="log"),data=data.frame)`
gam9	Gamma	thin shrinkage (ts)	`gam9<-gam(y ~ s(cov1,bs="ts") + s(cov2,bs="ts") + ... + s(cov-n,bs="ts"),Gamma (link="log"),data=data.frame)`
gam10	Gamma	cubic reg spline (cr)	`gam10<-gam(y ~ s(cov1,bs="cr") + s(cov2,bs="cr") + ... + s(cov-n,bs="cr"),Gamma (link="log"),data=data.frame)`
gam11	Gamma	cubic shrinkage spline (cs)	`gam11<-gam(y ~ s(cov1,bs="cs") + s(cov2,bs="cs") + ... + s(cov-n,bs="cs"),Gamma (link="log"),data=data.frame)`
gam12	Gamma	cubic cyclic spline (cc)	`gam1<-gam(EC ~ s(cov1,bs="cc") + s(cov2,bs="cc") + ... + s(cov-n,bs="cc"),data=data.frame)`

You can see the change on parameter bs and Gamma(link="log"). You can rename all y and the cov based on your model.

After each of the GAM model, you can add the following lines (eg for gam1)

summary(gam1) gam.check(gam1) plot(gam1,pages=1) AIC(gam1) plot(gam1,residuals=T,pages=1)

Model summary

Then after you run all of the models you can summarise all the results with the following lines.

Anova test for family=Gaussian, link=identity, default

anovaGaussian<-anova(gam1,gam2,gam3,gam4,gam5,gam6,test="Chisq")

Anova test for family=Gamma, link=log

anovaGamma<-anova(gam7,gam8,gam9,gam10,gam11,gam12, test="Chisq")

AIC test for family=Gaussian, link=identity, default

AICGaussian<-AIC(gam1,gam2,gam3,gam4,gam5,gam6,test="Chisq")

AIC test for family=Gamma, link=log

AICGamma<-anova(gam7,gam8,gam9,gam10,gam11,gam12, test="Chisq")

Print the anova and AIC

print(anovaGaussian) print(anovaGamma) print(AICGaussian) print(AICGamma)

Or if you use xtable package, you can write the following lines to produce html table (or LaTeX table)

print(xtable(AIC1)) print(xtable(AIC2))

Choosing the best model

Then choose the best model by looking at the smallest AIC and anova result as additional model evaluation.

Thanks for visiting my blog.

(the md and pdf files are available at onlinewaterbook.wordpress.com)

This post is related to:

Hydromad exercise

Variogram exercise

Pairs exercise

#---
# 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  : 8 July 2014
#---

# This code is following http://www3.nd.edu/~mclark19/learn/GAMS.pdf

# Load library and data
require("mgcv")
data <- read.csv("data97utm.csv")

##########################
# PAIRS ANALYSIS #########
##########################

# Data group 1: Physical parameters
group1 <- data[,c("x","y","ec","elv","aq","ph","hard","tds","temp","eh","Q")]
pairs(group1,labels=colnames(group1),
      main="Physical parameter", 
      pch=21,bg=c('red','green3','blue','yellow')
      [unclass(data$aq)],
      upper.panel=NULL)
legend(x=0.6,y=0.8,levels(data$aq), 
       pt.bg=c('red','green3','blue','yellow'), 
       pch=21,bty='n',ncol=2)

# Data group 2: Cations (unit = ppm)
group2 <- data[,c("x","y","ec","elv","Ca","Mg","Fe","Mn","K","Na","Ca")]
pairs(group2,labels=colnames(group2),main="Cations", 
      pch=21,bg=c('red','green3','blue','yellow')
      [unclass(data$aq)],
      upper.panel=NULL)
legend(x=0.6,y=0.8,legend=levels(data$aq), 
       pt.bg=c('red','green3','blue','yellow'), 
       pch=21, ncol=2, bty ='n')

## Data group 3: Anions (unit = ppm)
group3 = data[,c("x","y","ec","CO3","HCO3","CO2","Cl","SO4","NO2",
      "NO3","SiO2")]
pairs(group3,labels=colnames(group3),main="Anions", 
      pch=21,bg=c('red','green3','blue','yellow')
      [unclass(data$aq)],upper.panel=NULL)
#par(xpd='TRUE')
legend(x=0.6,y=0.8,legend=levels(data$aq), 
       pt.bg=c('red','green3','blue','yellow'), 
       pch=21, ncol=2, bty ='n')

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

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

# GAM models (check, all predictors must be numeric)

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

## ols (k=10 default changed to k=5, to avoid smoothing error)
k1<-5 # k=10 (default)
gam11<-gam(ec ~ s(x,k=k1) + s(y,k=k1) + s(elv,k=k1) + 
             s(ph,k=k1) + s(hard,k=k1) + 
             s(tds,k=k1) + s(temp,k=k1) + s(eh,k=k1) + 
             s(Q,k=k1), data=group1)
k2<-3
gam12<-gam(ec ~ s(x,k=k2) + s(y,k=k2) + s(Ca,k=k2) + 
             s(Mg,k=k2) + s(Fe,k=k2) + s(Mn,k=k2) + 
             s(K,k=k2) + s(Na,k=k2), data=group2)
k3<-5
gam13<-gam(ec ~ s(x,k=k3) + s(y,k=k3) + s(CO3,k=k3) + 
             s(HCO3,k=k3) + s(CO2,k=k3) + s(Cl,k=k3) + 
             s(SO4,k=k3) + s(NO2,k=k3) + s(NO3,k=k3) +
             + s(SiO2,k=k3), data=group3)

## smoothing=thin plate smoothing
#k1<-5
gam21<-gam(ec ~ s(x,k=k1,bs="tp") + s(y,k=k1,bs="tp") + s(elv,k=k1,bs="tp") + 
             s(ph,k=k1,bs="tp") + s(hard,k=k1,bs="tp") + 
             s(tds,k=k1,bs="tp") + s(temp,k=k1,bs="tp") + s(eh,k=k1,bs="tp") + 
             s(Q,k=k1,bs="tp"), data=group1)
#k2<-3
gam22<-gam(ec ~ s(x,k=k2,bs="tp") + s(y,k=k2,bs="tp") + s(Ca,k=k2,bs="tp") + 
             s(Mg,k=k2,bs="tp") + s(Fe,k=k2,bs="tp") + s(Mn,k=k2,bs="tp") + 
             s(K,k=k2,bs="tp") + s(Na,k=k2,bs="tp"), data=group2)
#k3<-5
gam23<-gam(ec ~ s(x,k=k3,bs="tp") + s(y,k=k3,bs="tp") + s(CO3,k=k3,bs="tp") + 
             s(HCO3,k=k3,bs="tp") + s(CO2,k=k3,bs="tp") + s(Cl,k=k3,bs="tp") + 
             s(SO4,k=k3,bs="tp") + s(NO2,k=k3,bs="tp") + s(NO3,k=k3,bs="tp") +
             + s(SiO2,k=k3,bs="tp"), data=group3)

## smoothing=thin shrinkage 
#k1<-5
bsm<-"ts"
gam31<-gam(ec ~ s(x,k=k1,bs=bsm) + s(y,k=k1,bs=bsm) + s(elv,k=k1,bs=bsm) + 
             s(ph,k=k1,bs=bsm) + s(hard,k=k1,bs=bsm) + 
             s(tds,k=k1,bs=bsm) + s(temp,k=k1,bs=bsm) + s(eh,k=k1,bs=bsm) + 
             s(Q,k=k1,bs=bsm), data=group1)

#k2<-3
bsm<-"ts"
gam32<-gam(ec ~ s(x,k=k2,bs=bsm) + s(y,k=k2,bs=bsm) + s(Ca,k=k2,bs=bsm) + 
             s(Mg,k=k2,bs=bsm) + s(Fe,k=k2,bs=bsm) + s(Mn,k=k2,bs=bsm) + 
             s(K,k=k2,bs=bsm) + s(Na,k=k2,bs=bsm), data=group2)

#k3<-5
bsm<-"ts"
gam33<-gam(ec ~ s(x,k=k3,bs=bsm) + s(y,k=k3,bs=bsm) + s(CO3,k=k3,bs=bsm) + 
             s(HCO3,k=k3,bs=bsm) + s(CO2,k=k3,bs=bsm) + s(Cl,k=k3,bs=bsm) + 
             s(SO4,k=k3,bs=bsm) + s(NO2,k=k3,bs=bsm) + s(NO3,k=k3,bs=bsm) +
             + s(SiO2,k=k3,bs=bsm), data=group3)

# smoothing=cubic regression spline
#k1<-5
bsm<-"cr"
gam41<-gam(ec ~ s(x,k=k1,bs=bsm) + s(y,k=k1,bs=bsm) + s(elv,k=k1,bs=bsm) + 
             s(ph,k=k1,bs=bsm) + s(hard,k=k1,bs=bsm) + 
             s(tds,k=k1,bs=bsm) + s(temp,k=k1,bs=bsm) + s(eh,k=k1,bs=bsm) + 
             s(Q,k=k1,bs=bsm), data=group1)

#k2<-3
gam42<-gam(ec ~ s(x,k=k2,bs=bsm) + s(y,k=k2,bs=bsm) + s(Ca,k=k2,bs=bsm) + 
             s(Mg,k=k2,bs=bsm) + s(Fe,k=k2,bs=bsm) + s(Mn,k=k2,bs=bsm) + 
             s(K,k=k2,bs=bsm) + s(Na,k=k2,bs=bsm), data=group2)

#k3<-5
gam43<-gam(ec ~ s(x,k=k3,bs=bsm) + s(y,k=k3,bs=bsm) + s(CO3,k=k3,bs=bsm) + 
             s(HCO3,k=k3,bs=bsm) + s(CO2,k=k3,bs=bsm) + s(Cl,k=k3,bs=bsm) + 
             s(SO4,k=k3,bs=bsm) + s(NO2,k=k3,bs=bsm) + s(NO3,k=k3,bs=bsm) +
             + s(SiO2,k=k3,bs=bsm), data=group3)

# smoothing=cubic shrinkage version
bsm<-"cs"
#k1<-5
gam51<-gam(ec ~ s(x,k=k1,bs=bsm) + s(y,k=k1,bs=bsm) + s(elv,k=k1,bs=bsm) + 
             s(ph,k=k1,bs=bsm) + s(hard,k=k1,bs=bsm) + 
             s(tds,k=k1,bs=bsm) + s(temp,k=k1,bs=bsm) + s(eh,k=k1,bs=bsm) + 
             s(Q,k=k1,bs=bsm), data=group1)

#k2<-3
gam52<-gam(ec ~ s(x,k=k2,bs=bsm) + s(y,k=k2,bs=bsm) + s(Ca,k=k2,bs=bsm) + 
             s(Mg,k=k2,bs=bsm) + s(Fe,k=k2,bs=bsm) + s(Mn,k=k2,bs=bsm) + 
             s(K,k=k2,bs=bsm) + s(Na,k=k2,bs=bsm), data=group2)

#k3<-5
gam53<-gam(ec ~ s(x,k=k3,bs=bsm) + s(y,k=k3,bs=bsm) + s(CO3,k=k3,bs=bsm) + 
             s(HCO3,k=k3,bs=bsm) + s(CO2,k=k3,bs=bsm) + s(Cl,k=k3,bs=bsm) + 
             s(SO4,k=k3,bs=bsm) + s(NO2,k=k3,bs=bsm) + s(NO3,k=k3,bs=bsm) +
             + s(SiO2,k=k3,bs=bsm), data=group3)

# smoothing=cyclic cubic regression spline
# k1<-5
bsm<-"cc"
gam61<-gam(ec ~ s(x,k=k1,bs=bsm) + s(y,k=k1,bs=bsm) + s(elv,k=k1,bs=bsm) + 
             s(ph,k=k1,bs=bsm) + s(hard,k=k1,bs=bsm) + 
             s(tds,k=k1,bs=bsm) + s(temp,k=k1,bs=bsm) + s(eh,k=k1,bs=bsm) + 
             s(Q,k=k1,bs=bsm), data=group1)

k2<-3
gam62<-gam(ec ~ s(x,k=k2,bs=bsm) + s(y,k=k2,bs=bsm) + s(Ca,k=k2,bs=bsm) + 
             s(Mg,k=k2,bs=bsm) + s(Fe,k=k2,bs=bsm) + s(Mn,k=k2,bs=bsm) + 
             s(K,k=k2,bs=bsm) + s(Na,k=k2,bs=bsm), data=group2)

#k3<-5
gam63<-gam(ec ~ s(x,k=k3,bs=bsm) + s(y,k=k3,bs=bsm) + s(CO3,k=k3,bs=bsm) + 
             s(HCO3,k=k3,bs=bsm) + s(CO2,k=k3,bs=bsm) + s(Cl,k=k3,bs=bsm) + 
             s(SO4,k=k3,bs=bsm) + s(NO2,k=k3,bs=bsm) + s(NO3,k=k3,bs=bsm) +
             + s(SiO2,k=k3,bs=bsm), data=group3)

# Dropping "cc" model, causing error, don't have cyclic pattern


################## FAMILY = GAMMA #####################
## link=log, default smoothing
#k1<-5 # k=10 (default)
gam71<-gam(ec ~ s(x,k=k1) + s(y,k=k1) + s(elv,k=k1) + 
             s(ph,k=k1) + s(hard,k=k1) + 
             s(tds,k=k1) + s(temp,k=k1) + s(eh,k=k1) + 
             s(Q,k=k1), Gamma (link="log"), data=group1)
#k2<-3
gam72<-gam(ec ~ s(x,k=k2) + s(y,k=k2) + s(Ca,k=k2) + 
             s(Mg,k=k2) + s(Fe,k=k2) + s(Mn,k=k2) + 
             s(K,k=k2) + s(Na,k=k2), Gamma (link="log"), data=group2)
#k3<-5
gam73<-gam(ec ~ s(x,k=k3) + s(y,k=k3) + s(CO3,k=k3) + 
             s(HCO3,k=k3) + s(CO2,k=k3) + s(Cl,k=k3) + 
             s(SO4,k=k3) + s(NO2,k=k3) + s(NO3,k=k3) +
             + s(SiO2,k=k3), Gamma (link="log"), data=group3)

## smoothing=thin plate smoothing
#k1<-5
gam81<-gam(ec ~ s(x,k=k1,bs="tp") + s(y,k=k1,bs="tp") + s(elv,k=k1,bs="tp") + 
             s(ph,k=k1,bs="tp") + s(hard,k=k1,bs="tp") + 
             s(tds,k=k1,bs="tp") + s(temp,k=k1,bs="tp") + s(eh,k=k1,bs="tp") + 
             s(Q,k=k1,bs="tp"), 
             Gamma (link="log"), data=group1)
#k2<-3
gam82<-gam(ec ~ s(x,k=k2,bs="tp") + s(y,k=k2,bs="tp") + s(Ca,k=k2,bs="tp") + 
             s(Mg,k=k2,bs="tp") + s(Fe,k=k2,bs="tp") + s(Mn,k=k2,bs="tp") + 
             s(K,k=k2,bs="tp") + s(Na,k=k2,bs="tp"), 
             Gamma (link="log"), data=group2)
#k3<-5
gam83<-gam(ec ~ s(x,k=k3,bs="tp") + s(y,k=k3,bs="tp") + s(CO3,k=k3,bs="tp") + 
             s(HCO3,k=k3,bs="tp") + s(CO2,k=k3,bs="tp") + s(Cl,k=k3,bs="tp") + 
             s(SO4,k=k3,bs="tp") + s(NO2,k=k3,bs="tp") + s(NO3,k=k3,bs="tp") +
             + s(SiO2,k=k3,bs="tp"), 
             Gamma (link="log"), data=group3)

## smoothing=thin shrinkage 
#k1<-5
bsm<-"ts"
gam91<-gam(ec ~ s(x,k=k1,bs=bsm) + s(y,k=k1,bs=bsm) + s(elv,k=k1,bs=bsm) + 
             s(ph,k=k1,bs=bsm) + s(hard,k=k1,bs=bsm) + 
             s(tds,k=k1,bs=bsm) + s(temp,k=k1,bs=bsm) + s(eh,k=k1,bs=bsm) + 
             s(Q,k=k1,bs=bsm), 
             Gamma (link="log"), data=group1)

#k2<-3
bsm<-"ts"
gam92<-gam(ec ~ s(x,k=k2,bs=bsm) + s(y,k=k2,bs=bsm) + s(Ca,k=k2,bs=bsm) + 
             s(Mg,k=k2,bs=bsm) + s(Fe,k=k2,bs=bsm) + s(Mn,k=k2,bs=bsm) + 
             s(K,k=k2,bs=bsm) + s(Na,k=k2,bs=bsm),
             Gamma (link="log"), data=group2)

#k3<-5
bsm<-"ts"
gam93<-gam(ec ~ s(x,k=k3,bs=bsm) + s(y,k=k3,bs=bsm) + s(CO3,k=k3,bs=bsm) + 
             s(HCO3,k=k3,bs=bsm) + s(CO2,k=k3,bs=bsm) + s(Cl,k=k3,bs=bsm) + 
             s(SO4,k=k3,bs=bsm) + s(NO2,k=k3,bs=bsm) + s(NO3,k=k3,bs=bsm) +
             + s(SiO2,k=k3,bs=bsm), 
             Gamma (link="log"), data=group3)

# Family=gaussian, smoothing=cubic regression spline
#k1<-5
bsm<-"cr"
gam101<-gam(ec ~ s(x,k=k1,bs=bsm) + s(y,k=k1,bs=bsm) + s(elv,k=k1,bs=bsm) + 
             s(ph,k=k1,bs=bsm) + s(hard,k=k1,bs=bsm) + 
             s(tds,k=k1,bs=bsm) + s(temp,k=k1,bs=bsm) + s(eh,k=k1,bs=bsm) + 
             s(Q,k=k1,bs=bsm), 
             Gamma (link="log"), data=group1)

#k2<-3
gam102<-gam(ec ~ s(x,k=k2,bs=bsm) + s(y,k=k2,bs=bsm) + s(Ca,k=k2,bs=bsm) + 
             s(Mg,k=k2,bs=bsm) + s(Fe,k=k2,bs=bsm) + s(Mn,k=k2,bs=bsm) + 
             s(K,k=k2,bs=bsm) + s(Na,k=k2,bs=bsm), 
             Gamma (link="log"), data=group2)

#k3<-5
gam103<-gam(ec ~ s(x,k=k3,bs=bsm) + s(y,k=k3,bs=bsm) + s(CO3,k=k3,bs=bsm) + 
             s(HCO3,k=k3,bs=bsm) + s(CO2,k=k3,bs=bsm) + s(Cl,k=k3,bs=bsm) + 
             s(SO4,k=k3,bs=bsm) + s(NO2,k=k3,bs=bsm) + s(NO3,k=k3,bs=bsm) +
             + s(SiO2,k=k3,bs=bsm), 
             Gamma (link="log"), data=group3)

# smoothing=cubic shrinkage version
bsm<-"cs"
#k1<-5
gam111<-gam(ec ~ s(x,k=k1,bs=bsm) + s(y,k=k1,bs=bsm) + s(elv,k=k1,bs=bsm) + 
             s(ph,k=k1,bs=bsm) + s(hard,k=k1,bs=bsm) + 
             s(tds,k=k1,bs=bsm) + s(temp,k=k1,bs=bsm) + s(eh,k=k1,bs=bsm) + 
             s(Q,k=k1,bs=bsm), 
             Gamma (link="log"), data=group1)

#k2<-3
gam112<-gam(ec ~ s(x,k=k2,bs=bsm) + s(y,k=k2,bs=bsm) + s(Ca,k=k2,bs=bsm) + 
             s(Mg,k=k2,bs=bsm) + s(Fe,k=k2,bs=bsm) + s(Mn,k=k2,bs=bsm) + 
             s(K,k=k2,bs=bsm) + s(Na,k=k2,bs=bsm), 
             Gamma (link="log"), data=group2)

#k3<-5
gam113<-gam(ec ~ s(x,k=k3,bs=bsm) + s(y,k=k3,bs=bsm) + s(CO3,k=k3,bs=bsm) + 
             s(HCO3,k=k3,bs=bsm) + s(CO2,k=k3,bs=bsm) + s(Cl,k=k3,bs=bsm) + 
             s(SO4,k=k3,bs=bsm) + s(NO2,k=k3,bs=bsm) + s(NO3,k=k3,bs=bsm) +
             + s(SiO2,k=k3,bs=bsm), 
             Gamma (link="log"), data=group3)

# smoothing=cyclic cubic regression spline
# k1<-5
bsm<-"cc"
gam121<-gam(ec ~ s(x,k=k1,bs=bsm) + s(y,k=k1,bs=bsm) + s(elv,k=k1,bs=bsm) + 
             s(ph,k=k1,bs=bsm) + s(hard,k=k1,bs=bsm) + 
             s(tds,k=k1,bs=bsm) + s(temp,k=k1,bs=bsm) + s(eh,k=k1,bs=bsm) + 
             s(Q,k=k1,bs=bsm), 
             Gamma (link="log"), data=group1)

k2<-3
gam122<-gam(ec ~ s(x,k=k2,bs=bsm) + s(y,k=k2,bs=bsm) + s(Ca,k=k2,bs=bsm) + 
             s(Mg,k=k2,bs=bsm) + s(Fe,k=k2,bs=bsm) + s(Mn,k=k2,bs=bsm) + 
             s(K,k=k2,bs=bsm) + s(Na,k=k2,bs=bsm), 
             Gamma (link="log"), data=group2)

#k3<-5
gam123<-gam(ec ~ s(x,k=k3,bs=bsm) + s(y,k=k3,bs=bsm) + s(CO3,k=k3,bs=bsm) + 
             s(HCO3,k=k3,bs=bsm) + s(CO2,k=k3,bs=bsm) + s(Cl,k=k3,bs=bsm) + 
             s(SO4,k=k3,bs=bsm) + s(NO2,k=k3,bs=bsm) + s(NO3,k=k3,bs=bsm) +
             + s(SiO2,k=k3,bs=bsm), 
             Gamma (link="log"), data=group3)

# Same situation as GAUSSIAN family: Dropping "cc" model, causing error, don't have cyclic pattern

######### GAM EVALUATION ################
# Gaussian
AIC.gsdef<-AIC(gam11,gam12,gam13)
AIC.gstp<-AIC(gam21,gam22,gam23)
AIC.gsts<-AIC(gam31,gam32,gam33)
AIC.gscr<-AIC(gam41,gam42,gam43)
AIC.gscs<-AIC(gam51,gam52,gam53)
#AIC.gscc<-AIC(gam61,gam62,gam63) # dropped
print(AIC.gsdef) ; print(AIC.gstp) # lowestAIC=gam13(572.5448) and gam23(572.5448)
print(AIC.gsts) ; print(AIC.gscr) # lowestAIC=gam33(571.6253) and gam43(587.7679)
print(AIC.gscs) # lowest AIC=gam53(593.4830)
#print(AIC.gscc) # dropped

summary(gam13) 
# R-sq=0.897, GCV=6765, scale=3120, sigpar=all,ytrend
gam.check(gam13)

summary(gam23)
# R-sq=0.897, GCV=6765, scale=3120, sigpar=all,ytrend,-HCO3
gam.check(gam23)

summary(gam33)
# R-sq=0.899, GCV=6676.2, scale=3065.4, sigpar=all,xytrend
gam.check(gam33)

summary(gam43)
# R-sq=0.861, GCV=7623.9, scale=4213.8, sigpar=all,xytrend,-Cl,NO2
gam.check(gam43)

summary(gam53)
# R-sq=0.845, GCV=8920, scale=4695.3, sigpar=all,ytrend,-Cl,NO2
gam.check(gam53)

# Gamma
# using AIC
AIC.gmdef<-AIC(gam71,gam72,gam73)
AIC.gmtp<-AIC(gam81,gam82,gam83)
AIC.gmts<-AIC(gam91,gam92,gam93)
AIC.gmcr<-AIC(gam101,gam102,gam103)
AIC.gmcs<-AIC(gam111,gam112,gam113)
#AIC.gmcc<-AIC(gam121,gam122,gam123) # dropped
print(AIC.gmdef) ; print(AIC.gmtp) # lowestAIC=gam73(603.1337) and gam83(603.1337)
print(AIC.gmts) ; print(AIC.gmcr) # lowestAIC=gam93(601.2405) and gam103(598.4275)
print(AIC.gmcs) # lowestAIC=gam111(597.2866)
#print(AIC.gmcc) # dropped

# using summary and gam.check
summary(gam73)
# R-sq=0.66, GCV=0.20824, scale=0.13229, sigpar=all,ytrend,(-HCO3,Cl,SO4)
gam.check(gam73)

summary(gam83)
# R-sq=0.66, GCV=0.20824, scale=0.13229, sigpar=all,ytrend,(-HCO3,Cl,SO4)
gam.check(gam83)

summary(gam93)
# R-sq=0.664, GCV=0.19493, scale=0.12899, sigpar=all,ytrend,(-HCO3,Cl,SO4)
gam.check(gam93)

summary(gam103)
# R-sq=0.617, GCV=0.20011, scale=0.11909, sigpar=all,ytrend,(-HCO3,Cl,SO4,NO2)
gam.check(gam103)

summary(gam111)
# R-sq=0.602, GCV=0.17701, scale=0.12065, sigpar=all,ytrend,(-elv,ph,hard)
gam.check(gam111)

17 June 2014

[#R] How to convert lat-long to UTM system

[#R] How to convert lat-long coordinates to UTM (easting-northing)

1. Intro (`rgdal` installation on Mac)

This bit is part of my work in modeling the hydrology of Cikapundung Catchment. I am kind of (forced) to do the spatial analysis in R :-). But after several trial and errors, it's really not that hard.

So If you do spatial analysis using R, you would need to set your coordinate into easting and northing system. Many R packages in spatial section only read this projection system (correct me if I'm wrong). One of the is geoR package. I'll talk about this later, since I'm still figuring out how to use it to make variograms and trend maps.

Now, to do a coordinate conversion, you would need rgdal package. This is the basic R package that deals with spatial analysis. Just install it with this command:

install.packages('rgdal')

In my case, the installation went well in Mac (running Snow Leopard) and also Windows 7, but it needed a workaround to install it on my Linux machine (Ubuntu 13.10). After installation, call the package using

library(rgdal) or require(rgdal)

2. Installing `rgdal` on Ubuntu

As I have said, installing rgdal on Linux was not as smooth as installing the other spatial packages, say maptools, geoR, or sp. At least in my case.

It has something to do with libgdal1-dev and libproj-dev files and their version which is corresponding to the Linux kernel version.

I googled it and this terminal command did work on my machine running Ubuntu 13.10. See this post:

sudo add-apt-repository ppa:marutter/c2d4u

sudo apt-get update

sudo apt-get install r-cran-rgdal

Another way (that did not work for me, but it might work on you). According to the same link above, this command was re-ran after switching off the related update site in the Ubuntu's update manager database:

sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable

sudo apt-get update

sudo apt-get install libproj-dev

Another way suggested by Robin Lovelace was working on Ubuntu 13.04.

3. Converting coordinates

Then load your data, in my case was

data = read.csv('27-1997data.csv') or

data <- read.csv('27-1997data.csv')

I often use the first one with “=” sign to avoid me hitting “shift key”.

Then you must set your existing coordinate system. Because otherwise, R will treat your x and y as plain data, not coordinate. In this case we are going to convert the existing lat-long coordinates to UTM (easting-northing system).

Run this command:

cord.dec = SpatialPoints(cbind(data$long, -data$lat), proj4string=CRS("+proj=longlat"))

You can change the “cord.dec”, “data$long”, and “data$lat”, based on your dataset.

Then you are going to transform the CRS or Coordinate Reference System to UTM by setting the EPSG to 32748 for WGS 84, UTM Zone 48M, southern hemisphere. This is where Bandung is located.

cord.UTM <- spTransform(cord.dec, CRS("+init=epsg:32748"))

cord.UTM

4. Coordinate conversion on GIS app

I know that this post was mainly about R. But It does not hurt to share a system using QGIS and GRASS, posted by a fellow user Dmitri, from GRASS GIS Facebook Community.

So basically we will have two options:

Option 1

Use QGIS to reproject your data from WGS-84 (EPSG:4326) to UTM 48 zone (EPSG:34728). To do that in QGIS: (a) run QGIS, (b) load the dataset, © then use menu: Layer->Save As… (select new proection in the dialog). Then you can create new GRASS mapset using your reprojected data.

Option 2

Reproject the data using GRASS: (a) Create mapset in WGS-84 (b) Import the data, © create mapset in UTM 48 zone, (d) reproj data from the first mapset to the second (r.proj and v.proj commands).

5. The complete R codes.

Do not expect fancy code with me (LOL).

# loading library

require(rgdal)

## Loading required package: rgdal
## Loading required package: sp
## rgdal: version: 0.8-16, (SVN revision 498)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgdal/proj

data = read.csv("27-1997data.csv")
cord.dec = SpatialPoints(cbind(data$long, -data$lat), proj4string = CRS("+proj=longlat"))


# Setting existing coordinate as lat-long system
cord.dec = SpatialPoints(cbind(data$long, -data$lat), proj4string = CRS("+proj=longlat"))

# Transforming coordinate to UTM using EPSG=32748 for WGS=84, UTM Zone=48M,
# Southern Hemisphere)
cord.UTM <- spTransform(cord.dec, CRS("+init=epsg:32748"))
cord.UTM

## SpatialPoints:
##       coords.x1 coords.x2
##  [1,]    783993  10767111
##  [2,]    778949  10758896
##  [3,]    784001  10765562
##  [4,]    789969  10766259
##  [5,]    783890  10765672
##  [6,]    785760  10767453
##  [7,]    777100  10753022
##  [8,]    776305  10757112
##  [9,]    776193  10757222
## [10,]    785220  10765016
## [11,]    779750  10753921
## [12,]    790006  10759509
## [13,]    785175  10752843
## [14,]    780604  10759679
## [15,]    784444  10765343
## [16,]    777882  10751587
## [17,]    782471  10762234
## [18,]    782032  10761457
## [19,]    783456  10764121
## [20,]    786240  10760484
## [21,]    775304  10758102
## [22,]    785175  10752843
## [23,]    785175  10752843
## [24,]    785175  10752843
## [25,]    785175  10752843
## [26,]    782141  10761790
## [27,]    791358  10754868
## [28,]    776638  10756892
## [29,]    779824  10760892
## [30,]    783354  10762571
## [31,]    783804  10761135
## [32,]    779808  10763769
## [33,]    785175  10752843
## [34,]    785175  10752843
## [35,]    785379  10756053
## [36,]    785175  10752843
## [37,]    778412  10755906
## [38,]    785310  10748307
## [39,]    784556  10765233
## [40,]    784116  10764678
## [41,]    785175  10752843
## [42,]    768459  10755854
## [43,]    785175  10752843
## [44,]    771336  10755537
## [45,]    785257  10758266
## [46,]    785175  10752843
## [47,]    781714  10759021
## [48,]    779287  10757791
## [49,]    785247  10760147
## [50,]    790201  10764158
## [51,]    784267  10757264
## Coordinate Reference System (CRS) arguments: +init=epsg:32748
## +proj=utm +zone=48 +south +datum=WGS84 +units=m +no_defs
## +ellps=WGS84 +towgs84=0,0,0


# Plotting points
par(mfrow = c(1, 2))
plot(cord.dec, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95)
plot(cord.UTM, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95)

plot of chunk unnamed-chunk-1

01 May 2014

#ubuntu: CA certificate request on Uni wifi network

Today I was using my laptop running Ubuntu 13.10 to the Uni. When it tried to get in to the "Uni Sydney" wifi network, it requested the usual userid and password. After filled them in, then it started to ask for CA Certificate. A little bit confusion at first, because I didn't get it with my Mac. I'm not quite sure why it was happening with Ubuntu.
After googling a little bit, a CA (Certificate Authority) is:

In cryptography, a certificate authority or certification authority (CA), is an entity that issues digital certificates. The digital certificate certifies the ownership of a public key by the named subject of the certificate. This allows others (relying parties) to rely upon signatures or assertions made by the private key that corresponds to the public key that is certified. In this model of trust relationships, a CA is a trusted third party that is trusted by both the subject (owner) of the certificate and the party relying upon the certificate. CAs are characteristic of many public key infrastructure (PKI) schemes (http://en.wikipedia.org/wiki/Certificate_authority).

It turned out that we can create our own CA on Ubuntu. Just visit Ubuntu Official Help Site on this link. With a little bit Terminal script, you will be able to use the network.

11 April 2014

[ubuntu] QGIS up and running

Finally: QGIS on up and running on my Saucy

After trying out several tutorial, finally I can setup the QGIS to run on my machine with the following commands.

Be sure to check that there are no http://qgis.org/debian lines from your /etc/apt/sources.list

If you aren't sure, just check open Terminal Window. Then type sudo gedit /etc/apt/sources.list. Then check the last commands for any qgis lines. Delete them, then hit save and close.

Next type the following commands on your terminat prompt. Hit return on every line.

sudo apt-get install python-software-properties
sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable
sudo apt-get update
sudo apt-get install qgis python-qgis qgis-plugin-grass

It should be working fine.

Good luck.

(thanks to: askubuntu)

05 April 2014

[ubuntu] Auto restart problem

(www.inventorairconditioner.com) [ubuntu] Auto restart problem after shutdown command

[ubuntu] Auto restart problem after shutdown command

I have experienced the problem since I installed my machine with Ubuntu 13.10. If I hit the shutdown button, my laptop (Acer Aspire V5) executed a shutdown sequence. But after several seconds, it went on again. Then I found this website owned by +Thierry Nicola.

So I did the following step, as mentioned in the site: Open a terminal window

$ gksudo gedit /etc/rc.local.shutdown

the command will open the file /etc/rc.local.shutdown

Then I type this command in the file:

rmmod ehci_hcd

Next save the file.

Then try to shutdown your machine. It works fine in mine.

But strangely on my other laptop (Acer Aspire V5) as well as other same laptops I found on Google Discussion behave differently.

It goes well by installing Laptop Mode Tools, as shown in the following code

sudo apt-get update && sudo apt-get install laptop-mode-tools

Don't know why, don't want to know, as long as it works perfectly :-)

Pages