date()
## [1] "Thu Apr 18 20:26:35 2013"
Krige the Parana rainfall.
require(geoR)
## Loading required package: geoR
## Warning: package 'geoR' was built under R version 2.15.3
## Loading required package: sp
## Warning: package 'sp' was built under R version 2.15.3
## 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("parana")
summary(parana)
## $n
## [1] 143
##
## $coords.summary
## east north
## min 150.1 70.36
## max 768.5 461.97
##
## $distances.summary
## min max
## 1.0 619.5
##
## $data.summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 163 234 270 274 318 414
##
## $borders.summary
## east north
## min 138.0 46.77
## max 798.6 507.93
##
## $others
## [1] "loci.paper"
##
## attr(,"class")
## [1] "summary.geodata"
# The largest distance between any two points is 619.5
plot(parana)
plot(parana, trend = "1st")
plot(parana, trend = "2nd")
plot(variog(parana, trend = "1st", max.dist = 309.75))
## variog: computing omnidirectional variogram
plot(variog4(parana, trend = "1st", max.dist = 309.75))
## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogram
plot(variog(parana, trend = "2nd", max.dist = 309.75))
## variog: computing omnidirectional variogram
plot(variog4(parana, trend = "2nd", max.dist = 309.75))
## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogram
plot(variog(parana, trend = "2nd", max.dist = 309.75))
## variog: computing omnidirectional variogram
lines.variomodel(cov.model = "sph", cov.pars = c(150, 10), nug = 500, max.dist = 309.75,
col = "red")
lines.variomodel(cov.model = "exp", cov.pars = c(150, 10), nug = 500, max.dist = 309.75,
col = "green")
parana_sph.vm = likfit(parana, trend = "2nd", ini = c(150, 10), nug = 500, cov.model = "sph")
## kappa not used for the spherical correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## 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(parana_sph.vm)
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta0 beta1 beta2 beta3 beta4 beta5
## 401.168 0.118 -0.468 -0.001 0.000 0.001
##
## Parameters of the spatial component:
## correlation function: spherical
## (estimated) variance parameter sigmasq (partial sill) = 589
## (estimated) cor. fct. parameter phi (range parameter) = 21.4
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 109
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 21.44
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-666" "9" "1351" "1377"
##
## non spatial model:
## log.L n.params AIC BIC
## "-671" "7" "1356" "1376"
##
## Call:
## likfit(geodata = parana, trend = "2nd", ini.cov.pars = c(150,
## 10), nugget = 500, cov.model = "sph")
lines.variomodel(cov.model = "sph", cov.pars = c(589.4, 21.43), nug = 108.9,
max.dist = 309.75, col = "yellow")
parana_exp.vm = likfit(parana, trend = "2nd", ini = c(150, 10), nug = 500, cov.model = "exp")
## kappa not used for the exponential correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## 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(parana_exp.vm)
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta0 beta1 beta2 beta3 beta4 beta5
## 423.928 0.062 -0.636 0.000 0.000 0.001
##
## Parameters of the spatial component:
## correlation function: exponential
## (estimated) variance parameter sigmasq (partial sill) = 373
## (estimated) cor. fct. parameter phi (range parameter) = 77.5
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 381
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 232.3
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-660" "9" "1338" "1365"
##
## non spatial model:
## log.L n.params AIC BIC
## "-671" "7" "1356" "1376"
##
## Call:
## likfit(geodata = parana, trend = "2nd", ini.cov.pars = c(150,
## 10), nugget = 500, cov.model = "exp")
lines.variomodel(cov.model = "exp", cov.pars = c(372.6, 77.54), nug = 381.2,
max.dist = 309.75, col = "blue")
Looking at the two models, the blue exponential model looks significantly better and will be kept.
pred.grid = expand.grid(seq(137, 799, l = 200), seq(46, 508, l = 200))
parana_exp.ks = krige.conv(parana, loc = pred.grid, krige = krige.control(trend.d = "2nd",
trend.l = "2nd", obj.m = parana_exp.vm))
## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with mean given by a 2nd order polynomial on the coordinates
## krige.conv: Kriging performed using global neighbourhood
image(parana_exp.ks, col = terrain.colors(21), xlab = "longitude", ylab = "latitude")
contour(parana_exp.ks, add = TRUE, nlev = 21)
title(main = "Interpolated surface for Exponential Model")
xv_exp.wk = xvalid(parana, model = parana_exp.vm)
## xvalid: number of data locations = 143
## xvalid: number of validation locations = 143
## xvalid: performing cross-validation at location ... 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, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,
## xvalid: end of cross-validation
plot(xv_exp.wk$data, xv_exp.wk$predicted, xlab = "Observed Rainfall", ylab = "Predicted Rainfall at Observed Location",
asp = 1, col = "red", pch = 20)
abline(0, 1)
abline(lm(xv_exp.wk$predicted ~ xv_exp.wk$data), col = "red")
mean(xv_exp.wk$error^2)
## [1] 536.8
mean(abs(xv_exp.wk$error))
## [1] 17.9
The mean squared cross validated prediction error is 536.85 unit2 and the mean absolute cross validated prediction error is 17.9 units.