Compute variogram for the midwest ozone field
good<- !is.na(ozone2$y[16,])
x<- ozone2$lon.lat[good,]
y<- ozone2$y[16,good]
## variogram plots
par(mfrow=c(1,2))
look<-vgram(x,y, N=15, lon.lat=TRUE)
plot(look, pch=19)
## or some boxplot bin summaries
brk<- seq(0, 250,, (25 + 1) ) # will give 25 bins.
boxplotVGram(look, breaks=brk, plot.args=list(type="o"))
plot(look, add=TRUE, breaks=brk, col=4)
Kriging
# a 2-d example
# fitting a surface to ozone
# measurements. Exponential covariance, range parameter is 20 (in miles)
fit <- Krig(ChicagoO3$x, ChicagoO3$y, theta=20)
summary(fit) # summary of fit
## CALL:
## Krig(x = ChicagoO3$x, Y = ChicagoO3$y, theta = 20)
##
## Number of Observations: 20
## Number of unique points: 20
## Number of parameters in the null space 3
## Parameters for fixed spatial drift 3
## Effective degrees of freedom: 5.4
## Residual degrees of freedom: 14.6
## MLE sigma 3.733
## GCV sigma 4.049
## MLE rho 4.796
## Scale passed for covariance (rho) <NA>
## Scale passed for nugget (sigma^2) <NA>
## Smoothing parameter lambda 2.906
##
## Residual Summary:
## min 1st Q median 3rd Q max
## -6.4820 -1.6170 -0.5672 1.8010 7.4870
##
## Covariance Model: stationary.cov
## Covariance function is
## Names of non-default covariance arguments:
## theta
##
## DETAILS ON SMOOTHING PARAMETER:
## Method used: REML Cost: 1
## lambda trA GCV GCV.one GCV.model shat
## 2.906 5.394 22.452 22.452 NA 4.049
##
## Summary of all estimates found for lambda
## lambda trA GCV shat -lnLike Prof converge
## GCV 4.120 4.799 22.39 4.125 49.26 4
## GCV.model NA NA NA NA NA NA
## GCV.one 4.120 4.799 22.39 4.125 NA 4
## RMSE NA NA NA NA NA NA
## pure error NA NA NA NA NA NA
## REML 2.906 5.394 22.45 4.049 49.24 2
par(mfrow=c(2,2))
plot(fit) # four diagnostic plots of fit
# predict at data
predict(fit)
## [,1]
## [1,] 38.32080
## [2,] 38.83593
## [3,] 38.12672
## [4,] 38.76760
## [5,] 39.04168
## [6,] 39.90186
## [7,] 38.56459
## [8,] 39.53299
## [9,] 39.89821
## [10,] 39.65490
## [11,] 40.75214
## [12,] 40.19162
## [13,] 40.92551
## [14,] 40.57972
## [15,] 41.04454
## [16,] 40.43037
## [17,] 39.62459
## [18,] 39.41161
## [19,] 41.03544
## [20,] 40.96869
# predict on a grid (grid chosen here by defaults)
out<- predictSurface(fit)
surface(out, type="C") # option "C" our favorite
# predict at arbitrary points (10,-10) and (20, 15)
xnew<- rbind(c(10, -10), c(20, 15))
predict(fit, xnew)
## [,1]
## [1,] 38.28184
## [2,] 40.13606
# standard errors of prediction based on covariance model.
predictSE(fit, xnew)
## [1] 1.632105 3.054689
# surface of standard errors on a default grid
predictSurfaceSE(fit)-> out.p # this takes some time!
surface(out.p, type="C")
points(fit$x)
Kriging with Matern covariance function
## Using another stationary covariance.
# smoothness is the shape parameter for the Matern
fit <- Krig(ChicagoO3$x, ChicagoO3$y,
Covariance="Matern", theta=10, smoothness=1.0)
summary( fit)
## CALL:
## Krig(x = ChicagoO3$x, Y = ChicagoO3$y, Covariance = "Matern",
## theta = 10, smoothness = 1)
##
## Number of Observations: 20
## Number of unique points: 20
## Number of parameters in the null space 3
## Parameters for fixed spatial drift 3
## Effective degrees of freedom: 4.8
## Residual degrees of freedom: 15.2
## MLE sigma 3.807
## GCV sigma 4.129
## MLE rho 3.39
## Scale passed for covariance (rho) <NA>
## Scale passed for nugget (sigma^2) <NA>
## Smoothing parameter lambda 4.276
##
## Residual Summary:
## min 1st Q median 3rd Q max
## -6.7380 -1.7120 -0.5869 1.9360 7.7350
##
## Covariance Model: stationary.cov
## Covariance function is Matern
## Names of non-default covariance arguments:
## Covariance, theta, smoothness
##
## DETAILS ON SMOOTHING PARAMETER:
## Method used: REML Cost: 1
## lambda trA GCV GCV.one GCV.model shat
## 4.276 4.849 22.505 22.505 NA 4.129
##
## Summary of all estimates found for lambda
## lambda trA GCV shat -lnLike Prof converge
## GCV 5.868 4.423 22.45 4.182 49.27 4
## GCV.model NA NA NA NA NA NA
## GCV.one 5.868 4.423 22.45 4.182 NA 4
## RMSE NA NA NA NA NA NA
## pure error NA NA NA NA NA NA
## REML 4.276 4.849 22.50 4.129 49.26 1
