Just starting to work on the variogram analysis of the Cikapundung dataset using geoR package.
My code outline will be:
I'll post the explanation later on.
# ANALYSIS FOR 1997 DATA
# LOADING LIBRARY
require("geoR")
## Loading required package: geoR
## Loading required package: sp
## Loading required package: MASS
## --------------------------------------------------------------
## Analysis of geostatistical data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.7-4 (built on 2012-06-29) is now loaded
## --------------------------------------------------------------
#require("gstat")
#require("lattice")
#require("knitcitations")
# --------
# LOADING DATA
Data97 <- read.csv("data97utm.csv")
# ANALYSIS FOR 1ST GROUP PARAMETER: ELEVATION, EC, pH, HARD, TDS, TEMP
## Pairs analysis
### Insert Pairs Code
## Variogram analysis
#### ??? What is variogram? It's basically similar to autocorrelation
#### but in spatial arrangement.
## Binding necessary columns
Data97.1 <- Data97[,c("x","y","elv","ec","ph",
"hard","tds","temp")]
### Assigning EC data as geodata (column 8)
EC97 <- as.geodata(Data97.1,coords.col=3:4,data.col=8) #variogram for EC column
### Checking any duplicate coordinates
dup.coords(EC97)
## NULL
### Building and plotting variogram
Var1.EC97 <- variog(EC97,trend="1st")
## variog: computing omnidirectional variogram
Var2.EC97 <- variog(EC97,trend="2nd")
## variog: computing omnidirectional variogram
#### Plot side by side
plot.new()
par(mfrow=c(2,1))
plot(Var1.EC97,pch=19,col="blue",main="1st order variogram")
plot(Var2.EC97,pch=19,col="red",main="2nd order variogram")
#### Plot overlaying
plot.new()
par(mfrow=c(1,1))
plot(Var1.EC97,pch=19,col="blue",main="Variogram EC Data 1997")
par(new=TRUE)
plot(Var2.EC97,pch=19,col="red",xaxt="n",yaxt="n")
##### ??? Question: How to choose the 1st or 2nd order?
### Fitting 1st EC variogram
init.values1 <- expand.grid(seq(0,150,10), seq(0,150,10)) #from geoR example
init.values2 <- expand.grid(seq(0,1.5e+07,by=100000),seq(0,1000, by=100)) #from Farzina
olsFit11.EC97 <- variofit(Var1.EC97,ini=init.values1,fix.nug=T,wei="equal")
## variofit: covariance model used is matern
## variofit: weights used: equal
## variofit: minimisation function used: optim
## variofit: searching for best initial value ... selected values:
## sigmasq phi tausq kappa
## initial.value "0" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 2.75374439558799
## Warning: unreasonable initial value for sigmasq + nugget (too low)
olsFit12.EC97 <- variofit(Var1.EC97,ini=init.values2,fix.nug=T,wei="equal")
## variofit: covariance model used is matern
## variofit: weights used: equal
## variofit: minimisation function used: optim
## variofit: searching for best initial value ... selected values:
## sigmasq phi tausq kappa
## initial.value "0" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 2.75374439558799
## Warning: unreasonable initial value for sigmasq + nugget (too low)
#### ??? Warning on both Fits: unreasonable initial value for sigmasq + nugget (too low), why?
summary(olsFit11.EC97)
## $pmethod
## [1] "OLS (ordinary least squares)"
##
## $cov.model
## [1] "matern"
##
## $spatial.component
## sigmasq phi
## 0.3953 0.0000
##
## $spatial.component.extra
## kappa
## 0.5
##
## $nugget.component
## tausq
## 0
##
## $fix.nugget
## [1] TRUE
##
## $fix.kappa
## [1] TRUE
##
## $practicalRange
## [1] 0.000116
##
## $sum.of.squares
## value
## 0.7224
##
## $estimated.pars
## sigmasq phi
## 0.3953 0.0000
##
## $weights
## [1] "equal"
##
## $call
## variofit(vario = Var1.EC97, ini.cov.pars = init.values1, fix.nugget = T,
## weights = "equal")
##
## attr(,"class")
## [1] "summary.variomodel"
summary(olsFit12.EC97)
## $pmethod
## [1] "OLS (ordinary least squares)"
##
## $cov.model
## [1] "matern"
##
## $spatial.component
## sigmasq phi
## 0.3953 0.0000
##
## $spatial.component.extra
## kappa
## 0.5
##
## $nugget.component
## tausq
## 0
##
## $fix.nugget
## [1] TRUE
##
## $fix.kappa
## [1] TRUE
##
## $practicalRange
## [1] 0.000116
##
## $sum.of.squares
## value
## 0.7224
##
## $estimated.pars
## sigmasq phi
## 0.3953 0.0000
##
## $weights
## [1] "equal"
##
## $call
## variofit(vario = Var1.EC97, ini.cov.pars = init.values2, fix.nugget = T,
## weights = "equal")
##
## attr(,"class")
## [1] "summary.variomodel"
wlsFit11.EC97 <- variofit(Var1.EC97,ini=init.values1,fix.nug=T,wei="npairs")
## variofit: covariance model used is matern
## variofit: weights used: npairs
## variofit: minimisation function used: optim
## variofit: searching for best initial value ... selected values:
## sigmasq phi tausq kappa
## initial.value "0" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 337.932288121874
## Warning: unreasonable initial value for sigmasq + nugget (too low)
wlsFit12.EC97 <- variofit(Var1.EC97,ini=init.values2,fix.nug=T,wei="npairs")
## variofit: covariance model used is matern
## variofit: weights used: npairs
## variofit: minimisation function used: optim
## variofit: searching for best initial value ... selected values:
## sigmasq phi tausq kappa
## initial.value "0" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 337.932288121874
## Warning: unreasonable initial value for sigmasq + nugget (too low)
#### ??? Warning on both Fits: unreasonable initial value for sigmasq + nugget (too low), why?
#### ??? What's the different between equal and npairs
summary(wlsFit11.EC97)
## $pmethod
## [1] "WLS (weighted least squares)"
##
## $cov.model
## [1] "matern"
##
## $spatial.component
## sigmasq phi
## 0.4768 0.0000
##
## $spatial.component.extra
## kappa
## 0.5
##
## $nugget.component
## tausq
## 0
##
## $fix.nugget
## [1] TRUE
##
## $fix.kappa
## [1] TRUE
##
## $practicalRange
## [1] 0.000116
##
## $sum.of.squares
## value
## 48.05
##
## $estimated.pars
## sigmasq phi
## 0.4768 0.0000
##
## $weights
## [1] "npairs"
##
## $call
## variofit(vario = Var1.EC97, ini.cov.pars = init.values1, fix.nugget = T,
## weights = "npairs")
##
## attr(,"class")
## [1] "summary.variomodel"
summary(wlsFit12.EC97)
## $pmethod
## [1] "WLS (weighted least squares)"
##
## $cov.model
## [1] "matern"
##
## $spatial.component
## sigmasq phi
## 0.4768 0.0000
##
## $spatial.component.extra
## kappa
## 0.5
##
## $nugget.component
## tausq
## 0
##
## $fix.nugget
## [1] TRUE
##
## $fix.kappa
## [1] TRUE
##
## $practicalRange
## [1] 0.000116
##
## $sum.of.squares
## value
## 48.05
##
## $estimated.pars
## sigmasq phi
## 0.4768 0.0000
##
## $weights
## [1] "npairs"
##
## $call
## variofit(vario = Var1.EC97, ini.cov.pars = init.values2, fix.nugget = T,
## weights = "npairs")
##
## attr(,"class")
## [1] "summary.variomodel"
### Fitting 2nd EC variogram
init.values21 <- expand.grid(seq(0,150,10), seq(0,150,10)) #from geoR example
init.values22 <- expand.grid(seq(0,1.5e+07,by=100000),seq(0,1000, by=100)) #from Farzina
olsFit21.EC97 <- variofit(Var2.EC97,ini=init.values21,fix.nug=T,wei="equal")
## variofit: covariance model used is matern
## variofit: weights used: equal
## variofit: minimisation function used: optim
## variofit: searching for best initial value ... selected values:
## sigmasq phi tausq kappa
## initial.value "0" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 2.28545755614521
## Warning: unreasonable initial value for sigmasq + nugget (too low)
olsFit22.EC97 <- variofit(Var2.EC97,ini=init.values22,fix.nug=T,wei="equal")
## variofit: covariance model used is matern
## variofit: weights used: equal
## variofit: minimisation function used: optim
## variofit: searching for best initial value ... selected values:
## sigmasq phi tausq kappa
## initial.value "0" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 2.28545755614521
## Warning: unreasonable initial value for sigmasq + nugget (too low)
#### ??? Warning on both Fits: unreasonable initial value for sigmasq + nugget (too low), why?
summary(olsFit21.EC97)
## $pmethod
## [1] "OLS (ordinary least squares)"
##
## $cov.model
## [1] "matern"
##
## $spatial.component
## sigmasq phi
## 0.377 0.000
##
## $spatial.component.extra
## kappa
## 0.5
##
## $nugget.component
## tausq
## 0
##
## $fix.nugget
## [1] TRUE
##
## $fix.kappa
## [1] TRUE
##
## $practicalRange
## [1] 0.000116
##
## $sum.of.squares
## value
## 0.4374
##
## $estimated.pars
## sigmasq phi
## 0.377 0.000
##
## $weights
## [1] "equal"
##
## $call
## variofit(vario = Var2.EC97, ini.cov.pars = init.values21, fix.nugget = T,
## weights = "equal")
##
## attr(,"class")
## [1] "summary.variomodel"
summary(olsFit22.EC97)
## $pmethod
## [1] "OLS (ordinary least squares)"
##
## $cov.model
## [1] "matern"
##
## $spatial.component
## sigmasq phi
## 0.377 0.000
##
## $spatial.component.extra
## kappa
## 0.5
##
## $nugget.component
## tausq
## 0
##
## $fix.nugget
## [1] TRUE
##
## $fix.kappa
## [1] TRUE
##
## $practicalRange
## [1] 0.000116
##
## $sum.of.squares
## value
## 0.4374
##
## $estimated.pars
## sigmasq phi
## 0.377 0.000
##
## $weights
## [1] "equal"
##
## $call
## variofit(vario = Var2.EC97, ini.cov.pars = init.values22, fix.nugget = T,
## weights = "equal")
##
## attr(,"class")
## [1] "summary.variomodel"
wlsFit21.EC97 <- variofit(Var2.EC97,ini=init.values21,fix.nug=T,wei="npairs")
## variofit: covariance model used is matern
## variofit: weights used: npairs
## variofit: minimisation function used: optim
## variofit: searching for best initial value ... selected values:
## sigmasq phi tausq kappa
## initial.value "0" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 275.686474473004
## Warning: unreasonable initial value for sigmasq + nugget (too low)
wlsFit22.EC97 <- variofit(Var2.EC97,ini=init.values22,fix.nug=T,wei="npairs")
## variofit: covariance model used is matern
## variofit: weights used: npairs
## variofit: minimisation function used: optim
## variofit: searching for best initial value ... selected values:
## sigmasq phi tausq kappa
## initial.value "0" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 275.686474473004
## Warning: unreasonable initial value for sigmasq + nugget (too low)
#### ??? Warning on both Fits: unreasonable initial value for sigmasq + nugget (too low), why?
#### ??? What's the different between equal and npairs
summary(wlsFit21.EC97)
## $pmethod
## [1] "WLS (weighted least squares)"
##
## $cov.model
## [1] "matern"
##
## $spatial.component
## sigmasq phi
## 0.4412 0.0000
##
## $spatial.component.extra
## kappa
## 0.5
##
## $nugget.component
## tausq
## 0
##
## $fix.nugget
## [1] TRUE
##
## $fix.kappa
## [1] TRUE
##
## $practicalRange
## [1] 0.000116
##
## $sum.of.squares
## value
## 27.46
##
## $estimated.pars
## sigmasq phi
## 0.4412 0.0000
##
## $weights
## [1] "npairs"
##
## $call
## variofit(vario = Var2.EC97, ini.cov.pars = init.values21, fix.nugget = T,
## weights = "npairs")
##
## attr(,"class")
## [1] "summary.variomodel"
summary(wlsFit22.EC97)
## $pmethod
## [1] "WLS (weighted least squares)"
##
## $cov.model
## [1] "matern"
##
## $spatial.component
## sigmasq phi
## 0.4412 0.0000
##
## $spatial.component.extra
## kappa
## 0.5
##
## $nugget.component
## tausq
## 0
##
## $fix.nugget
## [1] TRUE
##
## $fix.kappa
## [1] TRUE
##
## $practicalRange
## [1] 0.000116
##
## $sum.of.squares
## value
## 27.46
##
## $estimated.pars
## sigmasq phi
## 0.4412 0.0000
##
## $weights
## [1] "npairs"
##
## $call
## variofit(vario = Var2.EC97, ini.cov.pars = init.values22, fix.nugget = T,
## weights = "npairs")
##
## attr(,"class")
## [1] "summary.variomodel"
plot.new()
par(mfrow=c(1,2))
plot(Var1.EC97, main="1st order EC Semivar",pch=19,col="blue")
lines(wlsFit11.EC97,lty=2,col="green")
lines(olsFit11.EC97,lty=3,col="blue")
lines(wlsFit12.EC97,lty=2,col="grey")
lines(olsFit12.EC97,lty=3,col="red")
plot(Var2.EC97, main="2nd order EC Semivar",pch=19,col="red")
lines(wlsFit21.EC97,lty=2,col="green")
lines(olsFit21.EC97,lty=3,col="blue")
lines(wlsFit22.EC97,lty=2,col="grey")
lines(olsFit22.EC97,lty=3,col="red")
#legend("topleft",c("ordinary least squares","weighted least squares"),lty=c(2,1), lwd=c(1,1), col=c("blue","black"))
####### REFERENCES
#R Core Team (2014). R: A language and environment for statistical computing. R Foundation for
#Statistical Computing, Vienna, Austria. URL http://www.R-project.org/.
#RIBEIRO JR., P.J. and DIGGLE, P.J. (2001) geoR: A package for geostatistical analysis. R-NEWS Vol 1, No 2. ISSN 1609-3631.