date()
## [1] "Fri Apr 12 10:43:00 2013"
Krige the Parana rainfall.
require(geoR)
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"
plot(parana)
plot(parana, trend = "1st")
plot(parana, trend = "2nd")
The maximum distance is: 619.5
plot(variog4(parana, trend = "1st", max.dist = 310), omni = TRUE)
## 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(variog4(parana, trend = "2nd", max.dist = 310), omni = TRUE)
## 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 = "1st", max.dist = 310))
## variog: computing omnidirectional variogram
plot(variog(parana, trend = "2nd", max.dist = 310))
## variog: computing omnidirectional variogram
Find the best model:
iv = c(600, 160)
summary(likfit(parana, ini = iv, cov.model = "exp", fix.nug = TRUE, message = FALSE))$likelihood$AIC
## [1] 1394
summary(likfit(parana, ini = iv, cov.model = "exp", trend = "1st", fix.nug = FALSE,
message = FALSE))$likelihood$AIC
## [1] 1340
summary(likfit(parana, ini = iv, cov.model = "exp", trend = "2nd", fix.nug = TRUE,
message = FALSE))$likelihood$AIC
## [1] 1347
summary(likfit(parana, ini = iv, cov.model = "exp", trend = "2nd", fix.nug = FALSE,
message = FALSE))$likelihood$AIC
## [1] 1338
summary(likfit(parana, ini = iv, cov.model = "sph", trend = "2nd", fix.nug = TRUE,
message = FALSE))$likelihood$AIC
## [1] 1350
summary(likfit(parana, ini = iv, cov.model = "sph", trend = "2nd", fix.nug = FALSE,
message = FALSE))$likelihood$AIC
## [1] 1337
summary(likfit(parana, ini = iv, cov.model = "exp", trend = "2nd", fix.nug = FALSE))
## 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 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 = iv, fix.nugget = FALSE,
## cov.model = "exp")
summary(likfit(parana, ini = iv, cov.model = "sph", trend = "2nd", fix.nug = FALSE))
## 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 of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta0 beta1 beta2 beta3 beta4 beta5
## 425.075 0.075 -0.637 0.000 0.000 0.001
##
## Parameters of the spatial component:
## correlation function: spherical
## (estimated) variance parameter sigmasq (partial sill) = 329
## (estimated) cor. fct. parameter phi (range parameter) = 160
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 403
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 160
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-660" "9" "1337" "1364"
##
## non spatial model:
## log.L n.params AIC BIC
## "-671" "7" "1356" "1376"
##
## Call:
## likfit(geodata = parana, trend = "2nd", ini.cov.pars = iv, fix.nugget = FALSE,
## cov.model = "sph")
Plotting the models to see how they fit:
plot(variog(parana, trend = "2nd", max.dist = 310))
## variog: computing omnidirectional variogram
lines.variomodel(cov.model = "exp", cov.pars = c(372.6, 77.54), nug = 381.2,
max.dist = 310, col = "red")
lines.variomodel(cov.model = "sph", cov.pars = c(329.2, 160), nug = 403.3, max.dist = 310,
col = "green")
Save the two best (lowest AIC) models:
modelE = likfit(parana, ini = iv, cov.model = "exp", trend = "2nd", fix.nug = FALSE,
messages = TRUE)
## 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.
modelS = likfit(parana, ini = iv, cov.model = "sph", trend = "2nd", fix.nug = FALSE,
messages = TRUE)
## 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.
pred.grid = expand.grid(seq(137, 799, l = 200), seq(46, 508, l = 200))
parana.ks = krige.conv(parana, loc = pred.grid, krige = krige.control(trend.d = "2nd",
trend.l = "2nd", obj.m = modelS))
## 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
parana.ks2 = krige.conv(parana, loc = pred.grid, krige = krige.control(trend.d = "2nd",
trend.l = "2nd", obj.m = modelE))
## 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.ks, col = rev(heat.colors(20)), xlab = "Longitude", ylab = "Latitude")
title(main = "Spherical Model w/ Nonzero Nugget")
contour(parana.ks, add = TRUE, nlevel = 20)
image(parana.ks2, col = rev(heat.colors(20)), xlab = "Longitude", ylab = "Latitude")
title(main = "Exponential Model w/ Non-Zero Nugget")
contour(parana.ks2, add = TRUE, nlevel = 20)
They look very similar.
xv.wk = xvalid(parana, model = modelS)
## 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.wk$data, xv.wk$predicted, xlab = "Observed Rainfall", ylab = "Predicted Rainfall at Observed Location (spherical)",
asp = 1, col = "red", pch = 20)
abline(0, 1)
abline(lm(xv.wk$predicted ~ xv.wk$data), col = "red")
mean(xv.wk$error^2)
## [1] 534.3
mean(abs(xv.wk$error))
## [1] 17.89
xv.wk = xvalid(parana, model = modelE)
## 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.wk$data, xv.wk$predicted, xlab = "Observed Rainfall", ylab = "Predicted Rainfall at Observed Location (exponential)",
asp = 1, col = "red", pch = 20)
abline(0, 1)
abline(lm(xv.wk$predicted ~ xv.wk$data), col = "red")
mean(xv.wk$error^2)
## [1] 536.8
mean(abs(xv.wk$error))
## [1] 17.9
The mean squared cross validated prediction error is 534 units squared and the mean absolute cross validated prediction error is 18 units.