variog{geoR}: 用來計算empirical variogram的function
variog(geodata, coords = geodata$coords, data = geodata$data,
uvec = "default", breaks = "default", trend = "cte",
lambda = 1, option = c("bin", "cloud", "smooth"),
estimator.type = c("classical", "modulus"),
nugget.tolerance, max.dist, pairs.min = 2,
bin.cloud = FALSE, direction = "omnidirectional",
tolerance = pi/8, unit.angle = c("radians","degrees"),
angles = FALSE, messages, ...)
breaks: 可以指定variogram的bins estimator.type: 可以選擇用一般的moment估計(“classical”)，或是Hawkins and Cressie's提出的Robust估計方法(“modulus”)
variofit{geoR}: 是用來fit一個variogram模型的，也就是造出一個跟empirical variogram很接近的theoretical model
variofit(vario, ini.cov.pars, cov.model,
fix.nugget = FALSE, nugget = 0,
fix.kappa = TRUE, kappa = 0.5,
simul.number = NULL, max.dist = vario$max.dist,
weights, minimisation.function,
limits = pars.limits(), messages, ...)
nugget: 加入nugget effect max.dist: fix range effect kappa: 這個值的選定，會覺得模型的使用。value of the smoothness parameter. Regarded as a fixed values if fix.kappa = TRUE or as a initial value for the minimization algorithm if fix.kappa = FALSE. Only required if one of the following correlation functions is used: “matern”, “powered.exponential”, “cauchy” and “gneiting.matern”. Defaults to 0.5.
likfit{geoR}: 用likelihood-based方法來估計空間相關模型，只限定在Gaussian model
likfit(geodata, coords = geodata$coords, data = geodata$data,
trend = "cte", ini.cov.pars, fix.nugget = FALSE, nugget = 0,
fix.kappa = TRUE, kappa = 0.5, fix.lambda = TRUE, lambda = 1,
fix.psiA = TRUE, psiA = 0, fix.psiR = TRUE, psiR = 1,
cov.model, realisations, lik.method = "ML", components = TRUE,
nospatial = TRUE, limits = pars.limits(),
print.pars = FALSE, messages, ...)
指令與variofit差不多 lik.method: 可以選擇maximum likelihood estimate (“ML”)或是restricted maximum likelihood estimate (“REML”)
Example: 這邊會採用四種較常用的估計方法(由於這些程式的估計結果過於繁雜，僅列出最後比較的圖形)
library(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
## --------------------------------------------------------------
data(s100)
vario100 <- variog(s100, max.dist=1)
## variog: computing omnidirectional variogram
ols <- variofit(vario100, ini=c(0.5,0.5), fix.nug=TRUE, wei="equal")
## variofit: covariance model used is matern
## variofit: weights used: equal
## variofit: minimisation function used: optim
summary(ols)
## $pmethod
## [1] "OLS (ordinary least squares)"
##
## $cov.model
## [1] "matern"
##
## $spatial.component
## sigmasq phi
## 1.1070 0.4006
##
## $spatial.component.extra
## kappa
## 0.5
##
## $nugget.component
## tausq
## 0
##
## $fix.nugget
## [1] TRUE
##
## $fix.kappa
## [1] TRUE
##
## $practicalRange
## [1] 1.2
##
## $sum.of.squares
## value
## 0.1025
##
## $estimated.pars
## sigmasq phi
## 1.1070 0.4006
##
## $weights
## [1] "equal"
##
## $call
## variofit(vario = vario100, ini.cov.pars = c(0.5, 0.5), fix.nugget = TRUE,
## weights = "equal")
##
## attr(,"class")
## [1] "summary.variomodel"
wls <- variofit(vario100, ini=c(0.5,0.5), fix.nug=TRUE)
## variofit: covariance model used is matern
## variofit: weights used: npairs
## variofit: minimisation function used: optim
summary(wls)
## $pmethod
## [1] "WLS (weighted least squares)"
##
## $cov.model
## [1] "matern"
##
## $spatial.component
## sigmasq phi
## 1.1902 0.4726
##
## $spatial.component.extra
## kappa
## 0.5
##
## $nugget.component
## tausq
## 0
##
## $fix.nugget
## [1] TRUE
##
## $fix.kappa
## [1] TRUE
##
## $practicalRange
## [1] 1.416
##
## $sum.of.squares
## value
## 28.68
##
## $estimated.pars
## sigmasq phi
## 1.1902 0.4726
##
## $weights
## [1] "npairs"
##
## $call
## variofit(vario = vario100, ini.cov.pars = c(0.5, 0.5), fix.nugget = TRUE)
##
## attr(,"class")
## [1] "summary.variomodel"
ml <- likfit(s100, ini=c(0.5, 0.5), fix.nug = TRUE)
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optimize.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optimize.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
ml
## likfit: estimated model parameters:
## beta sigmasq phi
## "0.777" "0.752" "0.183"
## Practical Range with cor=0.05 for asymptotic range: 0.5473
##
## likfit: maximised log-likelihood = -83.6
summary(ml)
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta
## 0.777
##
## Parameters of the spatial component:
## correlation function: exponential
## (estimated) variance parameter sigmasq (partial sill) = 0.752
## (estimated) cor. fct. parameter phi (range parameter) = 0.183
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (fixed) nugget = 0
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 0.5473
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-83.6" "3" "173" "181"
##
## non spatial model:
## log.L n.params AIC BIC
## "-126" "2" "256" "261"
##
## Call:
## likfit(geodata = s100, ini.cov.pars = c(0.5, 0.5), fix.nugget = TRUE)
reml <- likfit(s100, ini=c(0.5, 0.5), fix.nug = TRUE, lik.met = "REML")
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optimize.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optimize.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
summary(reml)
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: restricted maximum likelihood
##
## Parameters of the mean component (trend):
## beta
## 0.748
##
## Parameters of the spatial component:
## correlation function: exponential
## (estimated) variance parameter sigmasq (partial sill) = 0.847
## (estimated) cor. fct. parameter phi (range parameter) = 0.21
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (fixed) nugget = 0
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 0.6296
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-81.5" "3" "169" "177"
##
## non spatial model:
## log.L n.params AIC BIC
## "-125" "2" "254" "259"
##
## Call:
## likfit(geodata = s100, ini.cov.pars = c(0.5, 0.5), fix.nugget = TRUE,
## lik.method = "REML")
par(pty="s")
plot(vario100,main="Estimations of variogram model: S100 data")
lines(ols,lwd=2, lty=1,col=1)
lines(wls,lwd=2,lty=2,col=2)
lines(ml, lwd=2,lty=3,col=3)
lines(reml, lwd=2,lty = 4,col=4)
legend("topleft",c("OLSE","WLSE","MLE","RMLE"),lty=1:4,col=1:4,lwd=2,bty="n")